Many-body effects in the excitations and dynamics of trapped Bose-Einstein condensates
MMany-body effects in the excitations and dynamics of trappedBose-Einstein condensates
Ofir E. Alon ∗ , Raphael Beinke † , and Lorenz S. Cederbaum ‡ Department of Mathematics, University of Haifa, Haifa 3498838, Israel Haifa Research Center for Theoretical Physics and Astrophysics, University of Haifa, Haifa3498838, Israel Theoretische Chemie, Physikalisch-Chemisches Institut, Universit¨at Heidelberg, ImNeuenheimer Feld 229, D-69120 Heidelberg, GermanyJanuary 29, 2021 ∗ ofi[email protected] † [email protected] ‡ [email protected] a r X i v : . [ c ond - m a t . qu a n t - g a s ] J a n bstract This review explores the dynamics and the low-energy excitation spectra of Bose-Einstein condensates(BECs) of interacting bosons in external potential traps putting particular emphasis on the emerging many-body effects beyond mean-field descriptions. To do so, methods have to be used that, in principle, canprovide numerically exact results for both the dynamics and the excitation spectra in a systematic manner.Numerically exact results for the dynamics are presented employing the well-established multicongurationaltime-dependent Hartree for bosons (MCTDHB) method. The respective excitation spectra are calculatedutilizing the more recently introduced linear-response theory atop it (LR-MCTDHB). The latter theorygives rise to an, in general, non-hermitian eigenvalue problem. The theory and its newly developedimplementation are described in detail and benchmarked towards the exactly-solvable harmonic-interactionmodel. Several applications to BECs in one- and two-dimensional potential traps are discussed. Withrespect to dynamics, it is shown that both the out-of-equilibrium tunneling dynamics and the dynamics oftrapped vortices are of many-body nature. Furthermore, many-body effects in the excitation spectra arepresented for BECs in different trap geometries. It is demonstrated that even for essentially-condensedsystems, the spectrum of the lowest-in-energy excitations computed at the many-body level can differsubstantially from the standard mean-field description. In general, it is shown that bosons carrying angularmomentum are more sensitive to many-body effects than bosons without. These effects are present in boththe dynamics and the excitation spectrum. ontents ppendix A Mean-field theory 79 A.1 Single-orbital approaches . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 80A.1.1 The GP equation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 80A.1.2 The BdG equations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 81A.2 Multi-orbital approaches . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 82A.2.1 Time-dependent multi-orbital mean field (TDMF) . . . . . . . . . . . . . . . . . . . . . 82A.2.2 Linear-response best-mean-field (LR-BMF) . . . . . . . . . . . . . . . . . . . . . . . . . 83
Appendix B LR-MCTDHB in block-diagonal form 85Appendix C Further benchmarks of LR-MCTDHB 88Appendix D Many-particle variance: A sensitive quantity to many-body effects 90 Introduction
In 1924, Albert Einstein discovered that the spectral energy distribution in certain systems of massive particlesis different from the Maxwell-Boltzmann statistics. He was inspired by a work of Satayendra N. Bose whoderived the thermal energy distribution of black-body radiation solely based on assumptions from quantumtheory [1], which is slightly different to the previous semi-classical approach of Max Planck. Einstein realizedthe far-reaching consequences of this, namely that a composite system of massive particles at sufficiently lowtemperatures could behave as a matter wave due to all particles sharing the same single-particle state. Thus,he applied Bose’s idea to a single-component ideal gas [2, 3]. This paved the way towards what is known asBose-Einstein condensation. Consequently, all particles that obey the statistical distribution found by Boseand Einstein are termed bosons.It took approximately 70 years until the first gaseous Bose-Einstein condensate (BEC) was realized in anexperiment. The group of Eric A. Cornell and Carl E. Wieman produced a repulsive BEC of rubidium atoms[4], shortly followed by the group of Wolfgang Ketterle using repulsive sodium atoms in a magneto-opticaltrap [5]. The latter three scientists obtained the Nobel prize in physics for this achievement in 2001. Thefirst attractive BEC was realized by the group of Randall G. Hulet using lithium atoms [6]. In more recentexperiments, condensation was observed in a large variety of systems, e.g., for exciton-polariton systems [7, 8],magnons at room temperature [9], or, rather astonishingly because of the lack of a rest mass and a vanishingchemical potential, even for photons in a cavity [10].Over the past two decades, the experimental control of BECs has grown remarkably. This includes boththe interaction of the constituent particles by employing Feshbach resonances [11, 12, 13] as well as the trapgeometry and dimensionality [14]. This qualifies BECs for testing the fundamental principles of quantummechanics, or very recently even to simulate certain processes in astrophysics [15, 16, 17].Of particular interest is the understanding of the dynamics as well as the spectrum of excited states oftrapped BECs made of interacting bosons. Theoretical descriptions for these were carried out many yearsbefore they were finally realized experimentally. The most prominent approach to obtain the ground state aswell as the dynamical motion of condensates is the Gross-Pitaevskii (GP) equation [18, 19, 20, 21]. In thistheory it is assumed that actually all bosons in the system occupy the same single-particle state, making it amean-field theory. The GP equation has been successfully applied in many cases [22] and can be seen as theadequate starting point to solve many physical problems of BECs. Similarly, the linear-response (LR) theoryatop the GP equation, resulting in the famous Bogoliubov-de Gennes (BdG) equations [23, 24], represents thecommonly utilized approach to describe excitations of trapped ultracold bosons, e.g., for dark-bright solitonpairs [25], in rotating BECs with dipole-dipole interaction [26], or in a 2D hexagonal lattice with superimposedharmonic confinement [27], to name only a few recent applications. Furthermore, the BdG equations wereextensively used for the description of Cooper pairs in superconductors [28] or for the understanding of theBCS-BEC crossover connecting superconductivity and superfluidity in correlated fermionic systems [29].A crucial limitation of the GP mean-field theory is the fact that it ignores correlations between the bosons,and thus does not account for condensate depletion, i.e., the fraction of particles outside the condensed modedue to correlations and interactions, even at zero temperature. In a highly-cited review article from 1999,the condensate depletion is found to be “very small (less than 1%) in the presently available experimentalconditions” [30], and thus the usage of mean-field methods like the BdG approach has been consideredjustified. However, the many-body methods developed in the following years were capable to demonstratethe remarkable impact of depletion and even fragmentation of a BEC onto its dynamics and its excitationspectrum. One of these many-body approaches is the multiconfigurational time-dependent Hartree methodfor bosons [MCTDHB, or MCTDHB( M ) with M being the number of utilized single-particle states, calledorbitals], which has been introduced [31, 32], benchmarked [33], and applied successfully to various systemsin recent years. Being derived from a variational principle, its results can be improved systematically byincreasing the number of orbitals M . Additionally, its LR theory, termed LR-MCTDHB( M ), has beenintroduced recently [34, 35]. The latter represents a full many-body description of excitations in trappedBECs. 5his review has two main objectives. At first, the newly developed numerical implementation of LR-MCTDHB, capable of treating large systems also in D > M = 1, MCTDHB reduces exactly tothe GP equation. Correspondingly, LR-MCTDHB(1) yields the BdG equations. Hence, the commonly usedmean-field approaches are included in the utilized many-body theories as their simplest limiting cases.To address these goals, the review is structured as follows. First, in Section 2, the main theoretical conceptsof many-boson physics are presented. This includes the many-body Schr¨odinger equation, the frameworkof second quantization which is particularly useful for the many-body theory described in this work, theintroduction of the Dirac-Frenkel variational principle from which the central equations of motion can bederived, as well as the definition of Bose-Einstein condensation. Second, in Section 3, a detailed derivation ofMCTDHB is given. The essential ingredients like the ansatz for the many-boson wave function are discussed.Then, its LR theory is derived, which represents the central theory utilized in this work to calculate low-energy excitations of trapped bosonic systems on the many-body level. Third, in Section 4, the newlydeveloped numerical implementation of LR-MCTDHB is described in detail. Therefore, general technicalchallenges and their solutions are explained. Furthermore, a detailed description of the numerical algorithmsutilized is presented. Subsequently, in Section 5, the used implementations of MCTDHB and LR-MCTDHBare benchmarked against the exactly-solvable harmonic-interaction model, whose analytic solutions of theenergies of the ground and excited states are presented beforehand. Applications of both implementations arethen shown and their physics discussed in Section 6. Concerning the dynamics, the focus is twofold. First,applications dealing with the phenomenon of tunneling of BECs between wells inside potential traps in one andtwo spatial dimensions are discussed, demonstrating the many-body nature of quantum tunneling of ultracoldBECs. Second, several cases of dynamical fragmentation of initially coherent bosonic systems are presented.For instance, the tunneling of a BEC to open space is analyzed. Furthermore, the connection between angularmomentum, especially of quantized vortices, and the development of fragmentation is investigated. Concerningapplications of LR-MCTDHB, the many-body effects in the low-energy spectra of BECs in one-dimensionalmulti-well traps and lattices are shown. Moreover, as a 2D application, the impact of angular momentum onthe lowest-in-energy excitations of a rotating BEC is examined. Essential physical results and implicationsare provided. Finally, in Section 7, a summary and an outlook of possible constituent future research as wellas of further enrichment of the implementation of LR-MCTDHB is provided.Additionally, the review contains four appendices which broaden its scope and provide additional insight.In Appendix A, the commonly-used mean-field theories, the GP theory and the constituent LR theory leadingto the BdG equations, are described. A further mean-field approach for both the statics and dynamics oftrapped BECs that has been developed prior to the many-body theories of Section 3 is presented, called thetime-dependent multi-orbital mean field (TDMF) and the linear-response best-mean-field (LR-BMF). Thelatter can be seen as intermediate theories between the standard single-orbital mean-field approaches and thefull many-body treatment and have their own merits. Appendix B demonstrates that LR-MCTDHB can bewritten in block-diagonal form, which, in some cases, can be beneficial. Advantages and possible weaknessesare discussed. In Appendix C, additional results of the LR-MCTDHB benchmark are provided. Finally, inAppendix D, the many-body variance, a quantity that turns out to be very sensitive to many-body effects,is introduced and derived for the example of the center-of-mass (c.m.) position of a BEC. Its relevance isdiscussed in several examples. 6 General concepts of many-boson physics
In this section, the general concepts of describing trapped ultracold bosons are introduced. This includes therepresentation of the many-body Hamiltonian, the bosonic field as well as properties of the wave function fora system of identical bosons. Moreover, the definition of Bose-Einstein condensation and depletion which isused in the subsequent sections is presented.
The time-dependent Schr¨odinger equation for a generic system of N identical spinless bosons in D ≤ H Ψ( r , r , ..., r N ; t ) = i ∂∂t Ψ( r , r , ..., r N ; t ) (1)with the full many-body wave function Ψ( r , r , ..., r N ; t ) depending on the spatial coordinates { r i } of all N bosons and the time t .The Hamiltonian ˆ H reads ˆ H = N (cid:88) i =1 ˆ h ( r i ) + λ N (cid:88) j>i ˆ W ( r i , r j ) (2)with the single-particle term ˆ h ( r i ) = ˆ T ( r i ) + V ( r i ) (3)ˆ T ( r i ) = −
12 ∆ i = − ∂ ∂ r i (4)comprising the contributions of the kinetic energy ˆ T and the external potential V which is assumed to bereal. For simplicity, the terms are given in dimensionless units with (cid:126) = m = 1. The two-body operatorˆ W ( r i , r j ) represents the interaction potential of strength λ . In this work, interaction potentials which solelydepend on the distance r ij = | r i − r j | between two bosons are considered. In general, both V and ˆ W can betime-dependent. In this section, the formalism of second quantization is briefly explained, including the Fock state represen-tation of an N -boson state, the symmetry properties of the wave function, and the Hamiltonian in secondquantization. As a consequence of the spin-statistics theorem [36, 37], the wave function Ψ( r , ..., r N ) of N identical bosonsis a fully symmetrized N -particle state, i.e., it is symmetric under the permutation of two bosons withcoordinates r i and r j , i.e., Ψ( ..., r i , ..., r j , ... ) = +Ψ( ..., r j , ..., r i , ... ) . (5)In second quantization, this is achieved by using a Fock state (or number state) basis of the N -particleHilbert space H ( N ) for the wave function Ψ. This basis incorporates the previously mentioned symmetrizationrequirements as well as the commutation relations of identical bosons shown below.For a complete set of orthonormalized single-particle states, also referred to as orbitals, { φ i : i = 1 , ..., ∞} ,a general Fock state reads | n , n , ... (cid:105) = (cid:89) i (cid:20) √ n i ! (cid:16) ˆ b † i (cid:17) n i (cid:21) | vac (cid:105) (6)7here the integers { n i } denote the occupation numbers of the individual orbitals and | vac (cid:105) the particle vacuum.In Eq. (6), the creation and annihilation operators ˆ b † i and ˆ b i , which create or annihilate a particle in the state φ i , were introduced. The action of these operators on a Fock state is given byˆ b i | n , n , ..., n i , ... (cid:105) = √ n i | n , n , ..., n i − , ... (cid:105) ˆ b † i | n , n , ..., n i , ... (cid:105) = √ n i + 1 | n , n , ..., n i + 1 , ... (cid:105) and thus the number operator for any single-particle state φ i readsˆ n i = ˆ b † i ˆ b i . (7)For a system with a fixed amount of bosons N , the occupation numbers of the orbitals sum up to the totalnumber of particles N , i.e., N = (cid:88) i n i . (8)The creation and annihilation operators fulfill the conventional bosonic commutation relations (cid:104) ˆ b i , ˆ b k (cid:105) = 0 , (cid:104) ˆ b † i , ˆ b † k (cid:105) = 0 , (cid:104) ˆ b i , ˆ b † k (cid:105) = δ ik (9)which ensure the exchange symmetry of the Fock state basis in Eq. (6) and thus of the total wave functionΨ as required in Eq. (5). Utilizing the creation and annihilation operators of bosons for single-particle states from the previous section,one can define operators that create or annihilate a boson at position r at time t . These operators are theso-called field operators, ˆΨ † ( r ) = (cid:88) i ˆ b † i ( t ) φ ∗ i ( r , t ) (10)ˆΨ( r ) = (cid:88) i ˆ b i ( t ) φ i ( r , t ) , (11)where it is taken into account that the orbitals as well as the creation and annihilation operators can in generalbe time-dependent. Thus, ˆ b † ( t ) and ˆ b ( t ) can be expressed in terms of the field operators,ˆ b † i ( t ) = (cid:90) φ i ( r , t ) ˆΨ † ( r ) d r , ˆ b i ( t ) = (cid:90) φ ∗ i ( r , t ) ˆΨ( r ) d r (12)where it is assumed that the set { φ i ( r , t ) } consists of orthonormalized functions. In the following, the timeargument in the above quantities is suppressed. The field operators in Eqs. (10) and (11) fulfill the bosoniccommutation relations (cid:104) ˆΨ( r ) , ˆΨ( r (cid:48) ) (cid:105) = 0 , (cid:104) ˆΨ † ( r ) , ˆΨ † ( r (cid:48) ) (cid:105) = 0 , (cid:104) ˆΨ( r ) , ˆΨ † ( r (cid:48) ) (cid:105) = δ ( r − r (cid:48) ) . (13)The density operator is defined as ˆ ρ ( r ) = ˆΨ † ( r ) ˆΨ( r ) (14)such that the operator of the total number of bosons in the system is given byˆ N = (cid:90) d r ˆ ρ ( r ) . (15)8 .2.3 One- and two-body operators In second quantization, a general one-body operator ˆ A (1) , i.e., an operator acting on a single particle, reducesto the form ˆ A (1) = (cid:88) i,j a ij ˆ b † i ˆ b j (16)where the matrix elements a ij are given by a ij = (cid:90) φ ∗ i ( r ) ˆ A (1) φ j ( r ) d r (17)or, equivalently in bra-ket notation, a ij = (cid:104) i | ˆ A (1) | j (cid:105) . Similarly, the general form of a two-body operator, i.e.,an operator acting on two particles, is given byˆ A (2) = 12 (cid:88) i,j,k,l a ijkl ˆ b † i ˆ b † j ˆ b k ˆ b l (18)with a ijkl = (cid:90) (cid:90) φ ∗ i ( r ) φ ∗ j ( r (cid:48) ) ˆ A (2) φ k ( r ) φ l ( r (cid:48) ) d r d r (cid:48) , (19)or a ijkl = (cid:104) i, j | ˆ A (2) | k, l (cid:105) in short. The factor 1 / The Hamiltonian in Eq. (2) can now be expressed within the framework of second quantization. Usingthe bosonic field operators in Eqs. (10) and (11) as well as the general expressions for one- and two-bodyoperators, Eqs. (16) and (18), one obtainsˆ H = N (cid:88) i =1 ˆ h ( r i ) + λ N (cid:88) j>i ˆ W ( r i , r j ) ⇒ ˆ H = (cid:90) d r (cid:18) ˆΨ † ( r )ˆ h ( r ) ˆΨ( r ) + λ (cid:90) d r (cid:48) ˆΨ † ( r ) ˆΨ † ( r (cid:48) ) ˆ W ( r , r (cid:48) ) ˆΨ( r ) ˆΨ( r (cid:48) ) (cid:19) = (cid:88) i,j h ij ˆ b † i ˆ b j + λ (cid:88) i,j,k,l W ijkl ˆ b † i ˆ b † j ˆ b k ˆ b l (20)with the matrix elements h ij = (cid:90) φ ∗ i ( r ) ˆ h ( r ) φ j ( r ) d r (21) W ijkl = (cid:90) (cid:90) φ ∗ i ( r ) φ ∗ j ( r (cid:48) ) ˆ W ( r , r (cid:48) ) φ k ( r ) φ l ( r (cid:48) ) d r d r (cid:48) . (22)To keep the notation simple, the time argument of the quantities in Eqs. (20)-(22) is suppressed. The aboverepresentation of ˆ H from Eq. (20) will be used in the many-body formalism presented in Section 3. The working equations of all utilized many-body methods in this work can be derived from the Dirac-Frenkelvariational principle (DFVP) [38, 39, 40]. Like for any other variational method, it yields a condition for the9ynamical evolution for any parametrized ansatz of the wave function Ψ. Within the DFVP, the equationthat determines the dynamics is given by (cid:28) δ Ψ( t ) (cid:12)(cid:12)(cid:12)(cid:12) ˆ H − i (cid:126) ∂∂ t (cid:12)(cid:12)(cid:12)(cid:12) Ψ( t ) (cid:29) = 0 (23)which basically means that any allowed variation | δ Ψ( t ) (cid:105) is orthogonal to the residual wave function (cid:16) ˆ H − i (cid:126) ∂∂ t (cid:17) | Ψ( t ) (cid:105) .It is important to note that under certain conditions, the DFVP is fully equivalent to two other commonlyused time-dependent variational principles, which are the variational principle due to MacLachlan [41] andthe least action principle. See in this context Ref. [42]. The least action principle is based on the assumptionthat the wave function minimizes the action S = (cid:90) dt (cid:28) Ψ (cid:12)(cid:12)(cid:12)(cid:12) ˆ H − i (cid:126) ∂∂t (cid:12)(cid:12)(cid:12)(cid:12) Ψ (cid:29) , (24)meaning that its variation with respect to Ψ should vanish, i.e., δS = 0 . (25)Eqs. (24) and (25), in combination with normalization constraints, are used to derive the equations of motionfor MCTDHB. The full derivation of the latter theory is carried out in detail in Section 3.1. Furthermore, theleast action principle can also be employed to derive the standard mean-field equation for the ground stateand dynamics of a trapped BEC, i.e., the GP equation. Its derivation is presented in Appendix A.1.1. In 1956, Penrose and Onsager developed a theoretical criterion which is utilized in the definition of Bose-Einstein condensation and is based on standard quantities from statistical quantum mechanics [43]. It makesuse of the one-body reduced density matrix (RDM) of an N -boson system which can be defined via the firstorder correlation function ρ (1) ( r | r (cid:48) ) = (cid:68) Ψ (cid:12)(cid:12)(cid:12) ˆΨ † ( r (cid:48) ) ˆΨ( r ) (cid:12)(cid:12)(cid:12) Ψ (cid:69) = (cid:88) i,j ρ ij φ ∗ i ( r (cid:48) ) φ j ( r ) (26)with ρ ij = (cid:104) Ψ | ˆ ρ ij | Ψ (cid:105) , ˆ ρ ij = ˆ b † i ˆ b j (27)where the time argument is suppressed in all quantities. More details on RDMs are given in Appendix D. For r = r (cid:48) , Eq. (26) denotes the one-particle density. Diagonalizing the one-body RDM yields ρ (1) ( r | r (cid:48) ) = (cid:88) i n (1) i α ∗ i ( r (cid:48) ) α i ( r ) (28)where the eigenvalues (cid:110) n (1) i | n (1)1 ≥ n (1)2 ≥ ... (cid:111) and the eigenfunctions { α i ( r ) } are the natural occupation num-bers and the natural orbitals, respectively. Whenever possible, the superscript ’(1)’ of the natural occupationnumbers will be omitted throughout this work. According to Ref. [43], any N -boson system where onlythe largest eigenvalue is macroscopic, i.e., n = O ( N ), is said to be condensed. However, if more than oneeigenvalue of the one-body RDM are macroscopic, the system is said to be depleted/fragmented [44, 45]. Inthis work, the degree of depletion (or fragmentation) is quantified by f = 1 N (cid:88) i> n i (29)10hich is the aggregated occupation of all but the first natural orbital. Another classification that uses theabove relation distinguishes the terms depletion and fragmentation more clearly and was used in Ref. [46].In Section 6, the dynamics and excited states of BECs in different trap geometries in one and two spatialdimensions are studied. There, the degree of depletion/fragmentation will be an essential quantity to dis-tinguish the many-body results obtained from the methods described in the subsequent Section 3 from themean-field results obtained by using the theories described in the Appendices A.1.1 and A.1.2. In this section, the many-body approach utilized in this work to compute the ground state and the dynamics oftrapped BECs as well as their low-energy excitation spectra is presented. First, an introduction to MCTDHBis made in Section 3.1, followed by the application of LR theory atop it, termed LR-MCTDHB, in Section3.2.2.
In this section, a detailed derivation of the MCTDHB [or MCTDHB( M )] theory which has been originallyintroduced in Refs. [31, 32], is presented. It originates from the closely related MCTDH method [47, 48, 49]which is a widely used technique to compute the dynamics and excited states of nuclei in molecules. Furtherextensions of MCTDHB to the case of type conversion of particles, to mixtures of bosons and fermions, toHubbard-type Hamiltonians of BECs on lattices, as well as to BECs with internal degrees of freedom canbe found in Refs. [50, 51, 52, 53]. Other many-body theories for the ground state and excitations of BECswere utilized as well. One example is the density matrix renormalization group (DMRG) theory [54]. Itis particularly designed for applications in the field of one-dimensional condensed matter physics, but wasextended beyond this to quantum chemistry and quantum information. An early review can be found in [55].Furthermore, DMRG has been extended and applied to two-dimensional systems as well [56]. However, thisreview deals with MCTDHB and its LR theory.The ansatz for the many-body wave function to solve the many-boson Schr¨odinger equation, Eq. (1), isgiven by | Ψ( t ) (cid:105) = (cid:88) n C n ( t ) | n ; t (cid:105) , (30)which is a superposition of permanents {| n ; t (cid:105)} comprised of M single-particle orbitals { φ j ( r , t ) : 1 ≤ j ≤ M } and expansion coefficients { C n ( t ) } where n = ( n , . . . , n M ) t is a vector carrying the individual occupationnumbers of the orbitals for a given permanent. The summation in Eq. (30) includes all N conf = (cid:0) N + M − N (cid:1) possibilities to distribute N bosons onto M single-particle orbitals. It is important to note that both thepermanents and the coefficients are time-dependent, which is the main difference compared to the frequentlyused many-body approach of the full configuration interaction (FCI) method. In the latter theory, the amountof configurations for a given number M of single-particle basis functions is the same as for MCTDHB. However,these basis states are fixed in shape, i.e., they remain unchanged and only the coefficients are time-dependent.For MCTDHB, the single-particle orbitals are time-adaptive and, as shown below, the corresponding setof equations of motion (EOMs) is coupled to the set of EOMs of the coefficients. In terms of numericalconvergence, utilizing a time-adaptive instead of a shape-fixed basis dramatically improves the convergencetowards exact results for both the ground state and the dynamics of trapped BECs, as discussed in Section 5.The working equations of MCTDHB are determined by the least action principle of Eqs. (24) and (25),which explicitly reads 0 = δS = δ (cid:18)(cid:90) dt L ( t ) (cid:19) (31)11ith L ( t ) = (cid:28) Ψ( t ) (cid:12)(cid:12)(cid:12)(cid:12) ˆ H − i ∂∂t (cid:12)(cid:12)(cid:12)(cid:12) Ψ( t ) (cid:29) − (cid:88) i,j µ ij ( t )[ (cid:104) φ i | φ j (cid:105) − δ ij ] (32)and (cid:126) = 1. The time-dependent Langrange multipliers { µ ij ( t ) } account for the orthonormalization of theorbitals in time. Whenever possible, the time-argument of all quantities is suppressed in the following.Plugging in the ansatz for the wave function, Eq. (30), and the expression for the Hamiltonian, Eq. (20), intoEq. (32) results in L ( t ) = M (cid:88) k,q ρ kq (cid:34) h kq − (cid:18) i ∂∂t (cid:19) kq (cid:35) + 12 M (cid:88) k,s,q,l ρ ksql W ksql − i (cid:88) n C ∗ n ( t ) ∂C n ( t ) ∂t − (cid:88) k,q µ kq [ (cid:104) φ k | φ q (cid:105) − δ kq ] (33)where the elements of the one- and two-body RDMs in this case read ρ kq = (cid:88) n , n (cid:48) C ∗ n ( t ) C n (cid:48) ( t ) (cid:68) n ; t (cid:12)(cid:12)(cid:12) ˆ b † k ˆ b q (cid:12)(cid:12)(cid:12) n (cid:48) ; t (cid:69) (34) ρ ksql = (cid:88) n , n (cid:48) C ∗ n ( t ) C n (cid:48) ( t ) (cid:68) n ; t (cid:12)(cid:12)(cid:12) ˆ b † k ˆ b † s ˆ b q ˆ b l (cid:12)(cid:12)(cid:12) n (cid:48) ; t (cid:69) . (35)The variation of S with respect to φ ∗ k ( r , t ) yields δSδφ ∗ k ( r , t ) = M (cid:88) q =1 ρ kq (cid:18) ˆ h − i ∂∂t (cid:19) − µ kq + M (cid:88) s,l =1 ρ ksql ˆ W sl | φ q (cid:105) (36)with ˆ W sl ( r , t ) = (cid:90) φ ∗ s ( r (cid:48) , t ) ˆ W ( r , r (cid:48) ) φ l ( r (cid:48) , t ) d r (cid:48) . (37)According to the DFVP, this variation should vanish, which determines the EOMs of the orbitals given by M (cid:88) q =1 ρ kq ˆ h − µ kq + M (cid:88) s,l =1 ρ ksql ˆ W sl | φ q (cid:105) = i M (cid:88) q =1 ρ kq ∂∂t | φ q (cid:105) . (38)Taking the scalar product with (cid:104) φ i | in Eq. (38) results in the expression for the Lagrange multipliers, µ ki ( t ) = M (cid:88) q =1 ρ kq (cid:34) h iq − (cid:18) i ∂∂t (cid:19) iq (cid:35) + M (cid:88) s,l =1 ρ ksql W isql . (39)By reinserting them into Eq. (38), multiplying by the inverse of the one-body RDM { ρ − } ik from the left-handside (LHS), and summing over k , one arrives at i ˆ P ∂∂t | φ i (cid:105) = ˆ P ˆ h | φ i (cid:105) + M (cid:88) k,s,l,q =1 { ρ − } ik ρ ksql ˆ W sl | φ q (cid:105) (40)with ˆ P = − M (cid:88) j (cid:48) =1 | φ j (cid:48) (cid:105)(cid:104) φ j (cid:48) | (41)12eing a projection operator on the tangential space of span ( φ , ..., φ M ). Due to invariance properties of theansatz for the wave function in Eq. (30), the orbitals can be chosen without loss of generality to be orthogonalto their corresponding time-derivatives [47, 48], i.e., (cid:104) φ i | ˙ φ j (cid:105) = 0 ∀ i, j = 1 , ..., M. (42)The time evolution of the orbitals is therefore fully orthogonal and ensures the orthonormality of the orbitalsin time. Hence, one can further simplify Eq. (40) to read i ∂∂t | φ i (cid:105) = ˆ P ˆ h | φ i (cid:105) + M (cid:88) k,s,l,q =1 { ρ − } ik ρ ksql ˆ W sl | φ q (cid:105) , (43)and Eq. (43) is referred to as the orbitals’ EOM in the following.With respect to the coefficients, the least action principle δSδC ∗ n ( t ) = 0 (44)leads to (cid:88) n (cid:48) (cid:28) n ; t (cid:12)(cid:12)(cid:12)(cid:12) ˆ H − i ∂∂t (cid:12)(cid:12)(cid:12)(cid:12) n (cid:48) ; t (cid:29) C n (cid:48) ( t ) = i ∂∂t C n ( t ) . (45)Thus, the EOM of the coefficients is given by H ( t ) C ( t ) = i ∂∂t C ( t ) (46)with the matrix elements H nn (cid:48) = (cid:28) n ; t (cid:12)(cid:12)(cid:12)(cid:12) ˆ H − i ∂∂t (cid:12)(cid:12)(cid:12)(cid:12) n (cid:48) ; t (cid:29) (47)and C ( t ) being a vector whose entries are the N conf coefficients { C n ( t ) } . It is worth noting that H ( t ) ishermitian. Eq. (46) is referred to as the coefficients’ EOM in the following.One observes that the sets of working equations for the orbitals and coefficients are coupled to each other.With regard to the orbitals’ EOM, the coefficients appear in the expressions for the two-body RDM as well asin the inverse of the one-body RDM, see Eqs. (34), (35) and (43), respectively. With regard to the coefficients’EOM, the implicit dependence on the orbitals is incorporated in the matrix elements of Eq. (47) since thepermanents appearing therein are assembled by the orbitals.By employing imaginary time-propagation, i.e., by setting it → τ in the working equations (43) and (46),the MCTDHB theory reduces to the multiconfigurational Hartree for bosons (MCHB) method which has beenintroduced in Ref. [57]. In that way, the ground state of a given system can be computed very efficiently, aswill be shown in Section 5.2.It can be inferred from Eq. (46) that the evolution of C ( t ) is unitary and therefore its normalization isconserved in time. However, in comparison to the orbitals’ EOMs in Eq. (43), the coefficient vector C ( t )does not automatically propagate in an orthogonal manner. The latter can be achieved by introducing atime-dependent phase of the form C ( t ) → C ( t ) e − i (cid:82) t dt (cid:48) C † ( t (cid:48) ) H ( t (cid:48) ) C ( t (cid:48) ) (48)which, by inserting it into Eq. (46), yields i ∂∂t C ( t ) = P c H ( t ) C ( t ) (49)13ith the projector that maps onto the orthogonal space of the coefficients is given by P c = − C ( t ) C † ( t ) . (50)Including the projector P c is however redundant in the sense that the evolution of the coefficient vector isunitary already without it, see Eq. (46). On the contrary, this is not the case for the orbitals, where theprojector ˆ P that stems from the Lagrange multipliers in Eq. (31) ensures that the orbitals stay normalizedin time. Henceforth, only the latter projector is considered in the subsequent derivation of LR-MCTDHB(Section 3.2.2), whereas the coefficients’ projector P c is set to unity in the following without any loss ofgenerality.It is important to note that for a single orbital, i.e., M = 1, the MCTDHB theory reduces to the GPmean-field equation [18, 19, 20, 21], meaning that GP ≡ MCTDHB(1). The derivation of the latter theoryusing the variational principle is described in Appendix A.1.1.
Before carrying out the derivation of the many-body LR theory atop MCTDHB, a rather general perspectiveof applying LR to the time-dependent Schr¨odinger equation, Eq. (1), in order to calculate the excitationspectrum of a trapped BEC is discussed. The goal is to demonstrate that in general, a LR analysis atop anexact eigenstate of the unperturbed Hamiltonian ˆ H , typically the ground state, yields the exact excitationspectrum of the many-particle system. The following derivation is closely related to the one in Ref. [35].One considers the time-dependent Hamiltonianˆ H ( t ) = ˆ H + ˆ H ext ( t ) , (51)ˆ H ext ( t ) = ˆ f + ( r ) e − iωt + ˆ f − ( r ) e iωt (52)where a weak time-dependent external field with amplitudes ˆ f + and ˆ f − and oscillation frequency ω is addedto the stationary hermitian Hamiltonian ˆ H . The latter has eigenenergies E k , k ∈ N , and the ground-stateenergy is denoted by ε . The projected time-dependent Schr¨odinger equation with (cid:126) = 1,ˆ P Ψ ˆ H ( t )Ψ( t ) = i ˙Ψ( t ) , ˆ P Ψ = − | Ψ( t ) (cid:105)(cid:104) Ψ( t ) | , (53)which can be obtained by introducing the phaseΨ( t ) → Ψ( t ) e − i (cid:82) t dt (cid:48) (cid:104) Ψ( t (cid:48) ) | ˆ H ( t (cid:48) ) (cid:105) Ψ( t (cid:48) ) (54)to the system’s wave function Ψ( t ) and inserting it into Eq. (1), ensures that the time evolution is orthogonalin the sense that (cid:104) Ψ( t ) | ˙Ψ( t ) (cid:105) = 0 (55)at any time t . To solve Eq. (53), the ansatzΨ( r , t ) = e − ε t (cid:0) Ψ ( r ) + u ( r ) e − iωt + v ∗ ( r ) e iωt (cid:1) (56)with is employed. The ground state Ψ obtained by ˆ H Ψ = ε Ψ as well as the correction amplitudes u and v are assumed to be stationary. The attached time-dependent oscillations have the same frequency ω as theexternal perturbing field. Utilizing the orthogonality condition, Eq. (55), one obtains (cid:104) Ψ | u (cid:105) = 0 and (cid:104) Ψ | v ∗ (cid:105) = 0 , (57)14eaning that the correction amplitudes are orthogonal to the ground-state wave function and thusˆ P Ψ | u (cid:105) = | u (cid:105) and ˆ P Ψ | v ∗ (cid:105) = | v ∗ (cid:105) (58)for the operator projecting onto the tangential space of Ψ , i.e., ˆ P Ψ = − | Ψ (cid:105)(cid:104) Ψ | . Inserting the ansatz(56) into Eq. (53), one arrives atˆ P Ψ ˆ H ( u e − iωt + v ∗ e iωt ) + ˆ P Ψ ( ˆ f + e − iωt + ˆ f − e iωt )Ψ = ( ω + ε )( u e − iωt − v ∗ e iωt ) (59)where only terms linear to u and v are kept. This can be written in matrix form by ordering terms proportionalto e − iωt and e iωt , respectively. The result reads (cid:20)(cid:18) ˆ P Ψ ( ˆ H − ε ) 00 − ˆ P ∗ Ψ ( ˆ H , ∗ − ε ) (cid:19) − ω (cid:21) (cid:18) uv (cid:19) = (cid:18) − ˆ P Ψ ˆ f + Ψ ˆ P ∗ Ψ ˆ f − , ∗ Ψ (cid:19) (60)where complex conjugation to the lower equation has been applied. Defining the LR matrix L as L = (cid:18) ˆ P Ψ ( ˆ H − ε ) ˆ P Ψ − ˆ P ∗ Ψ ( ˆ H , ∗ − ε ) ˆ P ∗ Ψ (cid:19) (61)where redundantly the projectors ˆ P Ψ and ˆ P ∗ Ψ have been inserted from the right-hand side (RHS) to theupper and lower blocks, the homogeneous version of Eq. (60) is given by L (cid:18) uv (cid:19) = ω (cid:18) uv (cid:19) . (62)For the upper block, one can immediately deduce that the spectrum of corresponding eigenvalues is given by ω k = E k − ε , meaning that one obtains the eigenvalues of the unperturbed Hamiltonian ˆ H . Moreover, oneobtains ω k = − ( E k − ε ) for the lower block.This is an important result because it means that the LR due to an external perturbation, as described inEq. (52), yields, if applied to an exact eigenstate of ˆ H (not necessarily the ground state), the full spectrumof exact excitation energies. It is worth stressing the generality of the above derivation because no physicalproperties of the system were specified, i.e., whether ˆ H describes a single- or multi-particle system, bosons orfermions (or even mixtures), or identical or distinguishable particles. Thus, LR represents a generic methodfor obtaining the spectrum of excited states of a quantum system. Section 3.2.2, explicitly deals with themany-body linear response of systems consisting of identical bosons.This section is closed by referring to systems where no exact eigenstate of ˆ H , e.g., the ground state, isknown. This is naturally the most common case since the number of systems where analytic solutions areavailable is very small. Hence, in order to obtain accurate results for the excitation energies from a LR analysis,one is in need of a very accurate approximation of the system’s ground state. In other words, increasing thequality of the underlying ground-state approximation will increase the accuracy of the spectrum of excitedstates. From a numerical point of view, Section 5 shows that MCTDHB turns out to be a very powerfulapproach to calculate the ground state of a trapped BEC to very high accuracy, and that it is clearly superiorto other widely used methods like the GP mean-field equation or the FCI method. As a result, LR-MCTDHBleads to highly accurate excitation spectra, for which numerical evidence is presented in Section 5.3. In this section, LR-MCTDHB [or LR-MCTDHB( M )] is derived. To this end, a trapped N -boson system isconsidered and its ground-state orbitals and coefficients, { φ k : 1 ≤ k ≤ M } and C , are calculated by usingimaginary time-propagation of the MCTDHB working equations [57], Eqs. (43) and (46). Afterwards, the15ame time-dependent external perturbation as given in Eq. (52) is added to the Hamiltonian, such that thelatter is finally given by Eqs. (51). The EOMs are linearized according to this. The subsequent derivation ofthe LR equations follows Refs. [34, 35].The ans¨atze for the response orbitals and coefficients due to the external perturbation are φ k ( r , t ) ≈ φ k ( r ) + δφ k ( r , t ) (63) δφ k ( r , t ) = u k ( r ) e − iωt + v ∗ k ( r ) e iωt (64) C ( t ) ≈ e − i(cid:15) t [ C + δ C ( t )] (65) δ C ( t ) = C u e − iωt + C ∗ v e iωt (66)with 1 ≤ k ≤ M . The k -th response orbital is thus given by the zeroth-order orbital, typically the time-independent ground-state orbital φ k ( r ), plus a weak, time-dependent perturbation δφ k ( r , t ) with oscillationfrequency ω and stationary response amplitudes u k ( r ) and v k ( r ). The response coefficients consist of sim-ilar parts and are additionally multiplied with a time-dependent phase e − iε t to which it is referred below.The time-dependence is therefore essentially incorporated in the perturbation part. In the following, thelinearization of the orbitals’ and coefficients’ EOMs are discussed separately.For linearizing the orbitals’ EOMs, it is instructive to utilize the version in Eq. (38) given by i M (cid:88) q =1 ρ kq | ˙ φ q ( t ) (cid:105) = M (cid:88) q =1 [ ˆ Z kq − µ kq ( t )] | φ q ( t ) (cid:105) , (67)with the replacement ˆ Z kq = ρ kq [ˆ h + ˆ H ext ( t )] + M (cid:88) s,l =1 ρ ksql ˆ W sl . (68)Inserting Eq. (63) into Eq. (67), one can group the resulting expression in terms of perturbation orders. Inzeroth order, one obtains 0 = (cid:88) q [ ˆ Z kq − µ kq ] | φ q (cid:105) (69)with ˆ Z kq = ρ kq ˆ h + M (cid:88) s,l =1 ρ ksql ˆ W sl . (70)Here and in the following, quantities with the superscript ’0’ denote that they only contain zeroth-orderorbitals and coefficients. Eq. (69) represents the orbital part of the MCHB working equations, see Ref. [57]in this context. For the first order, the LHS of Eq. (67) yields i (cid:88) q ( δρ kq | ˙ φ q (cid:105) (cid:124)(cid:123)(cid:122)(cid:125) =0 + ρ kq | ˙ δφ q ( t ) (cid:105) ) = i (cid:88) q ρ kq | ˙ δφ q ( t ) (cid:105) . (71)The RHS is given by (cid:88) q [ ˆ Z kq − µ kq ] | δφ q ( t ) (cid:105) + (cid:88) q [ δ ˆ Z kq − δµ kq ( t )] | φ q (cid:105) (72)with δ ˆ Z kq = δ ρ kq ˆ h + (cid:88) s,l ρ ksql ˆ W sl + ρ kq ˆ H ext ( t )= δρ kq ˆ h + (cid:88) s,l ρ ksql δ ˆ W sl + (cid:88) s,l δρ ksql ˆ W sl + ρ kq ˆ H ext ( t ) . (73)16o compute the RHS further, the variation of the chemical potential δµ kq ( t ), given by δµ kq ( t ) = (cid:42) δφ q (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:88) j (cid:20) ˆ Z kj − iρ kj ∂∂t (cid:21) (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) φ j (cid:43) + (cid:42) φ q (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:88) j δ (cid:20) ˆ Z kj − iρ kj ∂∂t (cid:12)(cid:12)(cid:12)(cid:12) φ j (cid:43) Eq.(69) = (cid:88) j µ kj (cid:10) δφ q | φ j (cid:11) + (cid:42) φ q (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:88) j ρ kj ˆ H ext ( t ) (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) φ j (cid:43) + (cid:42) φ q (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:88) j δρ kj ˆ h + (cid:88) j,s,l ( ρ ksjl δ ˆ W sl + δρ ksjl ˆ W sl ) (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) φ j (cid:43) + (cid:42) φ q (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:88) j (cid:18) ˆ Z kj − iρ kj ∂∂t (cid:19)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) δφ j (cid:43) , (74)is employed. Thus, Eq. (72) can be written as (cid:88) q [ ˆ Z kq − µ kq ] | δφ q (cid:105) + (cid:88) q δρ kq ˆ h + ρ kq ˆ H ext ( t ) + (cid:88) s,l ( ρ ksql δ ˆ W sl + δρ ksql ˆ W sl ) | φ q (cid:105)− (cid:88) q (cid:88) j ( − µ kj ) (cid:10) φ q | δφ j (cid:11) + (cid:42) φ q (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:88) j ρ kj ˆ H ext ( t ) (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) φ j (cid:43) + (cid:42) φ q (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:88) j δρ kj ˆ h + (cid:88) j,s,l ( ρ ksjl δ ˆ W sl + δρ ksjl ˆ W sl ) (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) φ j (cid:43) + (cid:42) φ q (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:88) j (cid:18) ˆ Z kj − iρ kj ∂∂t (cid:19)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) δφ j (cid:43) | φ q (cid:105) (75)which, combined with the LHS, Eq. (71), leads to (cid:88) q ˆ P (cid:110)(cid:16) ˆ Z kq − µ kq (cid:17) | δφ q (cid:105) + (cid:16) δρ kq ˆ h + ρ kq ˆ H ext ( t )+ (cid:88) s,l δρ ksql ˆ W sl + (cid:88) s,l ρ ksql δ ˆ W sl (cid:17) | φ q (cid:105) (cid:111) = i ˆ P (cid:88) q ρ kq | δ ˙ φ q (cid:105) (76)where the projection operator ˆ P = − (cid:80) j (cid:48) | φ j (cid:48) (cid:105)(cid:104) φ j (cid:48) | projects on the tangential space of the stationary orbitals.From the invariance condition, Eq. (42), one can deduce0 = δ (cid:104) φ k | ˙ φ q (cid:105) = (cid:104) δφ k | ˙ φ q (cid:105) (cid:124)(cid:123)(cid:122)(cid:125) =0 + (cid:104) φ k | δ ˙ φ q (cid:105) = − iω (cid:104) φ k | u q (cid:105) e − iωt + iω (cid:104) φ k | v ∗ q (cid:105) e iωt ⇒ (cid:104) φ k | u q (cid:105) = (cid:104) φ k | v ∗ q (cid:105) = 0 ∀ ≤ k, q ≤ M (77)such that ˆ P | δ ˙ φ q (cid:105) = − M (cid:88) j (cid:48) =1 | φ j (cid:48) (cid:105)(cid:104) φ j (cid:48) | | δ ˙ φ q (cid:105) = | δ ˙ φ q (cid:105) (78)17hich makes the projection operator on the RHS of Eq. (76) redundant. One further utilizes δ ˆ W sl [ u l , u ∗ s , v ∗ l , v s ] = (cid:90) d r (cid:48) (cid:110) δφ ∗ s ( r (cid:48) , t ) ˆ W ( r , r (cid:48) ) φ l ( r (cid:48) )+ φ , ∗ s ( r (cid:48) ) ˆ W ( r , r (cid:48) ) δφ l ( r (cid:48) , t ) (cid:111) (79) δρ kq [ C u , C v , C ∗ u , C ∗ v ] = (cid:104) δ C ( t ) | ˆ b † k ˆ b q | C (cid:105) + (cid:104) C | ˆ b † k ˆ b q | δ C ( t ) (cid:105) (80) δρ ksql [ C u , C v , C ∗ u , C ∗ v ] = (cid:104) δ C ( t ) | ˆ b † k ˆ b † s ˆ b q ˆ b l | C (cid:105) + (cid:104) C | ˆ b † k ˆ b † s ˆ b q ˆ b l | δ C ( t ) (cid:105) (81)where the action of, e.g., ˆ b † k ˆ b q on | C (cid:105) results in a new vector of coefficients for the same Fock states, i.e.,ˆ b † k ˆ b q | C (cid:105) = (cid:88) n C n ( t )ˆ b † k ˆ b q | n ; t (cid:105) = (cid:88) n (cid:48) C n (cid:48) ( t ) | n (cid:48) ; t (cid:105) ≡ | C (cid:48) (cid:105) , (82)and inserts the above equations, together with Eqs. (52), (64), and (66), into Eq. (76). Collecting all termsproportional to e − iωt results in (cid:88) q ˆ P (cid:110)(cid:16) ˆ Z kq − µ kq (cid:17) | u q (cid:105) + (cid:16) δρ kq [ C u , C v ] ˆ h + ρ kq ˆ f + + (cid:88) s,l δρ ksql [ C u , C v ] ˆ W sl + (cid:88) s,l ρ ksql δ ˆ W sl [ u l , v s ] (cid:17) | φ q (cid:105) (cid:111) = ω (cid:88) q ρ kq | u q (cid:105) . (83)Similarly, collecting all terms proportional to e iωt yields (cid:88) q ˆ P (cid:110)(cid:16) ˆ Z kq − µ kq (cid:17) | v ∗ q (cid:105) + (cid:16) δρ kq [ C ∗ u , C ∗ v ] ˆ h + ρ kq ˆ f − + (cid:88) s,l δρ ksql [ C ∗ u , C ∗ v ] ˆ W sl + (cid:88) s,l ρ ksql δ ˆ W sl [ u ∗ s , v ∗ l ] (cid:17) | φ q (cid:105) (cid:111) = − ω (cid:88) q ρ kq | v ∗ q (cid:105) C.C. ⇐⇒ (cid:88) q ˆ P ∗ (cid:110)(cid:16) ˆ Z , ∗ kq − µ qk (cid:17) | v q (cid:105) + (cid:16) δρ qk [ C u , C v ] ˆ h ∗ + ρ qk ˆ f − , ∗ + (cid:88) s,l δρ qlks [ C u , C v ] ˆ W ls + (cid:88) s,l ρ qlks δ ˆ W ls [ u s , v l ] (cid:17) | φ , ∗ q (cid:105) (cid:111) = − ω (cid:88) q ρ qk | v q (cid:105) (84)where ’C.C.’ denotes complex conjugation of the entire equation. Here, it is used that the matrix { µ ij } ofthe Lagrange multipliers is hermitian. Furthermore, the two-body interaction potential is assumed to be realthroughout this work.For linearizing the coefficients’ EOMs, one obtains in zeroth order the stationary equation H C = ε C (85)with the matrix elements {H nn (cid:48) } = (cid:104) n | ˆ H | n (cid:48) (cid:105) . Eq. (85) is the working equation for the coefficients withinthe MCHB theory. In first order, one arrives at i ∂∂t δ C ( t ) = δ H C + (cid:0) H − ε (cid:1) δ C ( t ) (86)18ith δ H = (cid:88) k,q δh kq ˆ b † k ˆ b q + 12 (cid:88) k,s,q,l δW ksql ˆ b † k ˆ b † s ˆ b q ˆ b l (87)and δh kq = (cid:90) d r δφ ∗ k ( r , t )ˆ hφ q ( r ) + (cid:90) d r φ , ∗ k ( r )ˆ hδφ q ( r , t )+ (cid:90) d r φ , ∗ k ( r ) ˆ H ext ( r , t ) φ q ( r ) (cid:124) (cid:123)(cid:122) (cid:125) ≡ ( H ext ) kq (88) δW ksql = (cid:90) (cid:90) d r d r (cid:48) δφ ∗ k ( r , t ) φ , ∗ s ( r (cid:48) ) ˆ W ( r , r (cid:48) ) φ q ( r ) φ l ( r (cid:48) )+ (cid:90) (cid:90) d r d r (cid:48) φ , ∗ k ( r ) δφ ∗ s ( r (cid:48) , t ) ˆ W ( r , r (cid:48) ) φ q ( r ) φ l ( r (cid:48) )+ (cid:90) (cid:90) d r d r (cid:48) φ , ∗ k ( r ) φ , ∗ s ( r (cid:48) ) ˆ W ( r , r (cid:48) ) δφ q ( r , t ) φ l ( r (cid:48) )+ (cid:90) (cid:90) d r d r (cid:48) φ , ∗ k ( r ) φ , ∗ s ( r (cid:48) ) ˆ W ( r , r (cid:48) ) φ q ( r ) δφ l ( r (cid:48) , t ) . (89)Collecting all terms proportional to e − iωt yields ω C u = ( H − ε ) C u + (cid:88) k,q δh k,q [ v k , u q , ˆ f + ]ˆ b † k ˆ b q + 12 (cid:88) k,s,q,l δW ksql [ v k , v s , u q , u l ]ˆ b † k ˆ b † s ˆ b q ˆ b l C , (90)whereas collecting all terms proportional to e + iωt results in − ω C ∗ v = ( H − ε ) C ∗ v + (cid:88) k,q δh k,q [ u ∗ k , v ∗ q , ˆ f − ]ˆ b † k ˆ b q + 12 (cid:88) k,s,q,l δW ksql [ u ∗ k , u ∗ s , v ∗ q , v ∗ l ]ˆ b † k ˆ b † s ˆ b q ˆ b l C ⇐⇒ − ω C v = ( H , ∗ − ε ) C v + (cid:88) k,q δh q,k [ u k , v q , ˆ f − , ∗ ] (cid:16) ˆ b † k ˆ b q (cid:17) ∗ + 12 (cid:88) q,l,k,s δW qlks [ u k , u s , v q , v l ] (cid:16) ˆ b † k ˆ b † s ˆ b q ˆ b l (cid:17) ∗ C , ∗ . (91)After the linearization of both the orbitals’ and coefficients’ EOMs, one can combine the sets of LRequations by casting Eqs. (83), (84), (90) and (91) into matrix form, leading to( PL − M ω ) uvC u C v = M P R [ ˆ f + , ˆ f − ] (92)19ith u = ( u , ..., u M ) t , v = ( v , ..., v M ) t and the LR matrix L = (cid:18) L oo L oc L co L cc (cid:19) (93)where the submatrices refer to the couplings between the orbitals and the coefficients. In Eq. (92), P is aprojection matrix, M is a metric which contains the one-body RDM, and the RHS R contains the amplitudesof the external perturbation, ˆ f + and ˆ f − . In the following, the submatrices in Eq. (93) are explicitly derived.The (2 M )-dimensional square orbital matrix L oo consists of four blocks and reads L oo = (cid:18) L uoo L voo − ( L voo ) ∗ − ( L uoo ) ∗ (cid:19) . (94)The upper part can be obtained from Eq. (83) by collecting all terms proportional to u , yielding L uoo = ρ ˆ h − µ + Ω + κ (95)with ρ = { ρ kq } , µ = { µ kq } , Ω = { Ω kq } = (cid:80) s,l ρ ksql ˆ W sl and κ = { κ kq } = (cid:80) s,l ρ ksql ˆ K sl . Therein, theexchange interaction operator is defined asˆ K sl = (cid:90) d r (cid:48) φ ∗ s ( r (cid:48) ) ˆ W ( r , r (cid:48) ) ˆ P rr (cid:48) φ l ( r (cid:48) ) (96)where ˆ P rr (cid:48) permutes the coordinates r and r (cid:48) , i.e.,ˆ P rr (cid:48) f ( r ) = f ( r (cid:48) ) (97)for any function f ( r ). The action of ˆ K sl on an arbitrary function f ( r ) therefore readsˆ K sl f ( r ) = ˆ W sf φ l ( r ) . (98)Furthermore, collecting all terms proportional to v yields L voo = κ , κ = { κ kq } = (cid:88) s,l ρ kqsl ˆ K l ∗ s (99)for the upper right block of L oo . The lower two submatrices can be obtained from Eq. (84) and turn out tobe the negative complex conjugate of the upper submatrices.The upper right block of L , i.e., the 2( M × N conf )-matrix L oc , has a similar structure than L oo , meaningit consists of four separate blocks given by L oc = (cid:18) L uoc L voc − ( L voc ) ∗ − ( L uoc ) ∗ (cid:19) . (100)The upper two submatrices can again be extracted from Eq. (83) and yield, collecting all terms proportionalto C u , L uoc = (cid:88) q ˆ h (cid:0) ˆ ρ qk C (cid:1) † + (cid:88) s,l ˆ W sl (cid:0) ˆ ρ qlks C (cid:1) † φ q , (101)whereas collecting all terms proportional to C v results in L voc = (cid:88) q ˆ h (cid:0) ˆ ρ kq C (cid:1) t + (cid:88) s,l ˆ W sl (cid:0) ˆ ρ ksql C (cid:1) t φ q . (102)20he lower two submatrices can be extracted from Eq. (84), see Ref. [35].With respect to the 2( N conf × M )-matrix L co , one finds L co = (cid:18) L uco L vco − ( L vco ) ∗ − ( L uco ) ∗ (cid:19) (103)with L uco = (cid:88) k φ , ∗ k (cid:40)(cid:0) ˆ ρ kq C (cid:1) ˆ h + (cid:88) sl (cid:0) ˆ ρ klqs C (cid:1) ˆ W ls (cid:41) (104)from collecting all terms proportional to u q in Eq. (90) and L vco = (cid:88) k φ k (cid:40)(cid:0) ˆ ρ qk C (cid:1) ˆ h ∗ + (cid:88) sl (cid:0) ˆ ρ qskl C (cid:1) ˆ W sl (cid:41) (105)from collecting all terms proportional to v q in Eq. (91). It can thus be seen that( L uoc ) † = L uco , ( L voc ) t = L vco . (106)Furthermore, one can extract the (2 N conf )-dimensional square matrix L cc from Eqs. (90) and (91), whichgives L cc = (cid:18) H − ε cc cc − (cid:0) H , ∗ − ε (cid:1)(cid:19) (107)where cc is the ( N conf )-dimensional zero matrix.The projection operator P is given by P = (cid:18) P oo oc co cc (cid:19) (108)with P oo = (cid:18) ˆ P ˆP ∗ (cid:19) (109)and cc being the (2 N conf )-dimensional unit matrix. The matrices oc and co denote the 2( M × N conf )- and2( N conf × M )-dimensional zero matrices, respectively.The metric M reads M = (cid:18) ρ oo oc co cc (cid:19) (110)where ρ oo is given by ρ oo = (cid:18) ρ o o ρ , ∗ (cid:19) (111)and o being the ( M × M )-dimensional zero matrix.In the following, the response amplitudes and coefficients are written in bra-ket notation for reasons thatbecome clear below. In order to render Eq. (92) an eigenvalue equation, it is helpful to notice that only theterm proportional to the external frequency ω has no projector, which implies PM | u (cid:105)| v (cid:105)| C u (cid:105)| C v (cid:105) = M | u (cid:105)| v (cid:105)| C u (cid:105)| C v (cid:105) (112)21nd thus one may add a redundant projector PL ⇒ PLP to the left term of the LHS of Eq. (92). Multiplyingnow with M − / from the left yields (cid:16) M − / PLPM − / − ω (cid:17) M + / | u (cid:105)| v (cid:105)| C u (cid:105)| C v (cid:105) = M + / PR [ ˆ f + , ˆ f − ] (113)or, in a more compact notation, (cid:0) ¯ L − ω (cid:1) | u (cid:105)| v (cid:105)| C u (cid:105)| C v (cid:105) = R [ ˆ f + , ˆ f − ] (114)where ¯ L = M − / PLPM − / = (cid:18) ρ oo − / P oo L oo P oo ρ oo − / ρ oo − / P oo L oc L co P oo ρ oo − / L cc (cid:19) , (115)as well as (cid:0) u , v , C u , C v (cid:1) t = M + / ( u , v , C u , C v ) t and R = M + / PR were introduced. To keep thenotation simple, the bar over the individual quantities in Eq. (114) will be omitted in the following. For thehomogeneous case, i.e., with R = 0, one arrives at the central equation of this work, given by( L − ω ) | u (cid:105)| v (cid:105)| C u (cid:105)| C v (cid:105) = 0 (116)which is a standard eigenvalue problem of the LR matrix L . Its solutions, i.e., the set of eigenvalues { ω k = E k − ε } and eigenvectors { ( | u k (cid:105) , | v k (cid:105) , | C uk (cid:105) , | C vk (cid:105) ) t } where k ∈ N labels the excitations, can be used to solvethe inhomogeneous problem of Eq. (114) by employing the ans¨atze | u (cid:105)| v (cid:105)| C u (cid:105)| C v (cid:105) = (cid:88) k c k | u k (cid:105)| v k (cid:105)| C uk (cid:105)| C vk (cid:105) (117)and R = − (cid:88) k γ k | u k (cid:105)| v k (cid:105)| C uk (cid:105)| C vk (cid:105) (118)where the minus sign is chosen arbitrarily. One finally obtains (cid:88) k c k ( ω k − ω ) | u k (cid:105)| v k (cid:105)| C uk (cid:105)| C vk (cid:105) = − (cid:88) k γ k | u k (cid:105)| v k (cid:105)| C uk (cid:105)| C vk (cid:105) (119)which means that the expansion coefficients c k in Eq. (117) are given by c k = γ k ω − ω k . (120)22t is stressed that the response weights { γ k } are the only quantities that depend on the perturbing fieldsˆ f + and ˆ f − and thus on the actual shape of the perturbation. All other quantities like the LR matrix L ,the eigenenergies { ω k } and the eigenvectors { ( | u k (cid:105) , | v k (cid:105) , | C uk (cid:105) , | C vk (cid:105) ) t } are obtained from the homogeneouseigenvalue problem of Eq. (116) where the amplitudes ˆ f + and ˆ f − of the perturbing field do not appear.Multiplying Eq. (118) by (cid:16) (cid:104) u k | , −(cid:104) v k | , (cid:104) C uk | , −(cid:104) C vk | (cid:17) from the left, taking into account the orthonor-mality conditions for the perturbed amplitudes and coefficients given by (cid:104) u k | u k (cid:48) (cid:105) − (cid:104) v k | v k (cid:48) (cid:105) + (cid:104) C uk | C uk (cid:48) (cid:105) − (cid:104) C vk | C vk (cid:48) (cid:105) = δ kk (cid:48) (121) (cid:104) v k | u k (cid:48) , ∗ (cid:105) − (cid:104) u k | v k (cid:48) , ∗ (cid:105) + (cid:104) C vk | C uk (cid:48) , ∗ (cid:105) − (cid:104) C uk | C vk (cid:48) , ∗ (cid:105) = 0 , (122)leads to the response weights γ k = − (cid:16) (cid:104) u k | , −(cid:104) v k | , (cid:104) C uk | , −(cid:104) C vk | (cid:17) R = − (cid:16) (cid:104) u k | , −(cid:104) v k | , (cid:104) C uk | , −(cid:104) C vk | (cid:17) − ( ρ ) / ˆ P ˆ f + | φ (cid:105) ( ρ , ∗ ) / ˆ P ∗ ˆ f − , ∗ | φ , ∗ (cid:105)− (cid:80) i,j (cid:104) φ i | ˆ f + | φ j (cid:105) ˆ b † i ˆ b j | C (cid:105) (cid:80) i,j (cid:104) φ j | ˆ f − , ∗ | φ i (cid:105) (cid:16) ˆ b † i ˆ b j (cid:17) ∗ | C , ∗ (cid:105) = (cid:104) u k | ˆ f + ( ρ ) / | φ (cid:105) + (cid:104) v k | ˆ f − , ∗ ( ρ , ∗ ) / | φ , ∗ (cid:105) + (cid:88) i,j (cid:104) φ i | ˆ f + | φ j (cid:105)(cid:104) C uk | ˆ b † i ˆ b j | C (cid:105) + (cid:88) i,j (cid:104) φ j | ˆ f − , ∗ | φ i (cid:105) (cid:68) C vk (cid:12)(cid:12)(cid:12)(cid:16) ˆ b † i ˆ b j (cid:17) ∗ (cid:12)(cid:12)(cid:12) C , ∗ (cid:69) (123)where the vector | φ (cid:105) = (cid:0) | φ (cid:105) , ..., | φ M (cid:105) (cid:1) t collects the stationary ground-state orbitals. Each response weight γ k quantifies how strong the k -th excitation contributes to the overall response of the system to an externalperturbation with amplitudes ˆ f + and ˆ f − . For instance, it may happen that for a chosen shape of theperturbing external fields, some excited states do not respond at all and thus appear to be transparent.The response density of the k -th excited state due to the perturbation can in general be expressed as∆ ρ k ( r ) = ∆ ρ ko ( r ) + ∆ ρ kc ( r ) (124)where ∆ ρ ko ( r ) and ∆ ρ kc ( r ) denote the orbitals’ and coefficients’ contribution, respectively. Those are given by∆ ρ ko = (cid:104) φ | ( ρ ) / (cid:0) | u k (cid:105) + | v k , ∗ (cid:105) (cid:1) + (cid:0) | u k (cid:105) + | v k , ∗ (cid:105) (cid:1) † ( ρ , ∗ ) − / ρ | φ (cid:105) (125)and ∆ ρ kc = M (cid:88) i,j =1 φ , ∗ i φ j (cid:16) (cid:104) C | ˆ b † i ˆ b j | C uk (cid:105) + (cid:104) C vk, ∗ | ˆ b † i ˆ b j | C (cid:105) (cid:17) . (126)To summarize, the many-body LR theory based on the MCTDHB working equations, termed LR-MCTDHB,was introduced and derived. For a given N -bosons system in a potential trap V and a two-body interactionˆ W , the central equation in order to obtain the excitation spectrum is Eq. (116). Although this equationappears to be rather straightforward to solve, it in general represents a demanding task because firstly, L isnot hermitian and secondly, it can become extremely large in size when the number of bosons N and/or thenumber of orbitals M grows. In Appendix B, it is shown that the LR matrix can be brought to block-diagonalform by using a complex transformation, which halves the dimensionality of the eigenvalue problem. However,it involves matrix-matrix products of the individual blocks of L , which can be numerically very expensive. Ithence remains a system-dependent question whether it is beneficial to perform this transformation or not.23n the next section, Section 4, details on the numerical implementation of LR-MCTDHB are presented.First, general rather technical challenges of the eigenvalue problem in Eq. (116) are discussed before explainingthe code structure and the numerical methods utilized. Finally, it is once again stressed that for M = 1, LR-MCTDHB reduces to the (particle-conserving) BdG equations, i.e., BdG ≡ LR-GP ≡ LR-MCTDHB(1). Thederivation of the BdG theory is described in Appendix A.1.2. Since the BdG is commonly used to calculateexcitations of trapped BECs, its results for certain applications (see Section 6) are compared to the many-bodyresults for
M > In this section, information about the LR-MCTDHB implementation are presented. Before dealing with thetechnical and numerical details in Sections 4.2 and 4.3, a brief overview about the general structure of thecode is given in Section 4.1, together with a description of technical challenges that were to be solved.
As a starting point, the structure of the numerical implementation, consisting of several thousand lines of code,is explained. The entire task of calculating the energies of excited states and the corresponding correctionamplitudes of a trapped interacting BEC, i.e., solving the eigenvalue problem of Eq. (116), is separated intotwo major parts. Those are (i) the construction of the LR matrix ¯ L (simply denoted by L in the following) asdescribed in Eq. (115) and (ii) finding a desired amount of eigenpairs, i.e., eigenvalues and the correspondingeigenvectors, of the low-energy spectrum. To this end, two running modes were implemented. Mode 1 is theconstruction mode. As an input, it requires a stationary state of a system of N bosons in a trap V withinteraction ˆ W , where the bosons can occupy M single-particle orbitals. Typically, this stationary state is theground state of the system, and it can therefore be obtained from a previous MCTDHB calculation usingimaginary time propagation. The output of mode 1 is the LR matrix L . Afterwards, the code can be restartedin mode 2, the diagonalization mode, which takes L as an input and calculates a desired amount of eigenpairsin the low-energy part of the excitation spectrum. Furthermore, the response weights and response densities,described in Section 3.2.2, are calculated.The main challenge of both the construction and the eigenvalue calculation is the memory consumptionof L , and, especially in mode 1, the size of intermediate matrices that appear during the construction, likethe projector P or the individual blocks L oo , L oc , L co and L cc . Although the memory consumption might bea minor issue for systems in 1D with small values of N and M , it can become a massive problem for largersystems. To give an example, one can consider a system with N = 100 particles and M = 4 orbitals. Insuch a case, already the coefficient matrix L cc has a dimensionality of 353702, which means it would needapproximately 4 terabytes of memory to store all matrix elements, assuming those are in general complexnumbers of 32 bytes. However, the available working memory (RAM) of compute nodes on modern highperformance computation clusters is typically much smaller, being of the order of 128 or 256 gigabytes. Onepossible solution of this problem is to store all (large) matrices that appear during the construction of L ,as well as L itself, in a sparse matrix storage format. There are multiple possibilities to do that, e.g., thecompressed-sparse-row (CSR) or compressed-sparse-column (CSC) format (see, e.g., Ref. [58] for additionalformats).To obtain eigenpairs of L by running the code in mode 2, one is in need of a numerical method thatis capable of dealing with sparse and large non-hermitian matrices. The method of choice is the implicitlyrestarted Arnoldi method (IRAM), implemented in the ARPACK numerical library [59]. In the currentimplementation, the parallel version, termed PARPACK, is utilized [60]. The IRAM is presented in detail inSection 4.3. The PARPACK routines have a significant advantage over eigenvalue routines from other librarieslike ScaLAPACK [61] because they do not require the matrix that one wishes to diagonalize as an input, noteven in any sparse format. It is only required to hand over the result of matrix-vector multiplications (matvecs)as input, and the actual calculation of these matvecs takes place outside of any PARPACK routine, meaning24hat the implementation of the matvecs is completely independent. It can be carried out separately, utilizingroutines from other numerical libraries, e.g., the Intel Math Kernel Library (MKL) [62], that calculate matvecsof large and sparse matrices in a very efficient manner. Moreover, calculating matvecs it is a process that canbe easily parallelized which further improves the code performance.Due to the large amount of operations in both running modes, especially for large N and M in two- orthree-dimensional space, sequential runs of the code would require very long runtimes. Therefore, a fullyMPI-parallelized implementation of the code was carried out for both running modes. Whereas this was acomparatively straightforward task for mode 2 – since only the matvecs needed to be parallelized – it turned outto be a demanding challenge with respect to mode 1. A parallelization scheme that is a compromise betweentaking care of the least possible memory consumption on the one hand, and distributing the workload amongall participating MPI processes (PEs) as equal as possible on the other hand, is employed. In Section 4.2, theoverall construction process is explained in more detail. In this section, details on the implementation of the LR matrix construction are presented. Figure 1 showsthe structure of the code and its main steps. As input, the code requires the stationary state to which the LRanalysis should be applied, usually the ground state of a trapped BEC. In particular, one needs to hand overthe set of orbitals { φ i : 1 ≤ i ≤ M } as well as the corresponding coefficients C . Afterwards, intermediatematrices are calculated which are essential to built up the individual blocks of L as described in Eq. (93).This includes the kinetic energy operator ˆ T kin , the direct and exchange interaction operators ˆ W sl and ˆ K sl , aswell as the projection operator P . All of these calculations are done in parallel, and the final result is storedin a sparse storage format on the master process PE 0.Then, the matrices L oo and L oc are calculated. Whereas L oo is a 2( M · N grid ) square matrix with N grid denoting the number of grid points, L oc is a 2( M · N grid × N conf ) matrix. It is stressed that constructingthese two matrices is usually the most time-consuming part of the entire calculation of L due to the largeamount of operations involved. On the one hand, for problems where the grid is large (e.g., in 2D and 3D)and M ≥
7, the construction of L oo represents the essential part of the calculations. On the other hand,if the number of particles is large such that N conf = O (10 ) or higher, the construction of L oc becomes themost demanding task. Additionally, both matrices need to be multiplied by ρ oo − / P oo from the left, whichfurther increases the amount of operations significantly. Because of this, both the construction of L oo and L oc are MPI-parallelized, as well as the subsequent multiplication with the metric and projector matrices. As forthe preliminary calculation described above, the final results of ( ρ oo − / P oo L oo ) and ( ρ oo − / P oo L oc ) arestored on PE 0 in a sparse storage format.Next, the generation of L starts. Therefore, the problem is separated into the construction of the upper andlower part. With respect to the upper blocks, this means that only the upper left block needs to be multipliedby P oo ρ oo − / from the right, whereas the upper right block is already fully calculated and can be writtento the output file directly. To efficiently calculate the upper-left block, a row-by-row parallelization scheme isemployed where each PE calculates a different row at a time. Therefore, PE 0 distributes all necessary matrixelements of ( ρ oo − / P oo L oo ) and P to the other PEs. Finally, all PEs send back their calculated rows (i.e.,the non-zero elements and their coordinates) to PE 0, which writes the results to the output file.Concerning the construction of the lower part, again a row-by-row parallelization scheme is utilized. Forthe lower-left block, each PE generates a different row at a time. However, in the case of L cc , all PEs workon the same row at a time. Upon finalization, the results are collected by PE 0 and written to the output file,again in a sparse storage format.This section is ended by briefly discussing the issue of small but non-zero matrix elements. As alreadystated, the matrices produced are typically sparse, i.e., most elements are zero. Moreover, many other elementscan be very small but nonzero. Their impact on the subsequent eigenvalue calculation, which is described inthe next sections, is very often negligible. However, for a problem where the dimensionality of L is very large,those elements can become challenging in terms of memory consumption. Thus, it is important to define a25 nput: { φ i } , C Preliminary calculations: ˆ T kin ˆ W sl ˆ K sl P • MPI-parallelized • On exit stored on PE 0
Calculate L oo L oc to do matrix multiplications ρ oo − / P oo L oo ρ oo − / P oo L oc • MPI-parallelized • On exit stored on PE 0
Row-by-row generationof upper blocks of L • Each PE generates differentrow • Written to file by PE 0
Row-by-row genera-tion of lower blocks of L • For lower left block: EachPE generates different row • For L cc : All PEs share workof each row • Written to file by PE 0
Output: L Figure 1: Schematic diagram of the MPI-parallelized implementation for the construction of L . The parallelization scheme shares the amountof calculations equally between all PEs. Before L is generated, preliminary calculations of importantintermediate matrices are performed. To minimizethe memory consumption, those are only stored onPE 0 in a sparse storage format. For the row-by-rowconstruction of the upper and lower blocks of L , PE0 distributes all necessary information to the otherPEs, collects the results afterwards and writes themto a file. See text for details. τ that is used to minimize the amount of non-zero elements which need to be stored. In thecurrent implementation, a matrix element is treated as zero and is therefore not stored if its magnitude issmaller than τ . For very large matrices, the value of τ has to be defined such that on the one hand not toomany matrix elements are set to zero in order to obtain numerically accurate results for the eigenvalues andeigenvectors, but, on the other hand, it has to be large enough because too many very small elements cansignificantly decrease the performance of the numerical method for calculating eigenpairs of L . The latternumerical procedure will be explained in the next section. In order to calculate the lowest eigenvalues of a given LR matrix, the IRAM is employed. In the following,the main ideas of this approach are presented by first discussing the basic ingredients (the Arnoldi method,the QR-iteration, and polynomial filtering) before dealing with the IRAM explicitly. In short, it generatesa Krylov subspace of a given matrix whose basis vectors are used to build a much smaller matrix. Theeigenvalues of this matrix are also approximate eigenvalues of the original matrix. After each iteration of theIRAM, the convergence of a set of wanted eigenvalues is checked, and if it is still not converged, the entireprocess is restarted with an initial vector where the directions of unwanted eigenvalues are projected out.Thus, the convergence with respect to the wanted eigenvalues is gradually enhanced. Below, the individualsteps of the IRAM are explained in more detail.
The first step is to calculate an orthogonal projection of a sparse, complex, in general non-hermitian matrix A ∈ C n × n . Therefore, the Arnoldi algorithm is used which generates an orthonormal basis (ONB) of the k -dimensional subspace given by K k ( A, v ) = span( v, Av, ..., A k − v ) (127)where v ∈ C n is an initial vector and k ≤ n , typically even k (cid:28) n . K k is called the k -dimensional Krylovsubspace of A . The entire Arnoldi algorithm is given by Algorithm 1. The ONB is obtained by performing aGram-Schmidt-orthonormalization at each Arnoldi step, i.e., making the vector v j in the j -th step orthogonalto all previous basis vectors v i , i ≤ j, and normalizing it afterwards. The elements of the ONB are called theArnoldi vectors. If the resulting Arnoldi vector v j +1 of the j -th step is non-zero, its norm is stored in thematrix element h j +1 ,j of the matrix ˜ H k +1 ∈ C ( k +1) × k . If j ≤ k , the next Arnoldi step is initialized. After thestep for j = k , the iteration stops and returns the k -dimensional ONB of K k , plus an additional orthonormalvector v k +1 ∈ C n , together with the full matrix ˜ H k +1 . It is possible that for a certain step j < k , the normof the Arnoldi vector v j +1 is zero, meaning that one has found an ONB of an invariant Krylov subspace of A with dimensionality j . In this case, the Arnoldi algorithm finishes immediately.As an essential result, one obtains the Arnoldi decomposition of A given by AV k = V k H k + h k +1 ,k v k +1 e tk (128)where V k ∈ C n × k contains the Arnoldi basis as columns and H k ∈ C k × k is an upper Hessenberg matrix, i.e., { ( H k ) ij } = 0 if i > j + 1, built by the first k rows of ˜ H k +1 . The second term in Eq. (128), which containsthe matrix element h k +1 ,k , the additional Arnoldi vector v k +1 and the transpose of the k -th unit vector, i.e., e k = (0 , , ..., , t ∈ C k , takes care of the fact that the orthonormal projection of A , A ≈ V k H k V † k , (129)is in general not exact. However, for the case of having found an invariant Krylov subspace of A of lowerdimensionality j ≤ k , it becomes exact. It is straightforward to realize that the eigenvalues of H k are27 lgorithm 1 The Arnoldi algorithm v ← v, || v || = 1 for j = 1 , j ≤ k do z ← Av j for i = 1 , i ≤ j do h ij ← (cid:104) v i , z (cid:105) z ← z − h ij v i end for h j +1 ,j ← || z || if h j +1 ,j = 0 then quit else v j +1 ← z/h j +1 ,j end if end for approximate, or, in the case of an invariant subspace, exact eigenvalues of A. This is due to the cyclicpermutation symmetry of the determinant that defines the characteristic polynomial of A , i.e.,0 = det( A − λI ) Eq. (129) ≈ det( V k H k V † k − λI )= det( V k [ H k − λI ] V † k )= det( V † k V k (cid:124) (cid:123)(cid:122) (cid:125) = I [ H k − λI ])= det( I ) (cid:124) (cid:123)(cid:122) (cid:125) =1 det( H k − λI )= det( H k − λI ) (130)with λ ∈ C and I denoting the identity matrix.If u i ∈ C k is an eigenvector of H k and θ i ∈ C the corresponding eigenvalue, the pair ( θ i , V k u i ) is called aRitz pair of A with Ritz value θ i and the associated Ritz vector V k u i ∈ C n . In order to estimate how accuratethis Ritz pair approximates an exact eigenpair of A , one can calculate the residual norm r defined as r ≡ || A ( V k u i ) − θ i ( V k u i ) || = || ( V k H k + h k +1 ,k v k +1 e tk ) u i − θ i ( V k u i ) || = || θ i ( V k u i ) + h k +1 ,k v k +1 e tk u i − θ i ( V k u i ) || = | h k +1 ,k | · | v k +1 | (cid:124) (cid:123)(cid:122) (cid:125) =1 ·| e tk u i | = | h k +1 ,k | · | e tk u i | , (131)i.e., r is the product of the magnitudes of the matrix element h k +1 ,k of ˜ H k +1 and the k -th element of u i . Thesmaller it is, the better does the Ritz pair ( θ i , V k u i ) approximate the exact i -th eigenpair of A .To summarize, the eigenvalue problem of a typically large matrix A can be mapped onto an eigenvalueproblem of a much smaller matrix H k whose eigenvalues are believed to be good approximations of theeigenvalues of A . Due to the Hessenberg form of H k , its eigenvalues can efficiently be calculated by utilizingthe QR-iteration, which is explained in the next section. Additional details on the Arnoldi algorithm can befound in Ref. [63]. 28 .3.2 The QR-iteration The QR-iteration is a technique to calculate the eigenvalues and eigenvectors of a given matrix A ∈ C n × n . Itis based on the so-called QR-decomposition which reads A = QR, Q, R ∈ C n × n (132)where Q is a unitary matrix, i.e., QQ † = I , and R is upper triangular. This decomposition exists for anymatrix [63]. Plugging the ansatz of Eq. (132) into the characteristic polynomial gives0 = det( A − λI )= det( QR − λI )= det( QQ † ( QR ) QQ † − QλIQ † )= det( Q ( RQ − λI ) Q † )= det( RQ − λI ) (133)which means that the eigenvalues of A are the same as the eigenvalues of the matrix RQ . Furthermore, A = QR ⇐⇒ Q † AQ = RQ, (134)i.e., RQ is just a unitary transformation of A .The QR-iteration is described in Algorithm 2. At each step k , the QR-decomposition of the current matrix A k is used to define the matrix of the next step as A k +1 = R k Q k = Q † k A k Q k . It can be shown that thissequence converges to a triangular matrix, i.e., lim k →∞ A k = ˜ A with { ˜ A ij } = 0 if i > j , and the eigenvaluesof a triangular matrix are simply given by its diagonal elements. One can show that the eigenvalues of A k arethe same as the ones of A by considering again the characteristic polynomial0 = det( A k − λI )= det( Q † k − A k − Q k − − λI )= det( Q † k − Q † k − A k − Q k − Q k − − λI )= ... = det( Q † k − Q † k − ...Q † A (cid:124)(cid:123)(cid:122)(cid:125) = A Q ...Q k − Q k − − λI )= det( A − λI ) (135)where the cyclic permutation symmetry of the determinant is used in the last step.Although the QR-iteration possibly requires a large amount of operations per step for a general matrix A ,it turns out to be of very low computational cost with respect to upper Hessenberg matrices like the matrix H k obtained from a k -step Arnoldi process described in Algorithm 1. In such a case, one can apply ( k − G i , i ∈ { , ..., k − } , which are all unitary operations such that the overall QR-decompositionof H k at each QR-step j is given by ( H k ) j = Q j R j = ( G k − ) j ... ( G ) j R j . (136)Further details on Gives rotations can be found in Ref. [64]. With respect to the IRAM, the QR-iterationturns out to be important not only because it is used to calculate the eigenvalues of H k , but also in terms ofrestarting the Arnoldi iteration with an initial vector that incorporates shifts that filter out the direction ofunwanted eigenvectors of H k . Those shifts are obtained by utilizing QR-steps, which are explained in detailin Section 4.3.4. In Ref. [63], further details on the QR-iteration are presented.29 lgorithm 2 The QR-iteration A ← A Calculate the QR-decomposition A = Q R for k = 1 until convergence do A k ← R k − Q k − if A k ”triangular enough” then quit else Calculate the QR-decomposition A k = Q k R k end if end for4.3.3 Polynomial filtering Before dealing with the IRAM in particular, a brief discussion on polynomial filtering techniques of the initialvector v of a Krylov subspace iteration, like the Arnoldi algorithm of Section 4.3.1, is realized. In general,the target of applying a filter to a vector v is to make the contribution of unwanted directions very small, oreven zero.To understand how polynomial filtering works in general, a square matrix A ∈ C n × n whose eigenpairs aregiven by the set { ( λ i , v i ) | ≤ i ≤ n } with λ i ∈ C and v i ∈ C n is considered. Suppose the eigenvectors form abasis. Each arbitrary vector v ∈ C n can then be written as v = (cid:80) nj =1 c j v j with complex expansion coefficients { c j } . Assuming that one wants to filter out the contribution of the k -th eigenvector v k , one can define p ( t ) ≡ t − λ k (137)and apply this very simple polynomial onto v : p ( A ) v = n (cid:88) j =1 c j p ( A ) v j = n (cid:88) j =1 c j ( A − λ k ) v j = n (cid:88) j =1 c j ( λ j − λ k ) v j = n (cid:88) j (cid:54) = k c j ( λ j − λ k ) v j . (138)If A is hermitian, the overlap with v k vanishes, i.e., (cid:104) v k | p ( A ) v (cid:105) = n (cid:88) j (cid:54) = k c j ( λ j − λ k ) (cid:104) v k | v j (cid:105) (cid:124) (cid:123)(cid:122) (cid:125) = δ kj = 0 , (139)because in that case the eigenvectors of A would be orthogonal. Thus, the vector p ( A ) v is orthogonal to v k ,and any additional power iteration of this vector with respect to A preserves this property, meaning that (cid:104) v k | A n p ( A ) v (cid:105) = n (cid:88) j (cid:54) = k c j λ nj ( λ j − λ k ) (cid:104) v k | v j (cid:105) (cid:124) (cid:123)(cid:122) (cid:125) = δ kj = 0 . (140)Therefore, the Arnoldi iteration, applied to an initial vector p ( A ) v , yields a Krylov basis which is orthogonalto v k , and also the orthogonal projection H k would not contain λ k in its spectrum of eigenvalues. If A is30 lgorithm 3 The implicitly restarted Arnoldi method (IRAM) Take an initial vector v ← v , || v || = 1 Perform a k -step Arnoldi iteration to get AV k = V k H k + f k e tk with f k := h k +1 ,k v k +1 Calculate the eigenvalues θ , ..., θ k of H k by using the QR-iteration while m wanted eigenvalues of H k not converged do Select p = k − m unwanted eigenvalues θ , ..., θ p Perform p QR-steps with unwanted eigenvalues as shifts:[ H ( p ) k , Q ] := QR [ H k , θ , ..., θ p ] , V ( p ) k := V k Q Set f m := h k +1 ,k Q kp v k +1 + h ( p ) m +1 ,m v ( p ) m +1 Define H m := H ( p ) k (1 : m, m ), V m := V ( p ) k (1 : n, m ) Perform p additional Arnoldi steps onto V m , H m and f m to obtaina k -step Arnoldi decomposition Calculate the eigenvalues of the new H k by using the QR-iteration end while non-hermitian, the overlap of p ( A ) v and v k is not strictly zero since the eigenvectors { v i } can be linearlydependent. However, within the IRAM, the overall contribution of v k gradually decreases upon increasingthe number of restarts with properly shifted initial vectors. Nevertheless, this shows that the LR-matrix L ,which is in general non-hermitian, leads to additional numerical effort which would be absent in the case ofa hermitian matrix.It is also possible to do multiple shifts at once by using a filter given by p ( t ) = N (cid:89) i =1 (cid:0) t − λ π ( i ) (cid:1) (141)where N ≤ n is the number of shifts and π ( i ) ∈ { , ..., n } ∀ i . In this case, the directions of N differenteigenvectors of A are filtered out. The above choices of p ( t ) do only reflect the simplest possibilities of filterpolynomials, but also more sophisticated ones like the Chebyshev polynomials can be chosen [64]. After the introduction of the three essential ingredients of the IRAM, namely the Arnoldi algorithm, theQR-iteration, and polynomial filtering, one can finally discuss the IRAM itself in detail. The full algorithm isgiven in Algorithm 3. At first, a k -step Arnoldi iteration is performed, yielding an Arnoldi decomposition of A as described in Eq. (128). Then, the eigenvalues ( θ , ..., θ k ) are calculated utilizing the QR-iteration. Fromthis set, p unwanted eigenvalues ( θ , ..., θ p ) are selected. If one for example seeks the eigenvalues of A thathave the smallest magnitude, an appropriate selection would be the p eigenvalues of H k that have the largestmagnitude. Afterwards, a sequence of p QR-shifts with the unwanted eigenvalues is performed. For the firstshift, the Arnoldi decomposition reads( A − θ I ) V k = V k ( H k − θ I ) + h k +1 ,k v k +1 e tk (142) ⇐⇒ ( A − θ I ) V k = V k Q R + h k +1 ,k v k +1 e tk (143) ⇐⇒ AV k Q = V k Q ( R Q + θ I ) + h k +1 ,k v k +1 e tk Q (144)where the QR-decomposition H k − θ I = Q R is used in Eq. (143). Setting H (1) k := R Q + θ I, V (1) k := V k Q , (cid:16) e (1) k (cid:17) t := e tk Q (145)31he next QR-shift can be applied:( A − θ I ) V (1) k = V (1) k ( H (1) k − θ I ) + h k +1 ,k v k +1 (cid:16) e (1) k (cid:17) t (146) H (1) k − θ I = Q R ⇐⇒ AV (1) k Q = V (1) k Q ( R Q + θ I ) + h k +1 ,k v k +1 (cid:16) e (1) k (cid:17) t Q (147)where one can define, similar to Eq. (145), H (2) k := R Q + θ I, V (2) k := V (1) k Q , (cid:16) e (2) k (cid:17) t := (cid:16) e (1) k (cid:17) t Q . (148)This procedure is continued until all p QR-shifts are applied.Multiplying Eq. (143) by e yields( A − θ I ) V k e = V k Q R e + h k +1 ,k v k +1 e tk e ⇐⇒ ( A − θ I ) v = V k Q r e ⇐⇒ ( A − θ I ) v = v (1)1 r (149)where r is the upper left element of R and v (1)1 is the first column of V k Q . This in turn means that thefirst vector of the new basis, i.e., v (1)1 , is proportional to ( A − θ I ) v . As discussed in the previous section,this polynomial filter reduces the contribution of the unwanted eigenvector corresponding to θ in v (1)1 , or,in the optimal case, makes its contribution even vanish completely. This holds for all successive vectors in V (1) k . After having applied the second QR-shift, one can show that v (2)1 ∼ ( A − θ I )( A − θ I ) v . Continuingin the same manner until the p -th QR-shift is applied yields v ( p )1 ∼ (cid:81) pi =1 ( A − θ i I ) v . Then, by taking thefirst m = k − p columns of V ( p ) k as the new Arnoldi basis V m , as well as taking the leading principal ( m × m )-matrix of H ( p ) k as the new orthogonal projection H m of A , one obtains an m -step Arnoldi decomposition AV m = V m H m + f m e tm with an error term given by f m = h k +1 ,k Q kp v k +1 + h ( p ) m +1 ,m v ( p ) m +1 (150)where it is used that the unitary matrix Q is of Hessenberg form. Based on this, one performs p additionalArnoldi steps yielding again a k -step Arnoldi decomposition. The essential difference to the previous k -stepArnoldi decomposition is that now all basis vectors in V k are approximately or even fully orthogonal to theeigenvectors of the associated unwanted eigenvalues. The accuracy of the eigenvalues and eigenvectors ofthe newly computed matrix H k can be estimated with the residual norm r introduced in Section 4.3.1, Eq.(131). If it is below a predefined tolerance for a desired amount of wanted eigenpairs, the algorithm finishes.Otherwise, the set of eigenvalues of H k is again separated into wanted and unwanted eigenvalues, and theprocedure of QR-shifting and refining the given Arnoldi basis restarts.By employing the above described shift strategy without restarting the Arnoldi iteration from the verybeginning, one obtains the same result as if an explicit restart of the entire k -step Arnoldi process with theinitial vector v ( p )1 had been performed. Because of this similarity, the overall algorithm is said to be implicitly restarted. The computational effort in terms of matvecs with A is given by ( k + p ) operations in both cases,i.e., for explicit and implicit restarting. However, implicit restarting turned out to be numerically more stable,and has additional advantages which are beyond the scope of this work.Figure 2 shows the general structure of the code for the matrix diagonalization utilizing the implementationof IRAM in the PARPACK library [60]. Essential for the performance is the efficient parallel implementationof the matvecs which are done outside of the PARPACK routine pznaupd . The latter carries out the IRAMiterations in parallel [60]. It is important to note that the IRAM, although it is particularly efficient for sparsematrices, would also work for dense matrices. Further details are presented in Refs. [64, 59, 63].32o summarize, the newly developed implementation of LR-MCTDHB was presented in this section andthe code was explained in detail. Furthermore, the IRAM, a very powerful method for computing eigenpairsof very large and sparse matrices, was described. In the next section, the code is benchmarked against anexactly-solvable model. Input: L Distribution of L • Each PE gets the sameamount of rows from theupper and lower blocks
Define initial vector v Arnoldi iteration • Call pznaupd (PARPACK routine thatperforms IRAM iterations)
Eigenvaluesconverged? Matvec z := Av i Calculate eigenvec-tors, response weightsand response densities Output noyes
Figure 2: Schematic diagram of the MPI-parallelized implementation for the diagonalizationof L . Different chunks of non-zero elements of L aredistributed among the PEs. The Arnoldi iterationis carried out by the routine pznaupd of ParallelARPACK (PARPACK) [60], which either finishesupon convergence of a desired amount of eigenvaluesor expects the user to perform a matvec and returnthe result to it. Finally, the eigenvalues, eigenvec-tors, as well as response weights and response den-sities are calculated and written to output files. Seetext for details. Comparison with an exactly-solvable model
After the presentation of the theoretical and numerical frameworks utilized, along with the structure of theimplemented code, the newly developed implementation of LR-MCTDHB is benchmarked in the following.Furthermore, since it builds the basis for the utilized LR theory, results from a benchmark of the MCTDHBimplementation used in this work [65], which makes use of a sophisticated mapping of bosonic operators inFock space [66], are presented as well. This implementation is now also available with a recently developedgraphical user interface and works on all common operating systems [67]. Other codes became availableduring the last years, e.g., the MCTDH-X software package [68] which contains a separate graphical userinterface. Furthermore, the multi-layer (ML-)MCTDHB approach [69, 70, 71], which is particularly suited forthe description of Bose-Bose or Bose-Fermi mixtures, is now implemented.As a reference system, the exactly-solvable harmonic-interaction model (HIM) is used. In the latter model,both the single-particle and interaction potentials are of harmonic type. This model system has been usedrecently to benchmark the fermionic counterpart of MCTDHB, termed MCTDHF [221], as well as the time-dependent restricted-active-space self-consistent-field theory for bosons (TD-RASSCF-B) [222] and its variantfor bosonic mixtures [223]. In the first part of this section, an introduction to the HIM and the analytic resultsof its energy spectrum are given (Section 5.1). Afterwards, an overview of the MCTDHB benchmark againstthe HIM, which has been carried out in great detail in Ref. [33], is presented (Section 5.2). This includes boththe ground state and the out-of-equilibrium quantum dynamics of a bosonic many-particle system. Finally, acomparison of the numerical results of the low-energy spectrum obtained from LR-MCTDHB and the analyticvalues, both for the 1D and 2D cases, is made (Section 5.3).
The HIM Hamiltonian contains the harmonic trap potential V ( r ) = 12 (cid:0) Ω x x + Ω y y + Ω z z (cid:1) , r = ( x, y, z ) t (151)where the boson mass m is set to unity. In D > { Ω i } are equal, whereas the latterrefers to the case where at least one frequency differs from the others.The interaction between the i -th and j -th boson is described asˆ W ( r i , r j ) = λ | r i − r j | (152)where a positive (negative) interaction strength λ denotes attraction (repulsion) between the bosons.As described in Ref. [72], the following coordinate transformation yields the separation of the relative andcenter-of-mass (c.m.) coordinates, Q k = 1 (cid:112) k ( k + 1) k (cid:88) i =1 ( r k +1 − r i ) , ≤ k ≤ N − Q N = 1 √ N k (cid:88) i =1 r i , (154)34here N is the number of bosons in the system. This leads to the Hamiltonian in the c.m. frame given byˆ H = ˆ H rel + ˆ H c.m. (155)ˆ H rel = 12 (cid:88) j = x,y,z N − (cid:88) k =1 ( p j,k + δ j Q j,k ) (156)ˆ H c.m. = 12 (cid:88) j = x,y,z ( p j,N + Ω j Q j,N ) (157)with momenta p j,N = (1 /i ) ∂ Q j,N of the c.m. coordinates, p j,k = (1 /i ) ∂ Q j,k of the relative coordinates, andthe parameters δ j = Ω j ± N | λ | where again the plus (minus) sign denotes attraction (repulsion) betweenthe bosons. For simplicity (cid:126) = 1 is assumed. It is instructive to interpret the Hamiltonian as a compositemodel system of ( N −
1) identical non-interacting particles in a harmonic confinement with trap frequencies δ = ( δ x , δ y , δ z ) t and the c.m. particle which moves in the harmonic confinement with the original trapfrequencies { Ω i } .The ground-state wave function in the c.m. frame takes the formΨ ( Q , ..., Q N ) = (cid:89) j = x,y,z (cid:18) δ j π (cid:19) N − (cid:18) Ω j π (cid:19) exp (cid:32) − δ j N − (cid:88) k =1 Q j,k − Ω j Q j,N (cid:33) , (158)i.e., it is basically a product of N harmonic oscillator (HO) ground-state wave functions in each spatial direction[72]. However, in the laboratory frame, the representation of the ground state is much more involved. Alreadyfor N = 2 bosons in 1D, it is an infinite sum of Hartree products of HO eigenfunctions instead of a singleproduct of two HO ground states in the c.m. frame [33]. Numerical convergence towards this state is thereforea very demanding challenge in the laboratory frame.The exact solution for the energy levels reads E ( { n j } , { m j } ) = (cid:88) j = x,y,z (cid:26)(cid:18) n j + N − (cid:19) δ j + (cid:18) m j + 12 (cid:19) Ω j (cid:27) (159)where the { n j } denote the quantum numbers of excitations of the relative coordinates and the { m j } denotethe quantum numbers of c.m. excitations in the x -, y - and z -directions. In particular, this means that theground-state energy is given by E ≡ ε = (cid:88) j = x,y,z (cid:26) ( N − δ j + 12 Ω j (cid:27) (160)whereas the energy distance between any excited state and the ground state reads ω ( { n j } , { m j } ) = (cid:88) j = x,y,z { n j δ j + Ω j m j } . (161)Eqs. (160) and (161) will be used for the comparison to the numerical results obtained from (LR-)MCTDHB.The degeneracies of excited states is discussed in [73]. It is worth noting that for any number of bosons N , there is no solution for an excitation corresponding to (cid:80) j = x,y,z n j = 1, i.e., a single-particle excitationwith respect to the relative coordinates. This is due to the even permutation symmetry for identical bosons,Eq. (5), and it can be shown by mathematical induction as follows. For the sake of simplicity, the proofis restricted to the 1D case, i.e., Q i = Q i , which does not limit the generality since in higher dimensions,the Hamiltonian decouples into the x -, y - and z -components, see Eqs. (156) and (157). The wave functioncorresponding to a first order excitation of the relative coordinates, denoted in the following by Ψ (cid:48) , obeysΨ (cid:48) ( Q , ..., Q N ) ∼ ∂∂Q k Ψ ( Q , ..., Q N ) , ≤ k ≤ N − . (162)35or N = 2, this yields Ψ (cid:48) ( Q , Q ) ∼ Q Ψ ∼ ( x − x )Ψ , (163)see Eq. (158). Because the wave function needs to be fully symmetric under the exchange x ↔ x , one findsΨ (cid:48) ∼ { ( x − x ) + ( x − x ) } Ψ = 0 , (164)meaning that there is no first order relative excited state for two bosons. Assuming that this is also true for N bosons, i.e., Ψ (cid:48) ( Q , ..., Q N ) ∼ ∂∂Q k Ψ ( Q , ..., Q N ) = 0 ∀ ≤ k ≤ N − , (165)one obtains for ( N + 1) bosons and K ≤ N Ψ (cid:48) ( Q , ..., Q N +1 ) ∼ ∂∂Q K Ψ ( Q , ..., Q N +1 ) ∼ K (cid:88) i =1 ( x K +1 − x i )Ψ = (cid:32) Kx K +1 − K (cid:88) i =1 x i (cid:33) Ψ = (cid:32) Kx K +1 − x K − K − (cid:88) i =1 x i (cid:33) Ψ = Kx K +1 − x K − ( K − x K + ( K − x K − K − (cid:88) i =1 x i (cid:124) (cid:123)(cid:122) (cid:125) ∼ Q K − Ψ ∼ K ( x K +1 − x K ) Ψ + Q K − Ψ (166)where the argument of Ψ is suppressed from the second line onward. For the second term in Eq. (166),one can make use of the assumption from Eq. (165) that the contribution of Q K − Ψ to Ψ (cid:48) vanishes aftersymmetrization of the wave function. Moreover, the first term is proportional to the relative coordinate of onlytwo bosons, for which it was already shown that its contribution to Ψ (cid:48) vanishes as well after symmetrization.Thus, one infers that also for ( N + 1) bosonsΨ (cid:48) ( Q , ..., Q N +1 ) = 0 (167)which ends the proof. At first, numerical results obtained from MCTDHB for the ground state of the 1D HIM for different particlenumbers and the exact values from Eq. (160) are compared. The index ’ x ’ is suppressed in the followingfor simplicity. The trap frequency is set to Ω = 1 and the two-particle interaction strength λ is adjustedsuch that the mean-field interaction parameter Λ = λ ( N −
1) is kept fixed at Λ = 0 .
5. Results are shown inTable 1, which originally appeared in Ref. [33] where it has been demonstrated that the MCTDHB methodis much more efficient and accurate for the description of the ground state and dynamics of a many-bosonsystem than the widely used FCI method. One observes that MCTDHB gives highly accurate results for the36 N = 10 N = 100 N = 10001 7.071067811865483 70.71067811865483 707.10678118654832 7.038769026303168 70.68016951747168 707.07643342573153 7.038350652406389 70.68012541218675 707.07642898718654 7.038348424909910 70.680125391745495 7.038348415349058 70.680125391737626 7.0383484153114947 7.038348415311018 ε exact Table 1: Ground-state energies of N = 10 , , and 1000 bosons in the 1D HIM with interaction parameter Λ = 0 . .
0. Numerical results are shown for different numbers of orbitals M . Exact results are obtained from Eq. (160).The table is adapted from Ref. [33]. ground-state energy for all considered particle numbers. Furthermore, one needs less self-consistent orbitals M in order to converge to the exact ground-state energy when N increases. This is anticipated, since in theHartree limit where N → ∞ for fixed Λ, the energy per particle converges to the GP value. What is howeversurprising is the degree of accuracy. As an example, for N = 1000 the relative error ( ε M =3 − ε ) /ε is lower than 10 − %. The fact that MCTDHB can very accurately assemble the ground state of a largemany-boson system with only very few orbitals has far-reaching consequences for the accuracy of excitationenergies obtained by LR-MCTDHB. As discussed in Section 3.2.1, the latter will be closer to the exact valuesif the ground state to which LR is applied is very accurately described. Further examples for the ground-stateenergies in the 1D and 2D HIM will be given in Section 5.3 and in Appendix C, where in particular theconvergence of excitation energies is discussed.With regard to the dynamics, comparisons of numerical and exact results were carried out in Ref. [33] for (i)the out-of-equilibrium quench dynamics and (ii) the dynamics due to time-dependent driving of the trappingfrequency and the interaction strength. Concerning (i), Fig. 3 compares the exact, FCI and MCTDHB resultsfor the oscillations of the density ρ of N = 2 bosons after a sudden quench of the interaction parameter fromΛ = 0 to Λ = 0 .
5. The initial state Ψ is therefore the ground state of the non-interacting system. In principle,the time-evolution is given by | Ψ( t ) (cid:105) = ∞ (cid:88) n =0 c j e − iE j t | α j (cid:105) , (168) c j = (cid:104) α j | Ψ(0) (cid:105) (169)where the {| α j (cid:105)} are an arbitrary basis of the two-particle Hilbert space and the { E j } the correspondingeigenenergies. One possible choice of this basis are the solutions of the time-independent problem. It remainsto compute the overlap integrals in Eq. (169). In practice, the infinite sum in Eq. (168) has to be truncated,and for the case described in Fig. 3, numerical convergence at each time step of the expansion in Eq. (168)is ensured by including the 60 lowest energy eigenstates of the time-independent problem. However, utilizingonly M = 4 time-adaptive orbitals, the MCTDHB calculations lead to a density oscillation that can no longerbe distinguished from the exact result. Already for M = 3 the results are highly accurate. If a fixed set oforbitals is used (FCI), even 8 orbitals do not lead to the same accuracy than MCTDHB(4). This remarkablyshows the numerical benefit of using a time-adaptive, self-consistent basis set.Fig. 4 shows the time-evolution of the c.m. energy (cid:15) ( t ) when both the trapping frequency and thetwo-body interaction strength are time-dependent and given by Ω( t ) = Ω (cid:112) f ( t ) with a constant Ω and λ ( t ) = λ (cid:104) − Ω Nλ f ( t ) (cid:105) . The function f ( t ) denotes a time-dependent driving frequency. In that way, the37rapping frequency of the ( N −
1) relative coordinates remains time-independent, i.e., δ = (cid:115) Ω [1 − f ( t )] + 2 N λ (cid:20) N λ f ( t ) (cid:21) = (cid:113) Ω + 2 N λ . (170)Thus, although in the laboratory frame the bosons move in a time-dependent trap and interact via a time-dependent potential, the same system in the c.m. frame appears to be much less complicated because only thetrapping potential of the c.m. particle is time-dependent due to f ( t ). Because the c.m. and relative motionsseparate, the time-dependent energy of the system is given by E = N − δ + ε ( t ) where ε ( t ) is determinedby the energy expectation value ε ( t ) = (cid:104) Ψ N ( t ) | ˆ H c.m. ( t ) | Ψ N ( t ) (cid:105) (171)of the c.m. particle. The state | Ψ N ( t ) (cid:105) is obtained from the effective single-particle Schr¨odinger equation i∂ t | Ψ N ( t ) (cid:105) = ˆ H c.m. ( t ) | Ψ N ( t ) (cid:105) . From Fig. 4, one can deduce that MCTDHB is capable to account for thedynamical evolution of ε ( t ) very accurately, even for longer times. Furthermore, one observes that for ahigher number of particles N a smaller number of time-adaptive orbitals M is required to obtain the samelevel of numerical accuracy. In contrast to that, the GP evolution of ε ( t ) only follows the exact curves of N = 10 and 50 bosons for a very short time up to t ≈
1, and becomes remarkably inaccurate afterwards. Forthe long-time behavior, it cannot even reproduce qualitatively the profile of ε ( t ).To sum up, the previous examples clearly indicate that the MCTDHB theory qualifies for describing theground state, the dynamics due to a sudden quench of system parameters as well as the dynamics due totime-dependent trapping and interaction potentials, and that it is clearly superior to the FCI many-body and Figure 3: Out-of-equilibrium dynamics of the one-particle density ρ ( x = 0) for N = 2 bosons after an interaction quench fromΛ = 0 to Λ = 0 .
5. Results are shown for different numbers of orbitals M . The trap frequency is Ω = 1 .
0. For the time scaleshown, the MCTDHB(4) dynamics cannot be distinguished from the exact dynamics. The FCI method is less accurate even for8 orbitals. All quantities are dimensionless. See text for details. The figure is taken from Ref. [33]. igure 4: Dynamics in the HIM due to a time-dependent driving f ( t ) of the trapping frequency Ω( t ) = Ω (cid:112) f ( t ) and theinteraction strength λ ( t ) = λ (cid:20) − Ω Nλ f ( t ) (cid:21) . The upper panel shows the complicated time-dependent profile of the drivingfrequency f ( t ). The lower panel shows the evolution of the c.m. energy ε ( t ) for different levels of MCTDHB( M ). For N = 10bosons, the curve of MCTDHB(7) can hardly be distinguished from the exact result for the shown time interval. The same holdstrue for M = 6 and N = 50 bosons. Even for less orbitals, MCTDHB yields an accurate description of the evolution of ε ( t ) upto t ≈
35 for N = 10 , M = 6 and t ≈
25 for N = 50 , M = 5. On the contrary, the GP results do only follow the exact curves forvery short times and clearly fail to describe the long-time behavior, even qualitatively. All quantities are dimensionless. See textfor details. The figure is taken from Ref. [33]. GP mean-field methods. The full benchmark, together with additional valuable details, can be found in Ref.[33].
In this section, the applicability of LR-MCTDHB to obtain highly accurate results for the energies of excitedstates of trapped multi-boson systems is discussed. As for the benchmark concerning the ground state anddynamics, the analytic results of the HIM from Eq. (161) are used to benchmark the numerical results. Inthe following, LR-MCTDHB is first benchmarked against the 1D HIM for repulsive bosons, i.e., with λ < x = Ω y is considered. In Appendix C, further benchmarksagainst the 1D HIM with attractive bosons, against the anisotropic case in 2D, as well as against the isotropicHIM in the rotating frame of reference are presented.For the benchmark against the 1D repulsive HIM, the system parameters are set to Ω = 1 for the trapfrequency and λ = − .
01 for the interaction strength, which yields a slight depletion of f ≈ .
03% for N = 10bosons and M = 4 orbitals. Calculations were carried out on a grid from [ − ,
9) with 128 grid points. Resultsfor the excitation energies relative to the ground-state energy, ω i = E i − ε , are shown in Table 2. There areseveral comments to be made. At first, it can be observed that the mean-field BdG calculation ( M = 1) missesseveral excited states in the low-energy spectrum. This is due to the fact that the BdG theory by construction39 = 1 M = 3 M = 4 ( m x , n x ) Exact ε ,
0) 4.524922 ω ,
0) 1.000000 ω ,
2) 1.788854 ω n/a 2.000005 2.000000 (2 ,
0) 2.000000 ω ,
3) 2.683282 ω n/a 2.788858 2.788854 (1 ,
2) 2.788854 ω n/a 3.000009 3.000000 (3 ,
0) 3.000000 ω ,
4) 3.577709 ω n/a 3.580111 3.577712 (0 ,
4) 3.577709 ω n/a 3.685038 3.681347 (1 ,
3) 3.683282 ω n/a 3.790078 3.788483 (2 ,
2) 3.788854 ω n/a 4.002097 3.999997 (4 ,
0) 4.000000
Table 2: Benchmark of the LR-MCTDHB implementation to the repulsive 1D HIM with N = 10 bosons. Shown are the ground-state energy ε and the energies ω i of the first few excitations for different values of M . The trapping frequency is Ω = 1 . λ = − .
01, yielding a marginal ground-state depletion of f ≈ .
03% for M = 4. Whereas theBdG ( M = 1) spectrum misses many states, the many-body spectra for M = 3 and 4 contain all low-lying excitations. Even thetwo-fold degenerate state (0 , ω and ω ). The assignment of the corresponding singleGP state to ω is arbitrary, and one would need to analyse for instance the response density to identify whether it correlateswith the 7th or 8th many-body excitation. The overall accuracy of excitation energies increases with the number of orbitals.Underlined digits denote deviations from the exact values from Eqs. (160) and (161). All quantities are dimensionless. See textfor more details. does not account for excitations were several particles at a time are excited from the ground-state manifold,so-called multi-particle excitations. For further details see Appendix A.1.1. The second observation is that thenumerical accuracy is enhanced if more self-consistent orbitals are included into the ground-state description.Also the two-fold degenerate excitation (0 , , M = 3 and4 orbitals. In addition, the convergence towards the exact ground-state energy is observed.Table 3 shows excitation energies obtained for the isotropic 2D HIM with trap frequencies Ω x = Ω y = 1and N = 100 repulsive bosons with interaction parameter λ = − . − , × [ − ,
9) grid with 64 ×
64 grid points. At first, one sees that the BdG approach yields the energy ofthe lowest c.m. excited state, which is two-fold degenerate, to perfect accuracy. All other states obtained arepure relative excitations. Higher c.m. excited states, as well as combinations of those with relative excitations,are missing in the BdG spectrum.In contrast to the BdG approach, LR-MCTDHB(2) gives access to more excited states. Before discussingthose, it is worth examining the two underlying ground-state orbitals. The first one is a Gaussian in the trapcenter, whereas the second one resembles a p x -orbital. The degree of depletion is marginal with f ≈ . M = 1, the spectrum contains the second-order c.m. excitation in the x -direction, as well as three out of six states corresponding to ω , which are excited states combining a first-orderc.m. excitation in the x -direction and relative excitations. The accuracy of levels that were already accessiblewithin the BdG approach has increased, e.g., for the first state of ω or the states of ω . Nevertheless, somestates have lost accuracy. The most remarkable example is the second first-order c.m. excitation of ω , whichis a c.m. excitation in the y -direction. Naturally the question arises why this is the case. The reason forthis is that the second ground-state orbital, as discussed above, resembles a p x -orbital. This means that theground-state description for M = 2 orbitals shows a preferred direction, which in this case is the directionalong the x -axis. Apparently, the consequence of this is that excitations in the y -direction either lose inaccuracy or remain inaccessible. A possible reason for the loss in accuracy is that, compared to the BdGcase, a small amount of probability weight of the first orbital has been transferred to the second orbital that40 i ( γ ) M = 1 M = 2 M = 3 ( m x , m y , n x + n y ) Exact ε (1) 89.554453 89.551373 89.548293 (0 , ,
0) 89.548292 ω (2) 1.000000 1.000000 1.000000 (1 , ,
0) 1.0000001.000000 0.999792 1.000000 (0 , ,
0) 1.000000 ω (3) 1.791089 1.788859 1.788859 (0 , ,
2) 1.7888541.791089 1.791089 1.788861 (0 , ,
2) 1.7888541.791089 1.791143 1.788859 (0 , ,
2) 1.788854 ω (3) n/a 2.000438 2.000439 (2 , ,
0) 2.000000n/a n/a 2.000657 (1 , ,
0) 2.000000n/a n/a 2.000439 (0 , ,
0) 2.000000 ω (4) 2.686634 2.683292 2.683226 (0 , ,
3) 2.6832822.686634 2.684452 2.683226 (0 , ,
3) 2.6832822.686634 2.685519 2.683312 (0 , ,
3) 2.6832822.686634 2.686634 2.683312 (0 , ,
3) 2.683282 ω (6) n/a 2.789191 2.787088 (1 , ,
2) 2.788854n/a 2.791314 2.787088 (1 , ,
2) 2.788854n/a 2.793711 2.789270 (1 , ,
2) 2.788854n/a n/a 2.789270 (0 , ,
2) 2.788854n/a n/a 2.791536 (0 , ,
2) 2.788854n/a n/a 2.791535 (0 , ,
2) 2.788854
Table 3: Benchmark of the LR-MCTDHB implementation to the isotropic HIM in 2D. Results are presented for N = 100 bosonsand different numbers of orbitals M . The trapping frequencies are Ω x = Ω y = 1 .
0, and the interaction strength is λ = − . f ≈ . M = 2 and f ≈ . M = 3. Shown are the energies of the ground state ε and the first 5 excited states, ω i = E i − ε , together with their multiplicities γ . In the BdG case, M = 1, only the first purec.m. excitation is found, together with all pure excitations of relative coordinates. For LR-MCTDHB(2), one obtains more c.m.excited states in the x -direction, together with combinations of excitations in the relative coordinates. For LR-MCTDHB(3), oneobtains the missing c.m. excited states in the y -direction and their combinations with excitations of the relative coordinates. Theoverall accuracy increases with M , with few exceptions in the spectrum for M = 2 which are discussed in the text. Underlineddigits denote deviations from the exact values from Eqs. (160) and (161). All quantities are dimensionless. See text for moredetails. effectively only covers the x -direction. It is stressed that rotating the second orbital by an arbitrary angle,e.g., π/ π/
2, leads to exactly the same spectrum, which means that it does not matter for the numericalresults which direction is preferred.Including M = 3 orbitals significantly enhances the obtained spectrum. One observes that both problemsof the spectrum for M = 2, i.e., the loss in accuracy for certain states as well as the absence of excitations in thelow-energy spectrum, are solved. The third ground-state orbital resembles a p y -orbital, i.e., it is perpendicularto the preferred direction of the M = 2 ground state. One can deduce from this observations a fundamentaldifference between the computation of the ground state of a given N -boson system and its excited states.Whereas increasing the number of orbitals always leads to a lower ground state energy, it can happen thatthe spectrum for ( M + 1) orbitals contains states that are less accurate than for only M orbitals. This is inparticular likely to happen in 2D and 3D because preferred directions may occur. As a result, in order toobtain accurate excitation spectra, one has to pay much more attention to the symmetry of the problem inhigher spatial dimensions.To summarize, it was demonstrated for the HIM that both the MCTDHB and LR-MCTDHB implementa-tions utilized in this work are well-suited to obtain accurate results for the spectrum and dynamics of trappedinteracting bosons. Furthermore, the many-body approaches clearly improve the GP and BdG mean-field41ethods. As mentioned above, further benchmarks are given in Appendix C. Since the advent of MCTDHB in 2007, the amount of scientific publications where it is utilized grew sub-stantially. As one of the first applications, the splitting of 1D repulsive, fully-coherent BECs by a potentialbarrier has been investigated [74, 75], showing that it can lead to fragmented condensates. Similar observa-tions were made by splitting a radially-symmetric BEC in 2D with a ring-shaped barrier [76]. Also the reverseis possible, meaning that two initially independent and thus fragmented BECs with no overlap in space cancollide, interfere and build up coherence [77]. Also for attractive BECs, the splitting due to a Gaussian-shapedbarrier results in the formation of a superposition of two distinct bosonic clouds, a so-called caton, which isnot describable at the GP level [78]. Moreover, it has been shown for the first time that coherent attractiveBECs in 1D can fragment even without a potential barrier, namely when the energy exceeds a threshold value[79]. Fragmentation persists and may even intensify when more orbitals are included in the calculation [80].In the attractive case full convergence is difficult to achieve as it requires more orbitals, see also [81]. Whenan attractive BEC propagates towards a potential barrier and scatters from it, the system also evolves to behighly fragmented although it was initially condensed [82, 83]. The degree of fragmentation depends on thedegree of transmission through the barrier. Even bright-soliton trains, i.e., stable multi-hump matter wavesof attractive bosons, whose dynamics were believed to be fully describable at the mean-field level, are shownto lose their initial coherence [84].The splitting of BECs has been further investigated with respect to the generation of number- or phase-squeezed states applying optimal control theory [85, 86, 87, 89, 88]. Optimal control has also been utilizedin 1D bosonic Josephson junctions (BJJ), on the one hand to steer the system parameters such that anenhancement of the Shapiro effect was observed [90], and on the other hand to drive a BEC from an initialstate to a target state at the quantum speed limit [91]. Moreover, it has been demonstrated that Mach-Zehnder interferometry with ultracold bosons is relatively robust against the nonlinear interaction betweenthe particles [92].A remarkable correspondence between the onset of wave chaos at the GP level and the onset of fragmenta-tion at the many-body level has been found for BECs that scatter from either shallow periodic or disorderedpotential landscapes [93, 94], but the result holds also for other potentials. Thus, the development of wavechaos, i.e., the exponential separation of initially close states, can be seen as an indication that a many-bodytreatment is necessary. The distance between states is measured utilizing the L norm between two Hilbertspace vectors.Of potentially high relevance is the demonstration how the positions of individual particles can be con-structed from the many-body wave function by simulating single shots [95]. The latter are connected toexperiments on trapped BECs where typically an image of the density represents a histogram of a single shot.Over the last years, MCTDHB has been applied to a growing number of systems, e.g., to study thebreathing dynamics in 1D harmonic traps due to a quantum quench [96] or the many-body effects in solitonicexcitations of ultracold BECs [97, 98]. Other examples are the correlated dynamics of a single atom coupledto a trapped BEC by collisions [99] or the examination of the structure of mesoscopic molecular ions [100],which contributes to the currently very active research in this field.Recently, MCTDHB has also been applied to dipolar BECs in a double well [101] and in 1D lattices[102, 103]. In addition, the coupling of one- and two-component BECs in a cavity to a radiation field has beenof interest [104, 105, 106]. One of the main results of the latter research was the finding of a fragmented super-radiant phase when the pump power of the laser exceeds a critical threshold. This fragmented superradiancecannot be explained by the Dicke model which assumes the BEC to be a simple two-level system.In the remaining part of this section, the focus is laid on several applications where the tunneling dynamicsin traps (Section 6.1) and the dynamical fragmentation of initially coherent BECs (Section 6.2) is investigated,both in 1D and 2D systems. Afterwards, all applications of LR-MCTDHB that exist up to now are presented42n detail (Section 6.3). Special emphasis for the 2D systems under consideration is laid on the impact ofangular momentum on the appearance of many-body effects in the dynamics and excitation spectra. Subsequently, the tunneling dynamics of BECs held in double-well traps in both one and two spatial dimensionsare discussed in this section. The main purpose is to apply MCTDHB to solve the many-boson Schr¨odingerequation and to discover novel phenomena that are not captured by the standard methods.
The quantum dynamics of a tunneling BEC in a 1D BJJ have constituted a large research field in thepast decades. Before the first experimental realizations [107, 108], theoretical predictions concerning Joseph-son oscillations and self-trapping, i.e., the suppression of tunneling between the wells, were obtained by a(multi-mode) GP description [109, 110, 111, 112, 113, 114], a path-integral formulation [115], a two-modeBose-Hubbard (BH) approach [116] or a Fock-space WKB method [117]. A review of various theoretical andexperimental results is given in [118]. A more recent study dealt with the tunneling dynamics in a controlledBJJ, where the spin state of an ionic impurity controls the tunneling between the wells [119]. The latterextends the findings on the dynamics of trapped ultracold BECs due to the presence of an ion [120, 121]. Fur-thermore, the orbital Josephson effect due to a time-dependent driving potential was investigated numericallyby employing MCTDHB, not only in a double-well system [122].The main objective of this section is however to study the out-of-equilibrium tunneling dynamics of asingle-component BEC in a BJJ from a many-body perspective, and the method of choice is MCTDHB. Asdescribed below, the latter approach allows for the observation of fundamentally new physics dealing with, e.g.,reduced self-trapping and dynamical fragmentation of initially fully-coherent bosonic clouds. A comparisonto other popular theoretical tools like the GP approach and the BH model is made. The findings presentedbelow were recently published in Refs. [123, 124, 125, 126], where further details can be found. Recently,results of a comparable study where both the BH model and MCTDHB were used to analyze the tunneling offew-boson systems in asymmetric double-wells, investigating the interaction blockade that isolates the motionof a single particle in the vicinity of others, became available [127].The trapping confinement considered here reflects two harmonic wells of the form V ± ( x ) = 12 ( x ± (172)which are connected by a cubic spline in the interval | x | ≤ . λ ˆ W ( x, x (cid:48) ) = λ δ ( x − x (cid:48) ). The repulsionstrength is given by the mean-field parameter Λ = λ ( N − P L ( t ) = 1 N (cid:90) −∞ ρ ( x, t ) dx (173)where ρ ( x, t ) describes the time-dependent one-body density, see Eq. (26), and N refers to the number ofparticles in the trap.Since the MCTDHB results are compared to the BH model below, it is instructive to define the left- andright-well orbitals, φ L and φ R , with the ground and the first excited state of the trap, denoted by φ g and φ u due to their gerade and ungerade symmetry. These two states are energetically well below the barrier. Theorbitals φ L and φ R are defined via φ L,R ( x ) = φ g ( x ) ± φ u ( x ) √ . (174)43he BH parameters for on-site interaction U and hopping J are derived from φ L and φ R by U = λ (cid:90) | φ L ( x ) | dx, J = − (cid:90) φ ∗ L ( x )ˆ h ( x ) φ R ( x ) dx (175)and utilized to define the BH interaction parameter Λ BH = U N/ (2 J ). The two-site BH Hamiltonian thenreads ˆ H BH = − J (ˆ b † L ˆ b R + ˆ b † R ˆ b L ) + U b † L ˆ b † L ˆ b L ˆ b L + ˆ b † R ˆ b † R ˆ b R ˆ b R ) (176)where the operators ˆ b ( † ) L,R annihilate (create) a boson in the left and right well, respectively. Within thetwo-mode BH approach, the eigenvalues and eigenfunctions of the one-body RDM ρ (1) = (cid:18) ρ LL ρ LR ρ RL ρ RR (cid:19) (177)with ρ LL = (cid:104) N, | e + i ˆ H BH t ˆ b † L ˆ b L e − i ˆ H BH t | N, (cid:105) (with (cid:126) = 1) and the other elements defined analogously, repre-sent the natural occupations and natural orbitals of the BH method. Therefore, the population of the secondnatural orbital can be used as a measure for fragmentation.In the following, the tunneling dynamics when the BEC is initially trapped in the left well, i.e., P L (0) = 1,will be explored for different repulsion strengths. Most importantly, for a two-mode GP description of thesystem, self-trapping of the bosons is predicted for Λ BH > Λ c = 2 when the magnitude of the initial imbalancebetween the occupations of the two wells, denoted by z = N L − N R N , is unity [109, 110, 111].Fig. 5 shows a comparison between the short-time dynamics of P L ( t ) for BECs with N = 20 (100)bosons obtained from MCTDHB (solid blue), BH (dashed magenta) and GP (dotted black) for two differentinteraction strengths where one is below and one is above the critical value Λ c . One observes that for the caseof 20 bosons with weak interaction parameter Λ BH < Λ c [panel (a)], the GP prediction is totally different thanthe BH and numerically exact MCTDHB predictions. According to GP, Rabi oscillations of the entire cloudbetween the wells occur, whereas the BH model and MCTDHB predict a density collapse after approximatelythree Rabi cycles, resulting in roughly 50% of the bosons in each well. The time for one Rabi cycle is closeto the theoretical prediction of t Rabi = π/J . Although the BH curve is clearly superior to the GP curve, italso deviates from the numerically exact result after half a Rabi cycle. In the case of 100 bosons [panel (b)],similar observations can be made. However, the density collapse occurs after a longer time.As mentioned above, for stronger repulsion with Λ BH > Λ c , tunneling should be entirely suppressed dueto self-trapping according to a two-mode GP description. It can be deduced from the exact dynamics of 20bosons [panel (c)] that indeed less bosons tunnel between the wells, meaning that the oscillation amplitudeis clearly smaller than for the weaker repulsion. As predicted by MCTDHB, at most 50 bosons tunnel intothe right well after half a Rabi cycle, and even less in the BH dynamics. Nevertheless, tunneling is notcompletely suppressed. This means that the effect of self-trapping, as predicted by the two-mode GP theory,is obviously reduced at the accurate many-body level. The deviations between the MCTDHB and BH resultsare more pronounced for stronger than for weaker repulsion. The dynamics of N = 100 bosons yields similarobservations [panel (d)]. In contrast to that, the GP theory predicts tunneling of the whole cloud.The reason why GP fails to describe the dynamics correctly is the development of fragmentation, alreadyon the short time scales shown. For the case of weak repulsion in the upper panels of Fig. 5, the initiallycondensed BEC with f = 10 − (10 − ) for N = 20 (100) evolves to be fragmented by 33% (26%) after threeRabi cycles. Interestingly, for stronger repulsion, the degree of fragmentation of the again initially fully-coherent BEC is less after three Rabi cycles than for the case of weaker repulsion, yielding approximately28% and 18% for 20 and 100 bosons, respectively. The predicted GP dynamics differ both quantitatively andqualitatively from the exact results. Most important, for all system parameters used in Fig. 5, GP predictsthat essentially all bosons tunnel back and forth between the wells. This is clearly not the case at the many-body level, neither for BH nor for MCTDHB. It is stressed that although the chosen system parameters are44 igure 5: Short-time dynamics of the occupation probability of the left well, P L ( t ), for (a) N = 20, Λ = 0 . N = 100,Λ = 0 .
152 (Λ BH < Λ c ) (c) N = 20, Λ = 0 . N = 100, Λ = 0 .
245 (Λ BH > Λ c ). The GP [dotted black], BH [dashedmagenta] and numerically exact results from MCTDHB( M ) [solid blue] are compared. For the weak interaction in panels (a)and (b), the GP theory predicts Rabi oscillations of the full density. At the many-body level, a collapse of these oscillationsoccurs after approximately three cycles. For stronger interaction in panels (c) and (d), GP still predicts that all bosons tunnelbetween the wells, whereas at the many-body level tunneling is clearly suppressed. However, MCTDHB and the BH model yielddifferent dynamics in all cases. The deviations become more pronounced for stronger repulsion. The inset shows the numericalconvergence with respect to M for 2, 4, 6 and 8 orbitals [(a) and (c)] and 2 and 4 orbitals [(b) and (d)]. All quantities aredimensionless. See text for details. The figure is taken from Ref. [123]. within the expected regime of validity for both the GP and BH approaches, both theories fail to describe theshort-time dynamics of the BECs.To further illustrate the effect of reduced self-trapping as well as the degree of fragmentation, Fig. 6 showsthe evolution of P L ( t ) and of the occupation of the first few natural orbitals n (1) i ( t ), both for MCTDHB andthe BH model. The utilized repulsion parameter exceeds the critical value significantly with Λ BH = 47 . . N = 10 (100). From the top panel, it can be inferred that the BH model predicts complete self-trappingfor both values of N , meaning that all bosons remain in the left well throughout the time-evolution. Onthe contrary, the predictions of MCTDHB show that indeed bosons tunnel from left to right, and one cananticipate that P L ( t ) approaches the long-time average of 0 .
5. With respect to the natural occupations, oneobserves that the BEC remains fully-condensed in time at the BH level, while the MCTDHB results show thatfragmentation sets in already during the first Rabi cycle. Moreover, more than ten orbitals are necessary toaccurately describe the dynamics. The different results are associated with a quick loss of coherence between45 igure 6: Upper panel: Dynamics of P L ( t ) for very strong repulsion. The interaction parameter Λ BH = 47 . .
4) for N = 10[solid blue] ( N = 100, solid green) is clearly above the critical value of Λ c = 2. While the BH predicts full self-trapping of thebosons [dashed magenta], the exact results show tunneling between the wells. Lower panel: Evolution of the natural occupations n (1) i ( t ). At the BH level, the system remains condensed throughout the propagation time, while it fragments at the exact many-body level where essentially four orbitals become macroscopically occupied. All quantities are dimensionless. See text for details.The figure is adapted from Ref. [123]. the bosons that is only described by MCTDHB and not at the BH level. In this context, see Ref. [123] formore details.Another interesting aspect is the symmetry of the two-mode BH Hamiltonian in Eq. (176) under theunitary transformation ˆ R = { ˆ b L → ˆ b L , ˆ b R → − ˆ b R } , leading to ˆ R ˆ H BH ( U ) ˆ R = − ˆ H BH ( − U ). A consequence ofthis is the equivalence of the occupation of the wells for attractive and repulsive interactions, i.e., P L,R ( t ; U ) = P L,R ( t ; − U ) . (178)Moreover, it can be shown that the eigenvalues of the BH one-body RDM in Eq. (177) do not depend on thesign of the interaction. Thus, P L ( t ) and the occupations n (1) i ( t ) will be the same for attraction and repulsionof the same magnitude. Naturally, the full many-body Hamiltonian of the system does not possess such asymmetry, and one can therefore expect different dynamics for attraction and repulsion between the bosons.The left panel of Fig. 7 shows the short-time dynamics of p L ( t ) for a BEC with N = 20 bosons withinteraction parameter | U | /J = 0 . | λ | = 0 . λ . On the time scale shown, the BH model mainly underestimatestunneling for the repulsive case, whereas it overestimates tunneling for the attractive case. The right panelof Fig. 7 shows the evolution of the first few natural occupation numbers for the same system as in the leftpanel, both at the BH and full many-body level. As expected, the BH model predicts the same evolutionfor n (1)1 ( t ) and n (1)2 ( t ), irrespective of the sign of the interaction parameter. In contrast to that, the results46 igure 7: BH versus exact dynamics for both attractive and repulsive BECs with system parameters N = 20, | U | J = 0 . | λ | = 0 . P L ( t ) of the left well. While the BH predicts the same evolution forattraction and repulsion, differences can be observed at the exact many-body level. Inset: Time-evolution of the ratio of P L ( t )for attraction and repulsion. Right panel: The same for the first two natural occupation numbers n (1) i ( t ). Again, the BH modelpredicts equivalent dynamics for attraction and repulsion, whereas the exact results show clear differences. All quantities aredimensionless. See text for details. The figures are taken from Ref. [124]. obtained from MCTDHB show that the evolution of the different occupations are distinct for the attractiveand repulsive cases. Compared to the exact results, the BH model mainly overestimates the fragmentation ofthe attractive BEC, whereas it mainly underestimates it for the repulsive BEC.The analysis of the dynamics in a 1D BJJ is closed by discussing the compelling feature of universalfragmentation. The latter means that all BECs with identical interaction parameter Λ fragment to the samevalue after the density has collapsed. Fig. 8 shows exactly this behavior for Λ = 0 .
152 and BECs with N = 100, 1000, and 10000 particles where initially all bosons are kept in the left well, i.e., P L (0) = 1. Fromthe upper panel, one sees that the density oscillations between the wells collapse after a few Rabi cycles,leading to the equal distribution of 50% of the bosons in each well. In general, the less bosons are contained inthe BEC, the faster the density collapses. Coming along with the damping of the oscillations, the BEC startsto fragment, reaching its maximal value after the collapse. Again, the less bosons are in the system, the fasterthe system fragments. Most importantly, irrespective of N , the final degree of fragmentation is the same.The fact that the degree of fragmentation appears to be universal is unexpected, since usually fragmentationstrongly depends on the number of particles in the BEC. Moreover, the degree of fragmentation upon the47 igure 8: Universality of the degree of fragmentation. Upper panel: Evolution of P L ( t ) for different numbers of bosons ( N = 100in blue, N = 1000 in red and N = 10000 in orange) where the interaction parameter Λ = λ ( N −
1) = 0 .
152 is kept fixed.Intermediate particle numbers ( N = 200 , , density collapse, denoted by f col , depends highly on the initial imbalance z of bosons in the left and rightwells. The larger z , the stronger does the system fragment in time. This holds true even in the limit ofvery weak interactions, Λ (cid:28)
1, meaning that there is no weak-interaction regime where GP yields the correctdynamics. It is important to note that the universality is not an artifact of any approximation, since theresult is obtained by solving the full many-boson Schr¨odinger equation numerically exact. Interestingly, thetwo-mode BH approach can be used to derive an analytic prediction of f col yielding results close to the exactones obtained by MCTDHB. Additional details on the analytic expression can be found in Ref. [126].To summarize, it was demonstrated that the out-of-equilibrium tunneling dynamics of BECs in a 1D BJJshow significant many-body features that are not accounted for at the GP mean-field level. This includesthe appearance of a density collapse in the oscillations between the wells as well as dynamical fragmentationwhich the GP theory cannot capture by construction. The degree of fragmentation is universal for fixedΛ = λ ( N −
1) for a wide range of initial states. Applying the commonly used BH model, it was shown that itis superior in certain aspects in comparison to the GP mean-field approach. However, it does also not predictthe numerically exact tunneling dynamics as obtained by MCTDHB, e.g., on the short-time dynamics wherethe interaction parameter is in the vicinity of its critical value or for very strong repulsion where it predictscomplete self-trapping of the bosonic cloud. However, the BH model can be utilized to derive an analytic48xpression for the universal degree of fragmentation that is in good agreement with the exact results.
After having shown the many-body nature of the tunneling dynamics of BECs in a 1D BJJ, the question ofhow a repulsive BEC behaves in a radially symmetric 2D system with a double-well structure is addressed.In particular, as for the previous section, the out-of-equilibrium tunneling dynamics are of interest, and it isstudied whether many-body effects can be observed. The results presented in this section were published inRef. [128], and additional details can be found therein.To study tunneling phenomena in 2D geometries has been of high interest in recent years, especiallyfor trapped vortices, i.e., macroscopic quantum states with definite angular momentum, characterized by adensity node at its core and phase discontinuities around it. Research concerning tunneling of vortices hasbeen carried out for example in 2D superfluids [129]. Others dealt with vortices tunneling through a Gaussian-shaped barrier [130], in pinning potentials [131, 132, 133], or between two Gaussian wells [134], to mentiononly a few. However, there has been no investigation of the tunneling dynamics of vortices in a radially-symmetric double well, in particular not at the many-body level. The latter geometry is very interesting tostudy many-body dynamics since the trap symmetry preserves the angular momentum of the bosonic cloud.It has been utilized recently to study many-body vortices which are separated in space due to the centralbarrier [135]. The following results, recently published in Ref. [128], especially deal with the impact of angularmomentum on the tunneling dynamics.Accurate many-body dynamics are obtained by solving the full many-body Schr¨odinger equation utilizingMCTDHB and compared to the corresponding mean-field dynamics from the GP theory.For the dynamics, the radial double well in which a BEC with N = 100 spinless bosons is confined reads V ( r ) = (cid:40) B e − r − R B ) + C e − . r − R C ) if r = ( x + y ) / ≤ R C C if r > R C (179)where B is the height of the central ring-shaped barrier located at the radius R B , and C is the height of thecrater wall at radius R C . A schematic illustration of V is shown in Fig. 9(a). The region enclosed by thebarrier is denoted by IN, the external rim, i.e., the region between the barrier and the crater wall, is denotedby OUT. For the numerical calculations, C = 200 . B = 1 . R C = 9 . λ ˆ W ( r i , r j ) = λ πσ e −| r i − r j | / σ , (180)which circumvents the regularization problems of the contact interaction [137]. The repulsion strength willagain be expressed in terms of the mean-field parameter Λ = λ ( N −
1) in the following. The effect of therange of the interaction on the degree of fragmentation has been investigated in detail for repulsive BECs in1D and 2D single wells [138], essentially showing that it becomes increasingly dependent on the density themore long-ranged it is.The analysis starts with the ground-state energy of two distinct cases, namely for (i) having all particlesin the IN region with closed external rim and (ii) having all particles in the OUT region with closed trapcenter. The notation E INΛ and E OUTΛ , where the superscript denotes whether the BEC is held in the IN orOUT region and the subscript denotes the repulsion strength Λ, is utilized. The ground-state energies fordifferent radial positions of the barrier R B are calculated using the GP equation. The corresponding resultsare shown in Fig. 9(b) for total angular momentum L = 0 and Fig. 9(c) for L = N , i.e., where each bosonon average carries one unit of angular momentum. One observes that for all shown parameters (Λ = 0 , . E E E E R R L=0 (b)
E/N R B E E E E R R L=N (c) E / N R B Figure 9: (a) Schematic plot of the radial double well.The IN region denotes the trap center, the OUT regiondenotes the outer annulus. Both are separated by thecentral ring-shaped barrier. (b) Ground-state energy perparticle of a BEC with N = 100 bosons and total angularmomentum L = 0 for different positions R B of the radialbarrier. The results are calculated at the mean-field GPlevel ( M = 1). The quantities E IN0 (short-dashed red)and E OUT0 (solid green) denote the energies for the INand OUT regions for interaction strength Λ = 0, whereas E IN2 (dash-double-dotted blue) and E OUT2 (dash-dottedmagenta) denote the energies for Λ = 2. The crossingpoints, i.e., the radii at which the energies for IN andOUT are alike, are located at R = 3 .
271 (Λ = 0) and R = 3 .
412 (Λ = 2), see the vertical black lines. Largerrepulsion between the bosons shifts the crossing point tolarger radii. (c) Same analysis as in (b) but for totalangular momentum L = N , meaning that each particlecarries on average one unit of angular momentum. Thelocation of the crossing points are R = 4 .
089 (Λ = 0)and R = 4 .
093 (Λ = 0 . L shifts them to larger radii, which isqualitatively similar as for increasing Λ. However, onewould need to highly increase the interaction strengthin order to shift the crossing points as much as addingone unit of angular momentum per particle does. Allquantities are dimensionless. See text for details. Figs.(b) and (c) are adapted from Ref. [128]. L = 0 and Λ = 0 , . L = N ) the curves for E INΛ and E OUTΛ intersect at certain radii which arereferred to as the crossing points or simply the crossings in the following. These are defined as the positionswhere the ground-state energies of the BEC in the IN and OUT regions are identical. Putting the barrieronto the crossing point is the 2D analogue of the symmetric double well in 1D. However, the IN and OUTregions are only energetically equivalent, whereas their geometries are totally different. It is found that boththe interaction strength Λ and the angular momentum L push the crossing point to larger radii once they areincreased. Though, the effect of L is significantly stronger since already for Λ = 0, the transition from L = 0to L/N = 1 leads to a shift of the crossing point from R B = 3 .
271 to R B = 4 . M = 4 orbitals are f ≈ − for L = 0and Λ = 2 . f ≈ − for L = N and Λ = 0 .
2, respectively. This means that in both cases, the systems’ground states are essentially condensed.To study the out-of-equilibrium dynamics for different values of Λ and L , the focus is laid on the situationwhere the barrier is located on the corresponding crossing points. First, the case of zero angular momentum,i.e., L = 0, is analyzed. As an initial state, the ground state of trapping the BEC in the OUT region withclosed central part is chosen. To trigger the dynamics, the system is quenched by suddenly opening the trapcenter, such that the entire potential is given by Eq. (179). It was found numerically that there is no differencein starting from either the IN or OUT region when the barrier is located at the crossing point. Fig. 10(a)shows the dynamical evolution of the occupation probability of the external rim, P OUT ( t ) = 1 − P IN ( t ) = 1 N (cid:90) R B 23, nearly thewhole cloud tunnels between IN and OUT. Then, the occupation of the second natural orbital starts togrow, accompanied by a damping of the tunneling oscillations. After approximately 12 cycles, the densityoscillations have collapsed and the particles are almost equally distributed between IN and OUT. The BECis now highly (two-fold) fragmented, with occupations of 56.7% and 41.1% of the first two natural orbitals α and α , respectively. In contrast to that, the GP dynamics do not show any damping of the oscillations,not even for the long-time behavior (not shown). Fig. 10(c) shows the real and imaginary parts as well asthe density of the GP orbital (top panels) and of the first two natural orbitals in the many-body case (middleand lower panels) at an instant in time at which the density oscillations have already collapsed. Whereasthere are no bosons in the trap center for the GP case, there is significant population of the IN region at themany-body level.Similar observations can be made for L = N and Λ = 0 . 2, presented in Fig. 10(b). The many-bodydynamics show that the density collapses and that approximately 50 bosons are occupying the IN and OUTregions finally. The BEC is two-fold fragmented, the occupation of higher natural orbitals is negligible. TheGP dynamics (see inset) do not account for the damped oscillations. With regard to the density of the GPand natural orbitals in Fig. 10(d), it can be seen that α and α clearly show the characteristic nodal vortexstructure in the center, whereas at the GP level there are no bosons in the IN region at the chosen instant oftime t = 183 . τ with τ = 66 . L = N and Λ = 2 . 0, the systemquickly becomes four-fold fragmented, i.e., at least four orbitals are necessary to account for the many-bodydynamics. For L = 0, the dynamics are essentially governed by only two orbitals, as can be seen from Fig.10(c). Thus, the weakly-interacting regime is apparently smaller for the case of L > 0, meaning that vorticesare more sensitive to many-body effects. In Section 6.3.3, this will be discussed again when excitations inrotating BECs are investigated.To summarize, the many-body nature of the tunneling dynamics of BECs in a 2D radial double well wasdemonstrated . The GP mean-field theory does not account for the collapse of the density oscillations betweenIN and OUT, and therefore fails to give an accurate description of the dynamics. Furthermore, the initially51 a) (b)(c) (d) Figure 10: (a) Time evolution of the occupation probability of the external rim, P OUT ( t ), for N = 100 bosons with L = 0 andΛ = 2 . 0. For M = 4, initially almost all bosons tunnel. The oscillation amplitude is damped over time and effectively collapsesafter 12 Rabi cycles with approximately half of the bosons in each part (solid red). The system evolves from being almost fullycondensed to two-fold fragmented (see solid green and dark blue curves for n and n ). The occupations of the third and fourthorbital are identical and negligible (solid light blue). On the contrary, the GP theory predicts undamped oscillations (dottedgray). (b) Same as in (a) but for L = N and Λ = 0 . 2. Again, the system evolves from being fully coherent to two-fold fragmented,and the particles finally occupy both the IN and OUT parts by roughly 50%. Mind the different time scales ( τ > τ ). Inset:Detailed view between 150 τ ≤ t ≤ τ . The GP dynamics again do not show the density collapse. (c) Real, imaginary, andabsolute value of the first two natural orbitals α and α at t = 18 . τ where the density is already collapsed. The GP case isshown in the top panels. For α and α , the density occupies both the IN and OUT regions. At the same time, the GP orbitalonly covers the external part. (d) Same as in (c) but for the dynamics described in (b). At the many-body level, the nodalstructure of the vortex state is clearly visible in the trap center. All quantities are dimensionless. See text for details. All figuresare adapted from Ref. [128]. condensed BECs fragment. It was found that bosons carrying angular momentum are more sensitive to many-body effects than non-rotating ones because fragmentation sets in already for weaker interaction strengths.In agreement with these findings, a recent work demonstrates analytically that a many-body approach isnecessary for describing the conservation of angular momentum in a 2D radially-symmetric system [139].52 .2 Dynamical fragmentation The obtained findings on the many-body tunneling process of initially trapped bosons to open space in 1D aresummarized in this section. The tunneling process of a single particle was understood already at the beginningof the 20th century [140, 141, 142, 143] and is still of current interest, e.g., with respect to the exit time andmomentum of an electron in the ionization process [144, 145]. Tunneling of a many-boson system, which hasalso been studied [146, 147, 148, 149], is naturally more complicated due to the interactions and correlationsbetween the constituent particles. The following results, published in Refs. [150, 151, 152, 153, 154], explainthe tunneling mechanism to open space from a many-body point-of-view utilizing MCTDHB. The main pointsthat are addressed deal with the questions whether the particles tunnel one by one or if several particles tunnelsimultaneously, and whether the entire process is indeed of many-body nature.To study the dynamics of the tunneling to open space, a system of N bosons confined in a harmonic trapof the form V ( x, t = 0) = x is considered. For t > 0, the trap is quenched such that it becomes open toone side. The region where the bosons are trapped initially and the open space are separated by a potentialbarrier. The full potential reads V ( x, t > 0) = θ ( x c − x ) 12 x + θ ( x − x c ) θ ( x c − x ) P ( x ) + θ ( x − x c ) T (182)where θ ( x ) denotes the Heaviside step function, P ( x ) denotes a third order polynomial that ensures a smoothconnection between the harmonic trap, the potential barrier and the threshold potential T whose role willbe discussed below. The coordinates x c and x c refer to connection points, where the former connectsthe harmonic trap and the barrier and the latter connects the barrier and the threshold. The bosons repeleach other via the contact interaction potential, and its strength is measured via the mean-field parameterΛ = λ ( N − t > N in the trap. Dependingon the threshold value T , certain states can become bound, meaning that no additional particles can escapefrom the trap via tunneling through the barrier. How the threshold can be utilized to control the tunneling-to-open-space dynamics will be discussed below.Important quantities for the analysis are the density ρ ( x, t ) [see Eq. (26)], as well as the density inmomentum space ρ ( k, t ), or simply the momentum distribution, which is obtained by performing a Fouriertransformation (FT) of the one-particle reduced density in real space. Furthermore, the degree of fragmenta-tion f , defined in Eq. (29), and especially its evolution in time, is of major interest. The above quantities arecomputed at the many-body level employing MCTDHB.The analysis starts by looking onto the shape of the momentum distribution in general. Fig. 12 shows ρ ( k, t ) for a system of N = 101 bosons with interaction parameter Λ = 0 . t < t < t < t . It consists of two parts, a Gaussian background and a peak structure. The former refers tothe particles that did not escape from the trap yet, since the initial real-space density is also a Gaussian. Onthe contrary, the peaks refer to the momenta of the already escaped particles. The heights of the peaks aretime-dependent, which can be seen from the growth of the dominant peak. The position of the peaks can beunderstood from a static model introduced below.The minimal total energy E TOT of the system is comprised of the energy of the bosons that remain inthe IN region, denoted by E HO ( λ , N IN ) indicating its dependence on the two-body repulsion strength λ andthe number of particles in the trap N IN , as well as on the minimal energy of the N OUT bosons that tunneledthrough the barrier into open space. Their minimal energy is the height of the threshold T . One thus obtains E TOT = E HO ( λ , N IN ) + N OUT T. (183)53 igure 11: Schematic plot of the potential V for t > ρ ( x, t = 0) is shown insolid blue. The IN and OUT regions are separated at the position of the maximum of the potential barrier (see solid red linesand labels). Colored horizontal lines denote the energies of states with different numbers of bosons N in the trap. The energydifferences between those states are denoted by the chemical potentials µ i . The value T denotes the height of the thresholdpotential in the OUT region. All quantities are dimensionless. See text for details. The figure is taken from Ref. [153].Figure 12: The momentum distribution ρ ( k, t ) for N = 101 bosons with repulsion strength Λ = 0 . t < t < t < t . The Gaussian background refers to the bosons remaining in the trap, whereas the peaks refer to the momentaof the emitted bosons. The height of the dominant peak is clearly growing in time. All quantities are dimensionless. See text fordetails. The figure is adapted from Ref. [152]. Assuming that the interaction among the emitted bosons is negligible, the kinetic energy of the i -th escapedboson is given by E kin ( T, µ i ) = µ i − T (184)where the { µ i } are chemical potentials that refer to the energy to bring an additional boson into the IN regionwith already N − i bosons. The corresponding momenta read k Ti = (cid:112) µ i − T ) (185)54here m = (cid:126) = 1 is assumed. It is stressed that this model describes the tunneling of a single boson. Figure 13: The tunneling to open space of N = 2 bosons. Left panel: Static consideration of the energy dependence of the possiblefinal states | N IN , N OUT (cid:105) on the threshold T . For T (cid:46) . 5, both particles are expected to tunnel, whereas for 0 . (cid:46) T (cid:46) . T , the bosons should remain trapped. Right panel: The obtained peaksin the momentum distribution ρ ( k, t = 600 , T ) for different T . For T ≤ . 4, two peaks appear at k T , , for T = 0 . . 6, onlyone peak appears. This is in agreement with the expected number of ejected bosons from the trap. For larger T , the peaks areshifted to smaller momenta since the bosons have to overcome a larger threshold. All quantities are dimensionless. See text fordetails. The figures are adapted from Ref. [153]. These model assumptions are tested against a system with N = 2 bosons and two-body repulsion strength λ = 1. The left panel of Fig. 13 shows the dependence of E TOT on the threshold value T for the threepossible final states | N IN , N OUT (cid:105) , given by | , (cid:105) , | , (cid:105) and | , (cid:105) . The utilized notation | N IN , N OUT (cid:105) shouldnot be confused with a Fock state and is just used to count the particles in the two regions IN and OUT.For T (cid:46) . 5, the lowest-in-energy state is | , (cid:105) , i.e., the final state where both bosons have tunneled toopen space. For 0 . (cid:46) T (cid:46) . 8, the state | , (cid:105) is the lowest, meaning that only one boson is expected totunnel in the long-time dynamics. For even higher T , tunneling is suppressed completely because the state | , (cid:105) corresponds to the lowest energy. This is in full agreement with the results shown in the right panelof Fig. 13. There, the peaks in the momentum distribution ρ ( k, t = 600), corresponding to the momenta k T and k T of the emitted particles, are presented for different values of the threshold T . The k T -peaks areplotted upside-down for better visibility. The first observation is that increasing T shifts the peaks to smallermomenta, which is due to the higher energy that is necessary to overcome the threshold. Secondly, two peaksappear up to T = 0 . 4, whereas for T = 0 . . T (cid:46) . . (cid:46) T (cid:46) . T (cid:46) . N , the chemicalpotential corresponding to a two-boson ejection would be µ b ≈ µ , which, for the case of T = 0, yields k b, ≈ √ k . This is simply not observed in the Fourier spectrum of ρ ( k, t ). Thus, it can be deduced that55 igure 14: Evolution of the natural occupations of the first few natural orbitals. Left panel: The case of zero threshold, T = 0,utilizing M = 4 orbitals. The interaction strength is λ = 0 . 3. For all particle numbers shown, the system evolves from beingcondensed to being multi-fold fragmented. Right panel: The occupation of the first two natural orbitals for N = 2 particlesand λ = 1, but for different heights of the threshold T . In general, the higher the threshold the later the development offragmentation sets in. Moreover, the system with larger T fragment less than the ones with smaller T . All quantities aredimensionless. See text for details. The figures are taken from Refs. [152, 153]. the bosons tunnel to open space one at a time. Nevertheless, as mentioned above, the individual processes ofthe ejection of single bosons are not independent, meaning that they quantum-mechanically interfere.Finally, it is discussed whether the overall process of tunneling to open space is indeed a many-bodyprocess or whether it can also be described at the mean-field level. To this end, the time-evolution of theoccupations of the first few natural orbitals for systems with different numbers of bosons and different heightsof the threshold is studied in Fig 14. In the weakly-interacting regime considered here, all initial states of thebosonic clouds are essentially condensed. For the case of zero threshold (left panel), one observes that for allnumbers of bosons considered the system becomes multi-fold fragmented. Especially for N = 4 bosons, atleast four orbitals become macroscopically occupied. With respect to the case of N = 2 particles and non-zerothresholds, the behavior is similar, meaning that the initially coherent BEC becomes highly fragmented. Ingeneral, the onset of fragmentation happens at later times when T is larger. Furthermore, the degree offragmentation is largest in the absence of any threshold.To summarize, it has been found that the tunneling-to-open-space dynamics of repulsive bosons froma harmonic trap are of many-body nature. With time, the system evolves from being coherent to becomehighly fragmented, and thus a mean-field description is inapplicable even in the weakly-interacting regime.The bosons tunnel one by one, meaning that no multi-boson tunneling occurs where more than a singleboson at the same time is emitted jointly from the trap. However, the individual tunneling processes are notindependent, they quantum-mechanically interfere. The threshold in the open space can be utilized to controlthe number of emitted particles. In particular, it can be used to suppress tunneling completely. A detailedanalysis on the coherence of the system, showing that the emitted bosons lose coherence with the bosons inthe trap and among themselves, can be found in Refs. [150, 151, 152, 153, 154], as well as additional details. A recently found vortex type in a rotating and repulsive 2D single-component BEC, termed phantom vortex,is presented in this section. Different to the common understanding of a quantum vortex that, among otherfeatures, has a density node (or core) and a phase discontinuity surrounding it, a phantom vortex cannotbe detected in the condensate density. Although they were not found in single-component condensates sofar, coreless vortices have been predicted earlier for multi-component condensates, e.g., spinor condensateswhere the core of one species, i.e., bosons with a certain spin, is filled by another species [155, 156]. Since56he experimental realization of gaseous BECs [4, 6, 5, 157, 158], the nucleation of vortices in rotating BECswas of interest and it was found that there is a critical rotation frequency below which no quantized vorticesappear [161, 159, 162, 163, 164, 160]. In contrast to that, phantom vortices appear also below this criticalrotation frequency. Instead of a node in the density, they rather manifest themselves in nodes of the densitiesof underlying natural orbitals. Thus, phantom vortices are a many-body object that cannot be described atthe mean-field level. Moreover, their occurrence is accompanied by the onset of fragmentation. The followingresults were published in Ref. [165]. Earlier results on phantom vortices in a stirred BEC where, e.g., thefragmentation entropy is analyzed, can be found in Ref. [166].The system under consideration is a BEC of N = 100 spinless bosons. They interact via the sameGaussian repulsion as in Eq. (180) with interaction parameter Λ = λ ( N − 1) = 17 . 1. They are confined in atime-dependent trap that is varied stepwise. It reads V ( r , t ) = 12 (cid:0) x ( t ) + y ( t ) (cid:1) + 12 η ( t ) (cid:0) x ( t ) − y ( t ) (cid:1) (186)where η is the parameter of the anisotropy given by the second term in Eq. (186). The x − and y − coordinatesare varied in time according to (cid:18) x ( t ) y ( t ) (cid:19) = (cid:18) cos( ωt ) sin( ωt ) − sin( ωt ) cos( ωt ) (cid:19) (187)where the frequency ω is set to 0 . 78. Initially, the trapping potential is isotropic and harmonic, meaning that η ( t = 0) = 0, and the corresponding ground state is computed both at the mean-field level using the GPequation as well as at the many-body level using imaginary time-propagation of the MCTDHB equations.Afterwards, a slight anisotropy is added by increasing η from zero to its maximal value of 0 . t = 80.The resulting elliptic trap is kept between 80 < t ≤ η is ramped down to zero again between300 < t ≤ t = 500, the isotropic harmonic trap is kept constant. The entireprocedure is shown in the upper panel of Fig. 15, together with representative densities of the BEC for theindividual intervals. These densities already suggest that at no point in time there is a vortex contained init, as discussed in more detail below.Of particular interest for the analysis will be, apart from the density ρ ( r ), the natural occupation numbers,denoted here by ρ (NO) i but with the same meaning as the n i in Eq. (28), and the densities | φ i | of the naturalorbitals (mind the different notation as compared to Section 3 where the φ i denote the working orbitals ofMCTDHB). Furthermore, the time-evolution of the expectation values of the orbital angular momenta, givenby ( L z ) ii = (cid:104) φ i | ˆ L z | φ i (cid:105) , and of the total angular momentum per particle, L z /N = N (cid:80) Mi,j =1 ρ ij ( L z ) ij , areanalyzed. The time-argument is suppressed in the above mentioned quantities. Computations are carried outfor orbital numbers M = 4 (many-body) and M = 1 (GP). Additionally, signatures of phantom vortices canalso be seen in the phase of the first-order correlation function g (1) ( r | r (cid:48) ; t ) [see Eq. (D.6) in Appendix D]defined by S g ( r | r (cid:48) ; t ) = arg[ g (1) ( r | r (cid:48) ; t )] (188)as well as in the phases of the natural orbitals given by S i ( r , t ) = arg[ φ i ( r ; t )] . (189)The discussion of the results is split into two parts. At first, the appearance of phantom vortices is analyzedat an instant in time when the trap is already isotropic again. Afterwards, the degree of fragmentation in thesystem is discussed.With respect to the first part, the total density and the densities of the natural orbitals at t = 450, i.e.,when the anisotropy is ramped down to zero again, are examined in Fig. 16. One observes that the totaldensity [panel (a)] does not exhibit any vortex. On the contrary, several phantom vortices appear in thedensities of the natural orbitals. The densities of the first and fourth natural orbitals contain a charge-1 and acharge-3 phantom vortex, respectively [panels (b) and (e)]. The density of the third orbital contains a pair of57 igure 15: Top panel: Time-evolution of the anisotropy parameter η ( t ) (red line) for the 4 distinct intervals discussed in the text,together with representative densities (blue). The dynamic protocol is divided into a period of ’ramp up’, ’maximal anisotropy’,’ramp down’ and ’isotropic’. Second panel: Time-evolution of the first 4 natural occupation numbers ρ (NO) i ( t ). Fragmentationsets in at the end of the ramp up of the anisotropy. Third panel: Time-evolution of the orbital angular momenta ( L z ) ii ( t )and the total angular momentum per particle, L z ( t ) /N . Fluctuations are maximal when fragmentation sets in. Bottom panel:Comparison between the time-evolution of the total energy E relative to the initial energy E and L z ( t ) /N . The GP theoryunderestimates both quantities, especially the angular momentum is approximately three times less than at the many-body level.All quantities are dimensionless. See text for details. The figure is taken from Ref. [165]. charge-1 phantom vortices [panel (d)]. The mechanism responsible for the nucleation of the phantom vorticesin φ and φ is different to the one responsible for the pair in φ . The former mechanism is called nodemutation, and is explained for the charge-1 phantom vortex in φ . This phantom vortex originates from thenode of the initial second orbital φ , which resembles a 2 p x orbital that has an angular node. When angularmomentum is added to the system due to the anisotropy, the lobes of the 2 p x orbital spread out and tendto close this angular node, leaving an elliptic density node in the center. After the trap is isotropic again, aphantom vortex is left at the center, which now appears in the density of φ since the first and second orbitalhave switched labels at t ≈ 220 because of their changing occupation numbers [see second panel from thetop in Fig. 15]. Thus, the remaining phantom vortex was nucleated from an initial node. A similar processhappens for the charge-3 phantom vortex in the density of φ . It originates from φ ( t = 0) which resembles a2 p y orbital. There, the angular node mutated first to three single phantom vortices that finally merged intoa larger charge-3 phantom vortex. The fact that it finally appears in the fourth natural orbital is due to thelabel switching of φ and φ during the ramp-up of the anisotropy.The mechanism responsible for the phantom-vortex pair in φ [panel (d)] is called slow orbital-orbitalvorticity transfer because of the following observation. At t ≈ φ . When time progressed, the two phantom vortices started to move towardsthe density edge and disappeared finally at t ≈ φ , which afterwards continued to move towards the center. At t ≈ igure 16: Top panels: Total density ρ ( r ) (a) and densities of the natural orbitals (b)-(e) at t = 450, i.e., when the trap is alreadyisotropic again. Although no vortex is observed in the total density, phantom vortices in the orbitals φ , φ and φ appear.Bottom panels: Phases of the first-order correlation function (f) and of the natural orbitals (g)-(j). Phase discontinuities areobserved for all phantom vortices. All quantities are dimensionless. See text for details. The figure is taken from Ref. [165].Figure 17: Time-evolution of the fragmentation for different numbers of particles N . The same dynamical protocol of thetrapping potential as in Fig. 15 is used. For all considered boson numbers, the occupation ρ (NO)1 decreases clearly. In general,fragmentation sets in later for higher particle numbers. Results are computed using M = 2 orbitals. The interaction parameteris kept fixed at Λ = 17 . 1. All quantities are dimensionless. See text for details. The figure is taken from the supplementarymaterial of Ref. [165]. between this pair is the same as before in the density of φ . Thus, vorticity has been transferred betweenthe two natural orbitals over a period much longer than the rotation period of τ = 2 π/ω ≈ 8. This processbetween the second and third orbital can also be seen from the correlated orbital angular momenta ( L z ) and ( L z ) in Fig. 15 for t ≥ S g ( r | 0; 450) in panel (f), one observesthe occurrence of several so-called ghost vortices at the density edge. The latter are distinct from phantomvortices because they do not enter the bulk density of the condensate and do also not contribute much to thetotal angular momentum.The dynamical fragmentation that occurs during the time-evolution, as well as its consequences for theangular momenta and the energy of the system [see again Fig. 15] is now analyzed. The initial state isessentially condensed ( ρ (NO)1 ≈ . < t ≤ ρ (NO)1 decreases significantlyto roughly 40%, and the second, third and fourth natural orbitals become macroscopically occupied. At thesame time, the orbital angular momenta vary the most and several of them reach their maximum values.Then, the occupations and angular momenta of φ and φ remain more or less constant. For the orbitals φ and φ , both quantities change slightly in a correlated manner, i.e., when ρ (NO)2 increases ρ (NO)3 decreases bythe same amount. This again reflects the slow orbital-orbital vorticity transfer between φ and φ discussedabove. With regard to the total angular momentum and energy, both quantities vary strongly during theperiods of maximal anisotropy and ramp down, and saturate only afterwards. The corresponding GP valuesare substantially different since they underestimate the system’s gain of energy and angular momentum. Inparticular for the latter quantity, it stays below 0 . 5, well below the threshold for the occurrence of vortices atthe mean-field level at unit angular momentum per particle. Interestingly, L z /N for M = 4 orbitals is greaterthan this threshold, but still no vortices are observed in the density ρ ( r ). This can be rationalized because onecan show analytically that a vortex in the density of a fragmented BEC, i.e., a fragmented vortex, requiresmore than unit angular momentum per particle [see the supplementary information of Ref. [165]]. The reasonis that a fragmented vortex originates from coincident phantom vortices, meaning that all natural orbitalsmust have a phantom vortex at the same position, and all of these phantom vortices need to have differentcharges such that the natural orbitals are orthogonal. As seen from panels (b)-(d) in Fig. 16, there are nocoincident phantom vortices, and thus the fragmented system does not show any vortices in the density.The onset of fragmentation was also observed for different numbers of bosons. Fig. 17 shows the evolutionof ρ (NO)1 for different N between 10 and 10 , utilizing M = 2 time-adaptive orbitals. The interaction parameterΛ = 17 . In this section, details on current applications of LR-MCTDHB are presented. The amount of scientificpublications where it is utilized is still rather small because the theory itself has been developed only afew years ago [34, 35]. In addition to that, the efficient MPI-parallelized numerical implementation of LR-MCTDHB was a very demanding task, and all the technical challenges going along with it, extensivelydescribed in Section 4, were solved very recently. It is stressed again that without such an implementation, avast majority of problems cannot be treated at all, especially in D > (a) (b) Figure 18: (a) Numerical convergence of the first three center-of-mass excited states (labeled by 1, 2 (cid:48) and 3 (cid:48) ) with respect to thenumber of orbitals M for N = 10 bosons in a harmonic trap with trap frequency ω ho = √ 2. The strength of the contact repulsionis λ N = 1. For M = 5 orbitals, the energies obtained from LR-MCTDHB coincide with the analytic predictions. The numericalimplementation used is the precursor of the implementation described and benchmarked in Sections 4 and 5. (b) Comparisonbetween the excitation spectra obtained from BdG (bottom) and LR-MCTDHB(2) (top) for N = 10 bosons in a shallow doublewell. The interaction strength is the same as in panel (a). Excitations labeled with primes refer to states solely obtained at themany-body level. The intensity of the response to a linear perturbation ( f + = f − = x ) is denoted by red triangles, whereas theintensity of the response to a quadratic perturbation ( f + = f − = x ) is denoted by green squares. The insets show the trappotential and the ground-state densities. One observes that the many-body low-energy spectrum contains a lot more states thanthe mean-field spectrum. All quantities are dimensionless. See the main text and Ref. [34], from which both figures are taken,for further details, especially for the discussion of the excitation type shown at the top of panel (b). The main purpose of Ref. [34] is to introduce LR-MCTDHB by giving a derivation of its main equations,presented in Section 3.2.2 of this work. In particular, the equations for the submatrices as described in Eq.(93) are developed for the case of contact interaction between the bosons. As an application, excitations forBECs in harmonic and shallow double-well traps are analyzed. Fig. 18(a) shows the numerical convergencefor the first three c.m. excited states (labeled by 1, 2 (cid:48) and 3 (cid:48) ) of N = 10 bosons in a harmonic trap with trapfrequency ω ho = √ M . The strength of the repulsion is given by λ N = 1.One can observe that the c.m. excitations are accurately obtained for M ≥ f + = f − = x, x ) are presented for N = 10 bosons in a shallow doublewell of the form V ( x ) = b/ πx/ 3) + ω x / b = 5. The most important observationby comparing the number of obtained states for BdG and LR-MCTDHB(2) is that the many-body spectrumcontains a lot more states. This in turn means that although the degree of depletion is only about 0 . M > As a first application of LR-MCTDHB, the many-body nature of excitations in a 1D harmonic trap as wellas in shallow and deep symmetric double wells is analyzed. The low-energy spectra obtained from the BdGequation at the mean-field level and from LR-MCTDHB at the many-body level are compared. Furthermore,it is investigated which excitations are involved in the dynamics due to certain quench scenarios. The followingresults where published in Ref. [167].Excited states are studied in a harmonic trapping potential with a Gaussian barrier given by V ( x ) = ax + b exp (cid:0) − cx (cid:1) (190)where a denotes the harmonic trap frequency and B denotes the barrier height. For simplicity, (cid:126) = m = 1is assumed. Furthermore, the trapping parameters are set to a = 1 / c = 1. The bosons interact viathe zero-ranged contact interaction potential. The interaction strength is again measured by the mean-fieldparameter Λ = λ ( N − 1) where λ is chosen to be positive to account for repulsion. Fig. 19 shows V ( x ) fordifferent barrier heights, together with the corresponding ground-state densities. For the purely harmonic trap( b = 0), the BEC is essentially condensed with n = 99 . b = 5) and deep( b = 10) double wells, the ground state is already fragmented by approximately 5% and 40%, respectively.One can further observe that for the deep well, the overlap between the densities in the left and right wells ismarginal, whereas there is still a clear overlap observable in the shallow case.Since the trap is reflection invariant, the excitations can be categorized in gerade ( g ) and ungerade ( u )symmetry. Fig. 20 shows the low-energy spectrum of N = 10 bosons for different barrier heights 0 ≤ b ≤ . 0. Excitation energies are given relative to the ground-state energy, i.e., ∆ E = E − E . Onereadily observes the formation of bands, which is even more pronounced for lower interaction strengths [seeFig. 2 in Ref. [167]]. The bands are labeled by integer numbers, where the zeroth band denotes the onethat is lowest in energy. Comparing the mean-field and many-body approaches, one can see that the BdGspectrum only contains 5 different excitations, whereas the spectrum computed with LR-MCTDHB shows amuch richer structure where a lot more excited states appear. In terms of numerical convergence, M = 2orbitals are sufficient to describe the spectrum up to ∆ E ∼ b ≥ 4. Forstates higher in energy, at least four orbitals are necessary.In the following, different shift and quench protocols are applied to trigger excitations from the low-energyspectrum. For gerade excitations, a quench of the harmonic confinement from a = 0 . → . x → x + x shift with x shift = 0 . 1. The gerade excitations are referred to as breathing modes, and the ungerade excitations arereferred to as dipole modes. The dynamics of the bosonic clouds for the two cases are shown in Fig. 21 on ashort time scale. In the left panel, dipole oscillations for a BEC in the deep double well are depicted, whichmeans that in each well, the bosons oscillate around the equilibrium position given by the minima of thewells. On the RHS, the breathing dynamics for N = 10 bosons in the purely harmonic trap is shown. Onecan observe that the density spreads to the left and right in a symmetric manner, and returns to the minimumof the well. The Fourier spectra of the time-evolution of certain quantities is used to extract the spectrumof involved excitations in the dynamics. This comprises ( a ) the expectation value (cid:104) x (cid:105) , ( b ) the variance perparticle N ∆ X , and ( c ) the density ρ at the position x = 1. The many-body variance is discussed in detail inAppendix D. In Ref. [167], the Fourier spectra for b = 0 , igure 19: Trap potentials (solid red) and ground-state densities (gray) for different barrier heights b . The ground states arecomputed with MCTDHB(2). The BEC consists of N = 10 bosons and the repulsion strength is Λ = 1 . 0. The occupationnumber n of the first natural orbital shows that increasing the barrier height leads to fragmentation. Vertical lines refer tothe lowest-in-energy excited states computed with LR-MCTDHB(2). The green potentials and the blue eigenvalues refer to theharmonic approximation of the right well. All quantities are dimensionless. See text for details. The figure is taken from Ref.[167]. Table 4 shows excitation and de-excitation energies where the obtained values from the BdG and LR-MCTDHB(2) theories (columns 2 and 3) are compared to the obtained excitation energies of the FT spectrafrom the shift and quench scenarios (columns 4-7). Moreover, several de-excitations, i.e., transitions betweena higher and a lower excited state, are presented (last 4 rows). Those refer to higher-order effects and arein principle not obtained directly by LR. However, since LR yields the exact excitation spectrum once it isapplied to an exact eigenstate of the full many-body Hamiltonian, the de-excitations can be extracted fromdifferences between the excitation energies of the LR spectrum. The excitation labels are shown in the firstcolumn. Many-body excitations, i.e., excitation that only appear in the many-body spectrum, are denoted byprimed labels.Up to ∆ E ∼ 4, the many-body approach yields ten excited states, whereas on the contrary the BdGspectrum only consists of five excitations. The energies of states that appear in both spectra differ slightlyfrom each other for 1 g , 1 u , 2 g and 2 u , but significantly for the lowest excitation 0 u . That demonstrates thatstrong deviations between the mean-field and many-body excitation energies can already occur for the lowestexcited states. Also the intensities of several states differ clearly from each other, especially for the loweststates 0 u , 1 g and 1 u . The intensities from the LR calculations are identical to the response weights describedin Eq. (123), with f + = f − = x for the ungerade and f + = f − = x for the gerade states.Columns 4 and 5 show the obtained peaks in the Fourier analysis of (cid:104) x (cid:105) and ρ ( x = 1) for the shiftscenario, both for the GP and many-body cases. The Fourier spectra are shown in the left panels of Fig. 22.63 igure 20: Low-energy excitation spectrum of N = 10 bosons with interaction parameter Λ = 1 . b .Many-body results computed with LR-MCTDHB(2) are denoted by solid lines (gerade excitations in green, ungerade excitationsin red). Results for M = 4 orbitals are given by dotted lines, and the corresponding BdG energies are given by open black circles.All states up to ∆ E ∼ Surprisingly, the shift also triggered a gerade excitation (2 g ) at the mean-field level, whereas at the many-bodylevel only ungerade states are observed. Moreover, no peak is observed for 0 u in the many-body dynamics,presumably because its intensity is too weak to be distinguished from the background noise. The obtainedpeak positions coincide with the predictions from BdG and LR-MCTDHB(2), respectively. With regard tode-excitations, the mean-field and many-body approaches yield a very small peak for the 2 u ↔ u transitionin the density spectrum.The Fourier analysis of the quench dynamics of N ∆ X and ρ ( x = 1) are shown in column 6 (GP) andcolumn 7 [MCTDHB(2)]. Only peaks corresponding to gerade excitations are observed. Both approachesyield peaks for 1 g and 2 g . At the many-body level, two additional peaks referring to 0 g (cid:48) and 1 g (cid:48) are obtained.The peak positions again accurately match the predictions from the static analysis. Most importantly, theintensities of these excitations in the variance spectrum are comparatively large, and therefore 0 g (cid:48) and 1 g (cid:48) arethe most promising many-body excited states to be detected in an experiment. Furthermore, no de-excitationpeaks are found in the GP spectra, while the variance spectrum shows small peaks for the 2 g ↔ g and1 g ↔ g (cid:48) transitions at the many-body level, coinciding with the values from LR-MCTDHB(2).In general, the sensitivity of the utilized operators to the detection of excited states can be seen as follows.The spectrum of the local operator ρ ( x = 1) is in particular useful to detect ungerade excitations in theshift scenario. The intensities of the peaks of many-body excitations are however weak. To detect geradeexcitation, nonlocal quantities, especially the variance, turn out to be more sensitive than the density. Inaddition, the variance spectrum yields intense peaks also for many-body excitations.To summarize, it has been demonstrated that from a methodological point of view, utilizing LR-MCTDHB64 igure 21: Many-body dynamics of the density of N = 10 bosons with interaction parameter Λ = 1 . M = 1 M = 2 Shift GP Shift MB Quench GP Quench MB0 u a )0 g (cid:48) b )1 g b ,0.93 c ) 2.092(0.33 b ,0.89 c )1 u a ,0.66 c ) 2.095 (0.93 a ,0.7 c )1 u (cid:48) a ,0.02 c )1 g (cid:48) b )2 g c ) 3.914 (0.07 c ) 3.934 (0.07 b ,0.1 c )2 u a ,0.23 c ) 3.943(0.03 a ,0.15 c )2 u (cid:48) a ,0.07 c )2 g (cid:48) g (cid:48) ↔ g g ↔ g b )2 u ↔ u c ) 1.849 (0.02 c )1 g ↔ g (cid:48) b ) Table 4: Table of excitation and de-excitation energies for N = 10 bosons with Λ = 1 . (cid:104) x (cid:105) (superscript a ), the variance per particle N ∆ X (superscript b ) and the real-space density ρ ( x = 1) (superscript c ), for both the shift andquench scenarios. GP denotes the mean-field and MB the many-body case. The intensities are given in brackets. All quantitiesare dimensionless. See text for more details. The table is adapted from Ref. [167]. in order to obtain the many-body excitation spectra of repulsive BECs in 1D harmonic and double-well trapgeometries is a promising and highly accurate technique. Its predictions coincide with the Fourier spectra ofthe time-evolution of different quantities like the local density or the position variance. From a physical pointof view, some of the obtained excited states corresponding to many-body states not included in standard65 igure 22: Fourier spectra of (cid:104) x (cid:105) , N ∆ X , and ρ ( x = 1) from the shift and quench dynamics in the deep double well. The GPspectra are denoted by solid black and the MCTDHB(2) spectra by dashed red lines. The intensities are normalized by the sumof the intensities of all constituent peaks, and excitations with less than 1% intensity are omitted. Especially the spectrum ofthe many-body variance consists of several peaks that are not observed in the corresponding GP spectrum. All quantities aredimensionless. See text for details. The figure is adapted from Ref. [167]. mean-field approaches yield intense responses especially in the variance spectra and are therefore likely to bemeasured in experiments. More details on the presented results can be obtained from Ref. [167]. In this section, many-body effects in the excitation spectra of weakly-interacting BECs in 1D lattices areinvestigated. The main objective is to demonstrate with numerically converged results that a marginaldepletion of the order of 1% is sufficient in order to observe clear differences between the lowest-in-energyexcitations computed at the mean-field and many-body levels. The results presented below were recentlypublished in Ref. [168].Studying dynamics and excitations of ultracold bosons in 1D optical lattices was of high interest in thepast, and the literature is extensive in this research field. Examples are the dynamics in the superfluid phase[169], nonlinear dynamics in periodic potentials [170] or in fluctuation-driven binary condensates [171]. Areview of earlier works can be found in Ref. [172]. Of particular relevance for the findings discussed in thissection are studies were the BdG theory does not agree with numerical and experimental observations. It hasbeen found that the oscillation frequency of a superfluid BEC in a large optical lattice with a superimposedharmonic confinement deviates clearly from the BdG prediction [173, 174]. In addition, the strong depletion66f a gaseous BEC in optical lattices cannot be fully explained by the BdG theory, not only in 1D [175].Optical lattices, in particular triple-well systems, have also been analyzed recently utilizing MCTDHB.This includes the out-of-equilibrium dynamics of rather small BECs consisting of a few particles only, inducedby quenching the interaction and the lattice depth [176, 177, 178] and periodic driving of the system [179, 180].Of particular interest was the analysis of cradle modes, dipole oscillations corresponding to the transport ofbosons over the barrier between adjacent wells. Furthermore, the many-body dynamics of BECs initiallyprepared in Mott-insulating states have been of interest [181]. The different phases of the BEC, i.e., thesuperfluid, the Mott-insulating, and the fermionized gas phase, can be detected from the many-body entropyand the Glauber correlation functions [182]. It was even shown for a BEC in a tilted triple well with hard-wallboundary conditions how the correlations between wells can be controlled by the well depth, the interactionstrength and the tilt [183]. Many-body effects in the ground state of a bosonic mixture in a more complex2D lattice where studied using the multi-layer version of MCTDHB [184].Many-body excitations in optical lattices were addressed as well, either by making a Fourier analysisof certain quantities [176] or by doing an exact diagonalization. There, the system’s wave function wasexpressed using a many-body ansatz based on a number state expansion including localized Wannier functions[179, 180, 177, 178]. The latter approach of exact diagonalization is possible for small systems with a fewbosons only.Below, the lowest-in-energy excitations of BECs in 1D lattice potentials with periodic boundary conditionsare investigated using LR-MCTDHB. This opens the possibility to study also larger systems with more wellsand particles. The utilized trap potential reads V ( x ) = V cos (cid:16) πl x (cid:17) (191)where V denotes the lattice depth and l denotes the distance between two neighboring lattice sites. Thedepth of the lattice is commonly expressed in terms of the recoil energy E R = (cid:126) k m with the lattice momentum k = πl . For simplicity, (cid:126) = m = l = 1 is assumed. The bosons repel each other via the contact interactionpotential and its strength is again expressed in terms of the mean-field parameter Λ = λ ( N − N = 10 bosons confined in a triple-well potential with periodic boundary conditions areconsidered. The low-energy excitation spectra for both a shallow ( V = 5 . ≈ . E R ) and a deep well( V = 50 . ≈ . E R ) are discussed. Results are presented in Figs. 23(a) and (b). In both wells, thespectra for the non-interacting BEC are structured in sequences of levels with different degeneracies. Thestates and levels can be associated with appropriate quantum numbers, envisioning that for Λ = 0 eigenstatesof a single particle in a periodic lattice are exact quasimomentum eigenstates. The eigenstates appear ina band structure that is more pronounced for deep lattices. For the shallow lattice in panel (a), the firstlevel consisting of two degenerate states denotes the excitations of putting one boson in either a state withquasimomentum (+1) or with quasimomentum ( − − 1) from the second single-particle band. To this end, the notation (cid:16) n (1)+1 , n (1) − ; n (2)+1 , n (2) − (cid:17) is introduced, where the first two entries in brackets denote the occupation of (+1) and( − 1) from the first single-particle band and the last two entries denote the occupation of (+1) and ( − 1) fromthe second single-particle band. As another example, level 2 contains the states (2 , 0; 0 , , 2; 0 , 0) and(1 , 1; 0 , > 0, this is no longer possible, as discussed below. In panel (b), only the firstsingle-particle band appears in the shown energy range because, as shown in Fig. 25, the first and secondband are strongly separated from each other.For non-zero repulsion in the shallow well in panel (a), the BECs are slightly depleted ( f = 1 . 1% forΛ = 4 . . . BdG(1)BdG(2) f Λ =4.0 =1.1% ω n n Λ =4.0 Λ =3.0 Λ =2.0 Λ =0.0 (a) Λ =2.0 =8.4% BdG(1) ω n n Λ =2.0 Λ =1.0 Λ =0.5 Λ =0.0 (b)(c) (d) Figure 23: Low-energy spectra of N = 10 bosons in a triple well computed with LR-MCTDHB(7). (a) Shallow case with latticedepth V = 1 . E R . For Λ = 0, the excitations form levels of degenerate states. Degeneracies are lifted for Λ > 0. TheBdG results only contain the two doublets denoted by BdG(1) and BdG(2). (b) Same as in (a) but in the deep lattice with V = 10 . E R . The BEC for Λ = 2 . . . − 1) ofthe first single-particle band are occupied by N (1) bosons. The difference grows with N (1) . Inset: Evolution of the relative erroras defined in the main text. (d) Same as in (c), but for the system in (b) with Λ = 2 . 0. All quantities are dimensionless. Seetext for details. The figures are adapted from Ref. [168]. non-interacting case vanishes as Λ increases, which makes it more difficult to identify the individual stateswith regard to the categorization in terms of quasimomentum states. The corresponding mean-field resultscontain only two doublets for the shallow and one doublet for the deep well, denoted by BdG(1) and BdG(2).All other states, representing multi-particle excitations where more than one boson at a time is excited fromthe ground state, are absent.Figs. 23(c) and (d) demonstrate in how far the BdG energies can be utilized to anticipate the positions ofthe missing multi-particle states. As mentioned above, this yields the exact energies for the missing states inthe non-interacting system. In particular, multi-particle excitations of putting in total N (1) bosons into thestates (+1) and ( − 1) of the first single-particle band are discussed. With regard to the BEC with Λ = 4 . ≤ N (1) ≤ 5, i.e., the accurate many-body values are lower. The discrepancy between mean-field and many-68 BdG(1)BdG(2) f Λ =1.0 =0.84% ω n n Λ =1.0 Λ =0.5 Λ =0.1 (a) (b) Figure 24: (a) Same as in Fig. 23a but in a lattice with 10 sites. The depletion grows faster with Λ than in the triple well,meaning that the system is even more sensitive for many-body effects. (b) Same as in Fig. 23c but for 10 sites and Λ = 1 . 0. Therelative error (inset) reaches similar values as for systems in the shallow triple well with comparable depletion. All quantities aredimensionless. See text for details. The figures are adapted from [168]. body results grows with N (1) . Especially the relative error, defined as E rel = (cid:12)(cid:12)(cid:12)(cid:12) ω BdG ( N (1) ) − ω MB ( N (1) ) ω MB ( N (1) ) (cid:12)(cid:12)(cid:12)(cid:12) where ω MB (cid:0) N (1) (cid:1) denotes the many-body energy with the largest distance to ω BdG (cid:0) N (1) (cid:1) = N (1) · BdG(1), growsquickly to more than 10%. It is emphasized again that the BEC is only slightly depleted by approximately 1%.For the fragmented BEC with Λ = 2 . N (1) = 1) the relative error is about 5%. Thismeans that in this case BdG does not only miss all multi-particle excited states in the low-energy spectrum,but it is also very inaccurate for the single-particle excitations.Fig. 24 shows a similar analysis for N = 10 bosons in a larger shallow lattice with 10 sites. In principle,the described many-body effects observed in the triple well also appear in this extended lattice. However,with respect to the degree of fragmentation, which is f = 0 . 84% for Λ = 1 . 0, one can deduce that these many-body effects set in even earlier than in the triple well because a comparable depletion is achieved already forweaker repulsion. Hence, an accurate many-body description is of even larger importance if the lattice size isincreased.The results shown in Fig. 25 demonstrate numerical convergence with respect to the number of self-consistent orbitals M utilized. In all cases, M = 7 orbitals clearly suffice for the low-energy excitations. Inparticular, for the middle panel describing a BEC in the deep triple well, one can see that firstly, M = 3orbitals already give very accurate values for the excitation energies, and, secondly, all multi-particle statesbuild up by the states (+1) and ( − 1) from the first single-particle band are obtained (65 in total). To remindthe reader, it is stressed again that the BEC is fragmented in this case, and nevertheless LR-MCTDHB iscapable to give converged results.So far, it was shown that the BdG energies for a BEC in a 1D lattice can be largely inaccurate, evenfor weakly-depleted condensates with f ≈ N → ∞ with fixed Λ, the BdGequation yields the exact energies of all single-particle excitations, and the energies of the multi-particlestates are given by multiples and sums of the BdG energies [185] (extending the findings of a previous workon homogeneous systems [186]). Whether this result is applicable also for BECs with contact interactionremains unclear. However, naturally the question arises how fast, if at all, the excitation spectrum of a BECwith arbitrary shape of the repulsion approaches the BdG energies upon increasing the number of bosons.69 10 sites, V =5.0, Λ =1.0 ω n n M=7M=9 0 5 10 15 20 25 10 20 30 40 50 60 Triple well, V =50.0, Λ =2.0 BdG(2) ω n n M=3M=5M=7 0 5 10 15 20 10 20 30 40 50 Triple well, V =5.0, Λ =4.0 ω n n M=5M=7M=9 Figure 25: Proof of numerical convergence for the systems discussed in Figs. 23 and 24. LR-MCTDHB(7) clearly suffices in allcases. For the BEC in the deep triple well [middle panel], all 65 single- and multi-particle excitations below the bottom of thesecond band [denoted by BdG(2)] are obtained. Already for M = 3, the results are highly accurate which is due to the largeseparation of the first two single-particle bands. All quantities are dimensionless. See text for details. The figures are adaptedfrom Ref. [168]. This question is discussed further in the subsequent section, Section 6.3.3, where excited states in rotatingBECs are analyzed.Finally, the zero-quasimomentum modes (ZQMs) are briefly discussed. The latter are defined in thiswork as the excitations where MOD( P, L ) = 0, with P denoting the total quasimomentum up to the second70 igure 26: (a) Energies of the zero-quasimomentum modes (ZQMs) for the system in Fig. 23(a) with Λ = 4 . 0, up to the top ofthe second single-particle band [denoted by BdG(3)]. All shown states are not obtained by the BdG mean-field approach. TheZQMs are categorized by the symmetry of their density responses, which is either gerade or ungerade on each lattice site. Mindthe different index n z , which enumerates only the ZQMs compared to the index n of the upper figures which enumerates allexcitations. (b) and (c): Real part of the density responses of the states (3 , 0; 0 , 0) and (0 , 3; 0 , 0) [panel (c)] as well as (4 , 1; 0 , , 4; 0 , 0) [panel (d)]. The vertical lines separate the lattice sites from each other. Notice the different scales. All quantitiesare dimensionless. See text for more details. The figures are taken from Ref. [168]. single-particle band given by P = n (1)+1 + n (2)+1 − n (1) − − n (2) − , (192)and L denoting the number of lattice sites. The spectrum of the ZQMs up to BdG(3), which marks the top ofthe second single-particle band, is shown in Fig. 26(a). The system parameters are the same as in Fig. 23(a)with Λ = 4 . 0. It is worth noting that all of the shown ZQMs are multi-particle excitations, meaning thatBdG does not account for a single ZQM in the low-energy spectrum although there exist more than 20. Byanalyzing the response densities of the ZQMs, it was found that they can be categorized by their symmetrywith respect to the individual wells, which is either gerade or ungerade. Specific examples are given in Figs.26(b) and (c). The ZQMs are most likely the easiest states to excite in an experiment where the BEC is subjectto either a quench of the potential or the interaction strength. For example, a simple quench scenario to excitethe gerade ZQMs would be a sudden change of the lattice depth. To this end, the ground state of a BEC with N = 10 bosons in the triple well with interaction parameter Λ = 4 . V = 4 . V = 5 . 0, and the subsequent time-evolution is investigated. Fig. 27shows the Fourier spectrum of the dynamical evolution of the position variance per particle, N ∆ X , obtainedfrom both MCTDHB(5) and GP. Shown is the energy range 5 . ≤ ω ≤ . 5. At the many-body level, 3distinct peaks are observed, whose positions correspond to the first, second, and fourth gerade excitations ofthe spectrum in Fig. 26(a). They can thus be labeled according to the occupation of (+1) and (-1) from thefirst and second single-particle band. The peaks corresponding to the two-particle excitations at ω ∼ . ω ∼ . 48 are more intense than the peak corresponding to the three-particle excitation at ω ∼ . Intensity ω M=5GP (1,1;0,0) (3,0;0,0) or (0,3;0,0)(1,0;0,1)or (0,1;1,0) Figure 27: Low-energy Fourier spectrum of the time evolution of the position variance per particle in x -direction, N ∆ X , dueto a quantum quench of the lattice depth from V = 4 . V = 5 . 0. The system parameters are the same as in Fig. 23(a)with Λ = 4 . 0. The many-body spectrum (dashed red), obtained from the evolution with M = 5 orbitals, shows two sharp peaksreferring to two-particle excitations and one less intense peak referring to a three-particle excitation. The first intense peak refersto the two-particle excitation (1 , 1; 0 , , 0; 0 , 1) [or (0 , 1; 1 , , 0; 0 , 0) [or (0 , 3; 0 , completely flat . The depletion of the BEC is f ∼ . 1% during the entire propagation time. All quantitiesare dimensionless. See text for more details. of f ∼ 1% is sufficient to observe these effects. In larger lattices, the sensitivity to many-body features evenincreases. It was further shown that a special type of many-body excited states (ZQMs) are easily accessibleby performing a simple quantum quench to the system. A promising quantity to detect these states is theposition variance, as already seen in the previous section. After the discussion of many-body excited states in 1D trapped BECs, an example in 2D, namely the case ofrotating BECs in an anharmonic, radially-symmetric potential crater, is examined. Of particular interest isthe weakly-interacting regime in which the ground-state depletion is marginal and one would naively expectmean-field approaches like the BdG theory to be valid and accurate. The results below were recently publishedin Ref. [187] and it is referred to this work for additional details.Rotating BECs were of high interest in the past. Research dealt with the occurrence of quantized vortices[161, 188, 162, 163, 164] or with the similarity to the fractional quantum Hall effect for charged particles inthe vicinity of a magnetic field [189, 191, 190, 192]. Moreover, low-lying excitations in such systems wereaddressed, for example in vortex lattices [26, 193, 194, 195]. Furthermore, investigations were made on the72ecay of the counter-rotating quadrupole mode [196], on Tkachenko modes [197, 198], on the twiston spectrum[199], or on excitations in anharmonic traps [200, 201]. In the latter works, the mean-field excitations werecalculated by utilizing the BdG theory. Computations employing a many-body description arerather rare.An examples is the analysis of the yrast spectra in a harmonic confinement obtained by exact diagonalizationusing the lowest Landau level approximation [203, 204, 202, 205].As mentioned above, accurate many-body results for the low-energy spectrum of rotating BECs in ananharmonic, radially-symmetric trapping potential are presented. Especially, the effect of the angular mo-mentum in the BECs’ ground states on the excitation energies is discussed. As for the examples in the earliersections, the BdG and many-body results are compared to each other. Figure 28: (a)-(e) Ground-state densities of rotating BECs with N = 10 bosons and interaction parameter Λ = 0 . l . Results are shown for MCTDHB(7), the corresponding GP densities look alike (not shown). For non-zero vorticity,the ground state is a single vortex whose core size grows with its charge l . (f) Low-energy excitation spectra atop the groundstates from the upper panels. Many-body results obtained by LR-MCTDHB(7) are denoted by red and blue squares, the BdGresults are denoted by black and blue lines. Once the ground state has non-zero vorticity, deviations between the mean-field andmany-body excitation energies are observed, and they grow with l . Moreover, the density of (many-body) excited states growsstrongly. All ground states are only slightly depleted, with a maximum of f ≈ . 38% for l = 4. Symbols colored in blue referto the (+1) excitation energies described in the main text. All quantities are dimensionless. See text for details. The figure isadapted from Ref. [187]. V ( r ) = (cid:40) C e − . r − R C ) if r ≤ R C C if r > R C (193)which is essentially the same as in Eq. (179) but without the central ring-shaped barrier. Here, the radius ofthe crater is set to R C = 3 . C = 200 . 0. Furthermore, the same Gaussianrepulsion as in Eq. (180) is employed.Excitations in this system are analyzed in the rotating frame of reference, where the Hamiltonian ˆ H in Eq.(2) contains an additional contribution accounting for the rotation around the z -axis. The whole Hamiltonianthus reads ˆ H rot = ˆ H − Ω rot ˆ L z (194)where Ω rot denotes the angular velocity around the z -axis and ˆ L z denotes the operator of the angular mo-mentum in the z -direction.A LR analysis is applied atop the ground states of N = 10 bosons with weak repulsion strength Λ = 0 . rot , and thus of the angular momentum per particle l , which is referred to as the vorticityor the charge below. The degree of condensation, N (1 − f ), ranges from 9 . 999 ( l = 0) to 9 . 962 ( l = 4) out of10 particles, meaning that the BEC is highly condensed for all considered parameter values. The upper panelsin Fig. 28 show the ground-state densities calculated with MCTDHB(7) for vorticities 0 ≤ l ≤ 4. Whereasfor l = 0, the density resembles a Gaussian in the center of the crater, the shape of the densities for l > l . The corresponding GP densities look alike (not shown).The lower panel in Fig. 28 compares the low-energy excitation spectrum calculated atop these groundstates, both for the mean-field BdG case (black and blue lines) and for LR-MCTDHB(7) (red and bluesquares). The excitations denoted by blue symbols are explained below. At first, the spectrum for l = 0 isexamined. Counting the number of excited states, one obtains 2 states from the BdG approach and 3 statesfrom the many-body approach. Moreover, the energies of the two mean-field excitations coincide with thecorresponding many-body results. Since the BdG theory only accounts for single-particle excitations, the thirdexcited stated obtained from LR-MCTDHB is an excitation where more than one particle is excited from theground state. In particular, it represents the state where two bosons are put into the orbital correspondingto the first excitation. The latter single-particle excitation is referred to as (+1) since it describes the caseof putting a boson from the ground-state mode with angular momentum l z = 0 to an orbital with l z = 1.Examining the spectra for l > 0, one observes two major changes. At first, the density of excitations in thelow-energy spectrum clearly grows with l , which can be deduced from the large number of many-body excitedstates entering. The number of mean-field excitations also grows (from 2 states for l = 0 to 4 states for l = 4),but very slightly compared to the number of many-body states. Therefore, the amount of excited states thatBdG misses is strongly growing with the vorticity of the BEC’s ground state. The second observation is thatthe values of the single-particle energies from BdG and LR-MCTDHB do not coincide any longer, as it wasthe case for l = 0. In fact, the entire BdG spectrum for vorticity l = 4 appears to be completely unrelated tothe corresponding many-body spectrum. The increasing gap between the (+1) results of the two approachesis clearly visible.Fig. 29 addresses the evolution of the energy of (+1), denoted by ω (+1) , with respect to the angular velocityΩ rot . Shown are the results obtained from the BdG theory, as well as many-body results for different numbersof bosons ( N = 10 , 100 and 1000). In the l = 0 sector, all energies coincide, indicating that for this case theBdG theory is sufficient. Henceforth, this obviously changes. With growing vortex charge of the ground state,the mean-field and many-body lines for (+1) separate further from each other. Most importantly, increasingthe number of bosons in the condensate only marginally changes this, meaning that the many-body andmean-field results for large N approach each other very slowly (see inset). As mentioned in the previoussection, it was proven analytically that for a trapped BEC with sufficiently weak and long-ranged repulsion,74 l=0 l=1 l=2 l=3 l=4 ω ( + ) Ω rot BdGN=10N=100N=1000 Figure 29: Excitation energies ω (+1) for different angular velocities Ω rot . Compared are the BdG energies and the many-bodyresults for different numbers of bosons N , obtained by LR-MCTDHB. The interaction strength is Λ = 0 . 5. The transition fromground-state vorticity l to l + 1 between two adjacent analyzed values of Ω rot are denoted by dotted vertical lines (mean-fieldin black and many-body in red). The higher the ground-state vorticity l , the larger is the separation between the BdG andmany-body energies. Inset: Magnified view for Ω rot = 3 . 8. All quantities are dimensionless. See text for details. The figure isadapted from Ref. [187]. the excitation energies in the Hartree limit converge towards the BdG energies [185]. The same is believed tohappen for a rotating BEC under certain conditions [206].In Fig. 30, the dependence of the difference ∆ N between the mean-field and many-body energies of (+1)on the number of bosons in the BEC is analyzed in more detail. The results are given relative to the energyfor N = 10 bosons, denoted by ∆ . For all shown vorticities, the gap becomes smaller with increasing N .However, it is obvious that even for BECs with an experimentally relevant number of bosons ( N = 1000)one is still far away from the predictions of the BdG theory. For l = 1, one obtains that the gap size ∆ is still about 97% of ∆ , and correspondingly higher for l > 1. Additionally, the shape of the exponentialfit curves to the numerical data does not show any sign of a steep descend when N is increased further. Itcan firstly be deduced that the excitation energies approach the BdG predictions for l > igure 30: Evolution of the energetic difference ∆ N between the mean-field and many-body results for ω (+1) for different particlenumbers N . The differences are shown relative to ∆ , i.e., the gap for N = 10 bosons. Solid lines denote exponential fits. Thegap becomes larger as l grows. Even for l = 1 and N = 1000, it is approximately 97% of ∆ , and the slopes of the fit curvesdo not indicate a sharp decent towards zero if N is increased further. All quantities are dimensionless. See text for details. Thefigure is taken from Ref. [187]. the ground state (see, e.g., Ref. [207]). Although the numerical results cannot give a definite answer for theHartree limit, they at least suggest that the regime in which one is clearly in need of an accurate many-bodytheory for the excitation energies covers the experimentally relevant regime of rotating BECs, with N beingof the order of 10 − bosons. A further discussion on the applicability of the GP and BdG theories in theHartree limit can be found in Appendix A.It is finally suggested, similar to the 1D systems discussed in Section 6.3.1, how a quench can be utilizedin order to excite breathing modes. To this end, the ground state of a BEC with N = 10 bosons and angularvelocity Ω rot = 3 . R C = 3 . 01, iscalculated. The obtained ground state has the same symmetry and vorticity l = 3 as the one for R C = 3 . rot = 3 . 2. Then, the crater is suddenly shrinked back to R C = 3 . N ∆ X ( t ) is carried out. Due to theradial symmetry, one could as well analyze the position variance in the y -direction. The result are shown inFig. 31. The many-body spectrum shows two distinct peaks, where the first corresponds to the (1 +1 , − )two-particle excitation, where one boson is excited to an orbital with l = 4 and one is excited to an orbital76 Ω rot =3.2 Λ =0.5 Intensity ω M=7GP (1 +1 ,1 -1 ) at ω =0.605Deexcitation: ω =10.6 -> 9.47 ∆ /N t Figure 31: Low-energy Fourier spectrum of the time evolution of the position variance in x -direction, N ∆ X , due to a quantumquench of the radius R C from 3 . 01 to 3 . 0. The number of bosons is N = 10. The many-body spectrum (red), obtained fromthe evolution with M = 7 orbitals, shows two sharp peaks in the considered energy range, whereas the GP spectrum (black) isflat. The first peak refers to the two-particle excitation (1 +1 , − ), whereas the second peak indicates a de-excitation betweentwo higher-lying states (see right inset). The left inset shows the evolution of N ∆ X from 0 ≤ t ≤ 30. The depletion of the BECstays below 0 . 4% during the entire propagation. All quantities are dimensionless. See text for details. with l = 2, such that the total angular momentum is unchanged. The energy that LR-MCTDHB(7) predictsfor this state is ω = 0 . 605 which coincides with the peak position of the Fourier analysis. The second peakat ω ≈ . 12 denotes a de-excitation between two states that are much higher in energy ( ω ≈ . → . ≤ ω ≤ 11 is shown. One sees that in the GP case, there isonly one peak at ω ≈ . 6, but no second one at ω ≈ . 47 as in the many-body spectrum. Thus, very differentspectra from the dynamics of the variance at the mean-field and many-body levels are obtained, and it isshown that the suggested quench leads to an excitation of a many-body breathing mode. It is stressed thatthe depletion stays below 0 . 4% for the entire propagation, which means that the BEC is essentially condensed.To summarize, the low-energy spectrum of a rotating BEC in an anharmonic, radially-symmetric cratershows significant many-body effects, even for weakly-interacting, highly-condensed systems. These effects can77e observed already for the lowest-in-energy, single-particle excited states, and remain even for mesoscopicand large BECs. Furthermore, a simple quench triggers a many-body breathing mode, which can be observed,e.g., by analyzing the Fourier spectrum of the position variance. This review presents recently developed many-body approaches capable to describe and compute accuratelythe ground state and dynamics as well as the energy spectrum of trapped ultracold bosons. It is relevant toemploy methods which are on the same footing. The dynamics is calculated by the well-known MCTDHBmethod, the ground state is obtained by imaginary-time propagation using the same method, and the energyspectrum of the system is then determined by applying the linear-response theory to the ground state computedvia MCTDHB. The resulting theory is named LR-MCTDHB.The MCTDHB theory is briefly discussed and more attention is paid to the less spread LR-MCTDHBtheory. As the implementation of the MCTDHB is well documented in the literature, only references are givenhere, while the newly developed LR-MCTDHB implementation which makes the theory widely applicable, ispresented here in some detail for the first time. The latter posed a major challenge due to both the complexityof the underlying theory and the accompanying technical difficulties like the vast memory consumption andthe necessity of parallelization. The code structure is explained in detail both with respect to the constructionof the linear-response matrix and its subsequent partial diagonalization which provides the desired excitationenergies. A comprehensive description of the procedure used to calculate the lowest eigenvalues is also pro-vided. The corresponding benchmarks of the LR-MCTDHB code (see also Appendix C) with an analyticallysolvable many-boson non-trivial model demonstrate its correct functioning and show that the theory providesnumerically exact results.The second main goal of this review is to present applications to many-body dynamics and excitations ofBECs in interesting scenarios, and, of course, to discuss the underlying physics. Several systems of trappedBECs in one and two dimensions are addressed in detail and the main attention is drawn to the detectionof many-body effects that are not describable by applying mean-field theory. Particularly emphasized arethe tunneling dynamics in double-well systems and the development of fragmentation as a function of timein initially coherent BECs. The applications shown in this work have demonstrated that a full many-bodytreatment of the dynamics of trapped interacting BECs is inevitable. For instance, it is exhibited that theout-of-equilibrium tunneling dynamics of initially coherent BECs in 1D and 2D double wells are of many-bodynature. After a few Rabi cycles, tunneling between the wells becomes suppressed at the many-body level,whereas unperturbed oscillations are predicted at the mean-field level. Surprisingly, the degree of fragmenta-tion in 1D is universal for constant values of the interaction parameter. Going to 2D, the sensitivity of thetunneling dynamics to many-body effects is even enhanced as soon as the bosons carry angular momentum.This can also be deduced from the larger amount of orbitals necessary in the limit of weak repulsion as com-pared to the case of bosons with zero angular momentum. Further applications deal with the nucleation ofphantom vortices, which are pure many-body objects not observable in the condensate density, as well as withthe tunneling of bosons to open space. In both cases, the initially coherent systems become fragmented andare, therefore, in need of an accurate many-body description.With respect to excitations, it is observed that the ground-state depletion has a strong impact on the low-energy spectrum. Even for marginal depletion where one might expect mean-field approaches to give accurateresults, substantial many-body effects like the shift of excitation energies and the appearance of complicatedmany-boson excitations occur. Applications in 1D explores the many-body effects in double- and triple-wellsystems as well as in larger lattices. In 2D, it is found that once the bosons carry angular momentum and theground state becomes a vortex, the deviations between the mean-field and many-body excitation spectra growstrongly. Moreover, this discrepancy tends to persist when the number of particles in the condensate increasesup to an experimentally relevant size. This has a potentially large impact on experiments on ultracold bosonssince the dynamics are mostly determined by the lowest-in-energy excitations. Hopefully, these findings and78he new insights into the physics of ultracold bosons will stimulate future experiments.Before closing, we would like to add several remarks concerning further general and technical developmentsof the methods which can further broaden their range of applicability. Until now the applications of the many-body linear-response theory were for excitations on top of the system’s ground state. We emphasize that theLR-MCTDHB can also be applied to any other state obtained via MCTDHB. In other words, it is possible tocompute excitations on top of any other state. This will help to address excitations which are difficult to reachfrom the ground state. Another general possibility lies in the formulation of LR-MCTDHB in block-diagonalform, which is potentially advantageous in calculating the low-energy spectrum, as discussed in Appendix B.Apart from these, further general extensions and developments of both the LR theory itself as well as ofits numerical implementation are worth of additional investigation. With regard to theoretical developments,an extension of LR-MCTDHB to BEC mixtures, i.e., compound systems of different types of bosons, is highlydesirable as their popularity has substantially grown [208, 184, 209, 210, 211, 212, 213]. The underlyingMCTDHB theory for bosonic mixtures, also including the possibility of internal degrees of freedom, has beenformulated [53, 70, 71], but a linear-response theory is completely missing. It turns out that the so-calledmulti-layer MCTDHB, introduced and described in Ref. [70], is particularly advantageous numerically. Tothis end, an interesting goal would be to apply the linear-response theory to multi-layer MCTDHB. Moreover,mixtures of bosons and fermions might become even more relevant [71, 214, 218, 215, 216, 217, 219], meaningthat also the equations of motion of the fermionic counterpart of MCTDHB, called MCTDHF (see, e.g., Ref.[220] and references therein, and the recently developed implementation described in Ref. [221]), should beutilized to develop the respective linear-response theory.Further technical developments which accelerate the computations are, of course, also welcome. Oneinteresting possible extension lies in the application of linear-response theory to the time-dependent restricted-active-space self-consistent-field (TD-RASSCF) method, introduced in Ref. [222]. The latter theory is similarto MCTDHB, but is an approximation to it. It is, however, capable of reducing systematically the amountof necessary coefficients significantly. This is achieved by restricting the possible excitations in the activespace of the single-particle orbitals according to certain conditions or symmetries of the physical problem athand. Thus, the dimensionalities of the linear-response matrices which would result when constructed basedon the linear-response theory on top of TD-RASSCF, could also be greatly reduced. Such an approach willsubstantially simplify the computation of the excitation spectrum and, moreover, enable the treatment ofeven larger systems, i.e., with more particles and if necessary higher numbers of orbitals.With regard to the technical development of the current implementation, the application of additionaltechniques to diagonalize the often huge linear-response matrix can be beneficial. A potentially very powerfulmethod may be the FEAST algorithm [224], which is a contour integration and density matrix-based method.It is available in a free software package [225] and also allows for high-level parallelization. In particular, itenables the specification of an energy interval in which the desired eigenvalues should be found. Hence, itmight become easier to explore not just the lowest but also intermediate eigenvalues.The above outlined avenues for future research demonstrate the large variety of perspectives, which thiswork will hopefully stimulate. Appendix A Mean-field theory In this Appendix, several mean-field approaches to describe the ground state and dynamics as well as theexcitation spectra of trapped BECs are presented. At first, the GP single-orbital mean-field theory (SectionA.1.1) and the LR theory atop it, yielding the famous BdG equations (Section A.1.2), are described. Thoseapproaches are the most common ones in the literature, and are expected to accurately describe the physics ofcoherent BECs. They are also used in this work to compare mean-field and many-body results. Afterwards,a time-dependent multi-orbital mean-field approach, called TDMF (Section A.2.1), and its corresponding LRtheory LR-BMF (Section. A.2.2), are discussed. The latter two theories are suitable for describing fragmentedBECs at the mean-field level. They can be understood as bridging theories between the GP/BdG theories79nd the full many-body descriptions given by MCTDHB and LR-MCTDHB, respectively. A.1 Single-orbital approaches A.1.1 The GP equation The GP equation is commonly used to obtain the ground state and the time evolution of a coherent BEC,i.e., where all N bosons occupy the same single-particle state φ ( r , t ). The many-boson wave function Ψ isthus given by a single Hartree product,Ψ( r , ..., r N ; t ) = φ ( r , t ) ... φ ( r N , t ) ≡ | N ; t (cid:105) . (A.1)Originally, the GP equation is derived assuming contact interaction between the bosons. However, to be moregeneral, the two-body interaction potential is taken to be of general form ˆ W ( r , r (cid:48) ).The GP equation can be derived from the least action principle defined in Eq. (25). The expectation valuethat appears in the action functional, one obtains (cid:28) N ; t (cid:12)(cid:12)(cid:12)(cid:12) ˆ H − i ∂∂ t (cid:12)(cid:12)(cid:12)(cid:12) N ; t (cid:29) = N h GP + λ N ( N − W GP − iN (cid:18) ∂∂t (cid:19) GP (A.2)where the quantities h GP and W GP are defined in Eqs. (21) and (22) utilizing φ ( r , t ). The time derivative isgiven by (cid:0) ∂∂t (cid:1) GP = (cid:82) φ ∗ ( r , t ) ∂∂t φ ( r , t ) d r . One arrives at0 = δSδφ ∗ ( r , t ) (cid:20) N h GP + λ N ( N − W GP − iN (cid:18) ∂∂t (cid:19) GP (cid:21) = N (cid:20) ˆ h + λ ( N − (cid:90) ˆ W ( r , r (cid:48) ) | φ ( r (cid:48) , t ) | d r (cid:48) − i ∂∂t (cid:21) φ ( r , t ) . (A.3)Thus, the time-dependent GP equation used in this work is thus given by (cid:20) ˆ h + Λ (cid:90) ˆ W ( r , r (cid:48) ) | φ ( r (cid:48) , t ) | d r (cid:48) (cid:21) φ ( r , t ) = i ∂∂t φ ( r , t ) (A.4)with the mean-field interaction parameter Λ = λ ( N − (cid:20) ˆ h + Λ (cid:90) ˆ W ( r , r (cid:48) ) | φ ( r (cid:48) ) | d r (cid:48) (cid:21) φ ( r ) = µφ ( r ) (A.5)where µ is the chemical potential. Whereas Eq. (A.4) describes the time-evolution of a fully-condensed BEC,Eq. (A.5) is used to determine the GP ground-state orbital φ ( r ).This section is closed with a brief discussion of the applicability of the GP theory in the Hartree limitof repulsive trapped BECs, meaning for N → ∞ with fixed interaction parameter Λ. It has been provenanalytically that in this limit, the GP theory provides the exact energy and density per particle of the system’sground state [226]. This, however, does not necessarily mean that the corresponding GP wave function, Eq.(A.1), coincides with the exact wave function in this limit. To this end, it has been demonstrated recentlythat the overlap of the exact and GP wave function in the Hartree limit can be clearly less than unity or evenbe vanishingly small [227]. Moreover, a subsequent study has shown the many-body character of the exactwave function in the Hartree limit by applying many-body perturbation theory to the mean-field Hamiltonian[228], which demonstrates that the mean-field and exact wave functions can be very different from each other.Moreover, as a further interesting result, it was found that the condensate depletion in the Hartree limitbecomes a constant. These findings potentially have severe consequences for the excitation spectrum of aBEC in or close to the Hartree limit, which is discussed at the end of Section A.1.2.80 .1.2 The BdG equations The LR theory atop the time-dependent GP equation, Eq. (A.4), is derived below. Therefore, a weaktime-dependent perturbation is added to the Hamiltonian, i.e., ˆ H ( r ) → ˆ H ( r ) + ˆ H ext ( r , t ). The perturbationˆ H ext ( r , t ) explicitly reads ˆ H ext ( r , t ) = f + ( r ) e − iωt + f − ( r ) e iωt (A.6)with the real amplitudes f + ( r ) and f − ( r ) and the frequency ω . The response of the system is calculatedaround the stationary GP ground-state orbital φ ( r ) obtained from Eq. (A.5). For the perturbed orbital φ ( r , t ), the ansatz √ N φ ( r , t ) = e − iµt (cid:16) √ N φ ( r ) + u ( r ) e − iωt + v ∗ ( r ) e iωt (cid:17) (A.7)is utilized. In Eq. (A.7), the quantities u ( r ) and v ( r ) denote small and time-independent response amplitudes.Plugging in Eqs. (A.6) and (A.7) into the time-dependent GP equation, Eq. (A.4), yields in zeroth order thetime-independent GP equation, Eq. (A.5). In first order, collecting all terms proportional to e − iωt leads to( ω + µ ) u ( r ) = √ N f + ( r ) φ ( r ) + ˆ H GP u ( r ) + Λ (cid:90) ˆ W ( r , r (cid:48) ) φ ( r (cid:48) ) v ( r (cid:48) ) φ ( r ) d r (cid:48) + Λ (cid:90) ˆ W ( r , r (cid:48) ) φ , ∗ ( r (cid:48) ) u ( r (cid:48) ) φ ( r ) d r (cid:48) (A.8)where the GP Hamiltonian ˆ H GP is defined asˆ H GP = ˆ h + Λ (cid:90) ˆ W ( r , r (cid:48) ) | φ ( r (cid:48) , t ) | d r (cid:48) . (A.9)Collecting all terms proportional to e iωt results in( µ − ω ) v ∗ ( r ) = √ N f − ( r ) φ ( r ) + ˆ H GP v ∗ ( r ) + Λ (cid:90) ˆ W ( r , r (cid:48) ) φ ( r (cid:48) ) u ∗ ( r (cid:48) ) φ ( r ) d r (cid:48) + Λ (cid:90) ˆ W ( r , r (cid:48) ) φ , ∗ ( r (cid:48) ) v ∗ ( r (cid:48) ) φ ( r ) d r (cid:48) . (A.10)Complex conjugation and multiplication by ( − 1) of Eq. (A.10) opens the possibility to cast Eqs. (A.8) and(A.10) in matrix form. The result reads( L BdG − ω ) (cid:18) u ( r ) v ( r ) (cid:19) = (cid:18) −√ N f + ( r ) φ ( r ) √ N f − , ∗ ( r ) φ , ∗ ( r ) (cid:19) (A.11)with the BdG matrix L BdG = (cid:18) A B − B ∗ − A ∗ (cid:19) (A.12)where the quantities A and B are given by A = ˆ H GP − µ + Λ (cid:90) d r (cid:48) ˆ W ( r , r (cid:48) ) φ , ∗ ( r (cid:48) ) φ ( r ) ˆ P rr (cid:48) (A.13)and B = Λ (cid:90) d r (cid:48) ˆ W ( r , r (cid:48) ) φ ( r (cid:48) ) φ ( r ) ˆ P rr (cid:48) . (A.14)where the operator ˆ P rr (cid:48) was introduced in Eq. (96). In the literature, Eq. (A.11) with contact interactionˆ W ( r , r (cid:48) ) = δ ( r − r (cid:48) ) and zero RHS, i.e., ( L BdG − ω ) (cid:18) u ( r ) v ( r ) (cid:19) = 0 , (A.15)81s referred to as the BdG equations [23, 24]. In this work, however, the interaction potential is in generaldifferent from the contact interaction potential. Moreover, the particle-conserving BdG equations [229, 230,231, 232], given by PL BdG P (cid:18) u ( r ) v ( r ) (cid:19) = ω (cid:18) u ( r ) v ( r ) (cid:19) (A.16)where the projector P = − | φ (cid:105)(cid:104) φ | removes any contribution of the ground-state orbital φ ( r ) from theresponse amplitudes u ( r ) and v ( r ), are used to calculate the low-energy excitation spectra at the mean-fieldlevel.From the above derivation one can readily see the close connection of the GP orbital and the BdG equations.As discussed in the previous section on the GP theory, the wave function of the system in the Hartree limitcan be very different from the GP wave function, although the latter gives the exact energy and density perparticle. Thus, it can be questioned whether the LR theory atop the GP wave function, resulting in the BdGequations, is in general reliable in the Hartree limit. Scientific works that prove the exactness of the BdGspectrum in this limit under certain conditions were discussed in Section 6.3.3. Nevertheless, the example ofa rotating BEC in an anharmonic trap, which does not fulfill the assumptions made in Ref. [206], serves asan example where at least in the experimentally relevant regime of N ≈ − bosons in the condensatethe lowest-in-energy excitations substantially deviate from the BdG predictions. A.2 Multi-orbital approaches A.2.1 Time-dependent multi-orbital mean field (TDMF) The TDMF theory, introduced in Ref. [233], is briefly described. It can be seen as an intermediate stepbetween the GP theory and MCTDHB because it allows for more than one single-particle orbital that thebosons can occupy, and is thus applicable to fragmented condensates. However, as a major difference to thefull many-body description in Section 3, the occupation of these single-particle orbitals is fixed. The ansatzfor the N -boson wave function is given byΨ( r , ..., r N ; t ) = ˆ S φ ( r , t ) φ ( r , t ) ... φ N ( r N , t ) (A.17)where the operator ˆ S accounts for the symmetrization of the wave function. In general, not all single-particleorbitals φ i need to be different. The only constraint is that the number of particles is fixed such that N = (cid:80) Mi =1 n i where M is the number of different orbitals in Eq. (A.17) and the n i are their fixed occupationnumbers. Applying the least action principle to S = (cid:90) dt (cid:104) Ψ | ˆ H − i ∂∂t | Ψ (cid:105) − M (cid:88) i,j =1 µ ij ( t )[ (cid:104) φ i | φ j (cid:105) − δ ij ] (A.18)with the time-dependent Lagrange multipliers µ ij ( t ) yields0 = δSδφ ∗ k ( r , t ) = n k ˆ h − i ∂∂t + λ ( n k − 1) ˆ J k + M (cid:88) l (cid:54) = k λ n l ( ˆ J l + ˆ K l ) | φ k (cid:105) = M (cid:88) j =1 µ kj ( t ) | φ j ( t ) (cid:105) ∀ k = 1 , ..., M (A.19)where ˆ J l ( r , t ) = (cid:90) d r (cid:48) φ ∗ l ( r (cid:48) , t ) ˆ W ( r , r (cid:48) ) φ l ( r (cid:48) , t ) (A.20)82nd ˆ K l ( r , t ) = (cid:90) d r (cid:48) φ ∗ l ( r (cid:48) , t ) ˆ W ( r , r (cid:48) ) ˆ P rr (cid:48) φ l ( r (cid:48) , t ) (A.21)are the direct and exchange interaction potentials, respectively. Multiplying Eq. (A.19) by (cid:104) φ l | gives µ kl ( t ) = n k h lk − (cid:18) ∂∂t (cid:19) lk + λ ( n k − 1) ˆ J k W lkkk + M (cid:88) j (cid:54) = k λ n j ( W ljjk + W ljkj ) (A.22)for the Lagrange multipliers. Plugging this into Eq. (A.19) results in i ˆP ∂∂t | φ k (cid:105) = ˆP ˆ h + λ ( n k − 1) ˆ J k + M (cid:88) l (cid:54) = k λ n l ( ˆ J l + ˆ K l ) | φ k (cid:105) ∀ k = 1 , ..., M (A.23)where the projector ˆP = − (cid:80) Mj =1 | φ j (cid:105)(cid:104) φ j | projects onto the tangential space of span( φ , ..., φ M ). To simplifythe above EOMs for the orbitals, the condition (cid:104) φ i | ˙ φ j (cid:105) = 0 ∀ i, j = 1 , ..., M (A.24)is invoked. This leads to the conservation of energy, and the projector on the LHS of Eq. (A.23) can beomitted, i.e., i ∂∂t | φ k (cid:105) = ˆP ˆ h + λ ( n k − 1) ˆ J k + M (cid:88) l (cid:54) = k λ n l ( ˆ J l + ˆ K l ) | φ k (cid:105) ∀ k = 1 , ..., M. (A.25)It is important to note that the condition in Eq. (A.24) is a manually-invoked restriction to the orbitals withinthe TDMF theory, whereas the same condition in MCTDHB appears naturally and rigorously. Furthermore,one observes that for M = 1, i.e., when only a single orbital is available, Eq. (A.25) reduces to the (projected)time-dependent GP equation, Eq. (A.4). The TDMF theory has been extended to describe spinor condensatesand Bose-Bose or Bose-Fermi mixtures as well, see in this context Ref. [234]. By imaginary-time propagationof the TDMF EOMs, one obtains the corresponding theory for the stationary equations, called the best meanfield for condensates [235]. Applications can be found in, e.g., Refs. [236, 237, 238, 239]. A.2.2 Linear-response best-mean-field (LR-BMF) A possible LR approach based on a multi-orbital mean-field description is the LR-BMF theory which hasbeen introduced in Ref. [240]. To derive it, the EOMs of the previously described TDMF theory, Eq.(A.19), are linearized for the case of adding a weak time-dependent external field to the Hamiltonian, i.e.,ˆ H → ˆ H + ˆ H ext ( t ). The perturbation ˆ H ext ( t ) is assumed to be of the same form as in Eq. (A.6). Pluggingthis into Eq. (A.19) yields (cid:18) ˆ Z i − i ∂∂t (cid:19) | φ i (cid:105) − M (cid:88) j =1 µ ij ( t ) n i | φ j (cid:105) = − ˆ H ext ( t ) | φ i (cid:105) (A.26)with the replacement ˆ Z i = ˆ h + λ ( n i − 1) ˆ J i + M (cid:88) l (cid:54) = i λ n l ( ˆ J l + ˆ K l ) . (A.27)83he Lagrange multipliers, after redefining them as µ ij ( t ) /n i → µ ij ( t ), can be expressed by µ ij ( t ) = (cid:104) φ j | ˆ Z i + ˆ H ext ( t ) | φ i (cid:105) (A.28)where the condition of Eq. (A.24) was used. The ansatz for the perturbed orbitals is given by φ i ( r , t ) ≈ φ i ( r ) + δφ i ( r , t ) (A.29)Plugging Eq. (A.29) into Eq. (A.26) leads to − ˆ H ext ( t ) | φ i (cid:105) = (cid:18) ˆ Z i − i ∂∂t (cid:19) | δφ i (cid:105) − M (cid:88) j =1 µ ij ( t ) | δφ j (cid:105) + λ ( n i − δ ˆ J i | φ i (cid:105) + M (cid:88) j (cid:54) = i λ n j ( δ ˆ J j + δ ˆ K j ) | φ i (cid:105) (A.30) − M (cid:88) j =1 δµ ij ( t ) | φ j (cid:105) with δ ˆ J i ( r , t ) = (cid:90) d r (cid:48) δφ ∗ i ( r (cid:48) , t ) ˆ W ( r , r (cid:48) ) φ i ( r (cid:48) ) + (cid:90) d r (cid:48) φ , ∗ i ( r (cid:48) ) ˆ W ( r , r (cid:48) ) δφ i ( r (cid:48) , t ) (A.31)and δ ˆ K i ( r , t ) = (cid:90) d r (cid:48) δφ ∗ i ( r (cid:48) , t ) ˆ W ( r , r (cid:48) ) ˆ P rr (cid:48) φ i ( r (cid:48) ) + (cid:90) d r (cid:48) φ , ∗ i ( r (cid:48) ) ˆ W ( r , r (cid:48) ) ˆ P rr (cid:48) δφ i ( r (cid:48) , t ) . (A.32)The variation of the Lagrange multipliers can be expressed as δµ ij ( t ) = − M (cid:88) l =1 µ il (cid:104) φ j | δφ l (cid:105) + (cid:104) φ j | δ ( ˆ Z i | φ i (cid:105) ) + (cid:104) φ j | ˆ H ext ( t ) | φ i (cid:105) (A.33)where the identity ˆ Z i | φ i (cid:105) = (cid:80) Ml =1 µ il | φ l (cid:105) and integration by parts was used in the first term. One finallyobtains ˆP ˆ Z i | δφ i (cid:105) − M (cid:88) j =1 µ ij | δφ j (cid:105) + λ ( n i − δ ˆ J i + M (cid:88) j (cid:54) = i λ n j ( δ ˆ J j + δ ˆ K j ) | φ i (cid:105) = − ˆP ˆ H ext ( t ) | φ i (cid:105) + i ∂∂t | δφ i (cid:105) (A.34)with the projector ˆP = − (cid:80) k | φ k (cid:105)(cid:104) φ k | . By first making the ansatz δφ i ( r , t ) = 1 √ n i (cid:0) u i ( r ) e − iωt + v ∗ i ( r ) e iωt (cid:1) (A.35)with stationary response amplitudes u i ( r ) and v i ( r ) and plugging this into Eq. (A.34) afterwards, it is againpossible, similar to derivation of the BdG equations in Section A.1.2, to group the terms of the resultingequations with respect to their linearity to either e − iωt or e iωt . As a result, one arrives at an eigenvalueequation of the form ( PL − ω ) (cid:18) | u (cid:105)| v (cid:105) (cid:19) = P (cid:18) − f + ( r ) | φ n (cid:105) f − , ∗ ( r ) | φ , ∗ n (cid:105) (cid:19) (A.36)84ith the vector | φ n (cid:105) = |√ n φ , ..., √ n M φ M (cid:105) , the 2( M × M ) square matrix L = (cid:18) Z − µ + A B − B ∗ − ( Z − µ + A ) ∗ (cid:19) , (A.37)and the projection matrix P = {P ij } = ˆP i = j ≤ M ˆP ∗ i = j > M i (cid:54) = j (A.38)of the same size. The matrix Z = diag( ˆ Z , ..., ˆ Z M ) is a diagonal matrix, whereas µ = { µ ij } contains theLagrange multipliers. The matrices A and B , for the case of the contact interaction potential, read A = { A ij } = (cid:40) λ ( n i − | φ i | i = j λ √ n i n j φ i φ , ∗ j i (cid:54) = j (A.39)and B = { B ij } = (cid:40) λ ( n i − φ i ) i = j λ √ n i n j φ i φ j i (cid:54) = j. (A.40)Since the term proportional to ω is the only one without a projector in Eq. (A.36), a redundant projectorcan be added to the RHS of PL , i.e., PL → PLP . The resulting homogeneous eigenvalue equation PLP (cid:18) | u (cid:105)| v (cid:105) (cid:19) = ω (cid:18) | u (cid:105)| v (cid:105) (cid:19) (A.41)is the central equation in the LR-BMF theory. A comparison to the BdG matrix L BdG in Eqs. (A.12)-(A.14)shows that for M = 1, LR-BMF reduces to the BdG theory, i.e., LR-BMF( M = 1) ≡ BdG. If one howevercompares the LR matrix of Eq. (A.37) to the one of LR-MCTDHB in Section 3.2.2, it can be seen thatthe latter is clearly more complicated because it also includes the couplings between orbitals and coefficients.Again, in both mean-field approaches presented in this Appendix, there is only one possible configuration, i.e.,the configuration | N ; t (cid:105) for the BdG case and the configuration | n , ..., n M ; t (cid:105) for the multi-orbital mean-fieldcase. The occupations of the underlying single-particle orbitals are fixed. In fact, the LR-BMF matrix PLP can be seen as the upper-left block L oo of the LR matrix in Eq. (93). Further details on the LR-BMF theory,together with an application to a BEC in a double well, can be found in Ref. [240]. Appendix B LR-MCTDHB in block-diagonal form In this Appendix, a complex transformation that essentially halves the dimensionality of the homogeneouseigenvalue problem of Eq. (116) is introduced. However, as shown below, it involves matrix-matrix multipli-cations of individual blocks of the LR matrix L . The latter products might be numerically expensive. It istherefore not in general beneficial to perform the transformation described below, which was first discussedin Ref. [241].One considers the transformation matrix Q = 1 √ (cid:18) θ − θ (cid:19) , Q − = 1 √ (cid:18) θ − θ (cid:19) , QQ − = Q − Q = (B.1)where the operation θ is defined via its action on any operator ˆ O given by θ ˆ O = ˆ O ∗ . A matrix of the form A = (cid:18) A B − B ∗ − A ∗ (cid:19) (B.2)85ehaves under the transformation in Eq. (B.1) as Q A Q − = 12 (cid:18) θ − θ (cid:19) (cid:18) A B − B ∗ − A ∗ (cid:19) (cid:18) θ − θ (cid:19) = (cid:18) A − BθA + Bθ (cid:19) . (B.3)Multiplying the result of Eq. (B.3) with itself yields the block-diagonal matrix (cid:0) Q A Q − (cid:1) (cid:0) Q A Q − (cid:1) = (cid:18) ( A − Bθ )( A + Bθ ) 00 ( A + Bθ )( A − Bθ ) (cid:19) . (B.4)This transforms the eigenvalue problem A α = λα into λ ( Qα ) = (cid:0) Q A Q − (cid:1) (cid:0) Q A Q − (cid:1) ( Qα ) ⇔ λ √ (cid:18) u + v ∗ u − v ∗ (cid:19) = (cid:18) ( A − Bθ )( A + Bθ ) 00 ( A + Bθ )( A − Bθ ) (cid:19) √ (cid:18) u + v ∗ u − v ∗ (cid:19) (B.5)with α = √ (cid:18) uv (cid:19) . Thus, the original eigenvalue problem of A can be split up into two smaller eigenvalueproblems of L (1) ≡ ( A − Bθ )( A + Bθ ) and L (2) ≡ ( A + Bθ )( A − Bθ ) with dim( L (1) ) = dim( L (2) ) = dim( A ) / L (1) and L (2) have the same eigenvalues, which are the squared eigenvalues of A .Interchanging the second and third row as well as the second and third column of the LR matrix in Eq.(115) which explicitly reads L = ρ − ˆ P L uoo ˆ P ρ − ρ − ˆ P L voo ˆ P ∗ ( ρ ∗ ) − ρ − ˆ P L uoc ρ − ˆ P L voc − ( ρ ∗ ) − ˆ P ∗ L v, ∗ oo ˆ P ρ − − ( ρ ∗ ) − ˆ P ∗ L u, ∗ oo ˆ P ∗ ( ρ ∗ ) − − ( ρ ∗ ) − ˆ P ∗ L v, ∗ oc − ( ρ ∗ ) − ˆ P ∗ L u, ∗ oc L uco ˆ P ρ − L vco ˆ P ∗ ( ρ ∗ ) − H − ε −L v, ∗ co ˆ P ρ − −L u, ∗ co ˆ P ∗ ( ρ ∗ ) − − ( H − ε ) ∗ (B.6) yields a matrix of the form given in Eq. B.2. Its submatrices are A = (cid:18) ρ − ˆ P L uoo ˆ P ρ − ρ − ˆ P L uoc L uco ˆ P ρ − H − ε (cid:19) (B.7)and B = (cid:18) ρ − ˆ P L voo ˆ P ∗ ( ρ ∗ ) − ρ − ˆ P L voc L vco ˆ P ∗ ( ρ ∗ ) − (cid:19) . (B.8)Mind the suppressed superscript ’0’ for the density matrices as compared to, e.g., Eq. (111). It is hencepossible to diagonalize either L (1) or L (2) with A and B defined in Eqs. (B.7) and (B.8) to obtain the squaredeigenvalues of the original LR matrix L . The advantage is, as mentioned above, that the dimensionality of L (1) and L (2) is only half of the dimensionality of L . Moreover, for two distinct eigenvalues ω and ω of L with ω ≥ ω > . ≡ ω − ω ≥ 0, the separation ∆ (2) in the spectrum of the squaredeigenvalues is given by ∆ (2) ≡ ω − ω = ( ω + ∆) − ω = 2 ω ∆ + ∆ = ∆(2 ω + ∆) ≥ ω ∆ > ∆ . (B.9)Thus, the separation grows in the spectra of L (1) and L (2) . Since in general, the IRAM converges in lessiterations when the eigenvalues are more separated from each other, eigenvalues of larger magnitude than 0 . L (1) or L (2) instead of L . On the contrary, foreigenvalues with 0 ≤ ω ≤ ω ≤ . 5, one finds that∆ (2) = ∆(2 ω + ∆)= ∆ ( ω + ω ) (cid:124) (cid:123)(cid:122) (cid:125) ≤ ≤ ∆ , (B.10)which means that the separation of these eigenvalues is decreased in the spectra of L (1) and L (2) . Therefore,although each iteration of the IRAM is less expensive for the smaller matrices, more iterations might benecessary to obtain converged eigenvalues with smaller magnitude than 0.5. A possible solution to thisproblem is to shift L → L + 0 . · , where is the 2( M + N conf )-dimensional unit matrix, before applying thetransformation in Eq. (B.1).It therefore remains a system-dependent question whether it is beneficial to use the full LR matrix L orone of the smaller matrices L (1) and L (2) to obtain the low-energy spectrum. One needs to take into accountthe computational cost of building the latter two matrices. Moreover, it is not in general true that, assumingthe blocks A and B are sparse matrices, their product is also sparse. If one is in the unfortunate situationthat the matrices L (1) and L (2) are much denser than the original matrix L , it is potentially also much moredemanding to converge to their lowest eigenvalues, although they are much smaller. It is therefore a priori notpossible to decide which approach is more efficient, and in principle, one would always need to first calculatethe (shifted) matrices L (1) and L (2) and analyze how dense they are before making an appropriate decision.All applications of LR-MCTDHB discussed in this work, however, avoided this preliminary analysis, and usedthe full LR matrix L to obtain the low-energy spectra. For additional details, it is referred to Ref. [241].87 ppendix C Further benchmarks of LR-MCTDHB Further results on the benchmark of the newly developed implementation of LR-MCTDHB against the HIMare presented in this Appendix. In particular, comparisons of the exact and numerical results for attractivebosons in 1D as well as for repulsive bosons in the anisotropic HIM in 2D are made. In addition, the numericalconvergence in the rotating frame of reference, which is important for the validity of the results presented inSection 6.3.3, is demonstrated.For the case of attractive bosons in the 1D HIM, the system parameters are chosen as Ω = 1 . λ = 0 . 13 for the two-body interaction strength. The number of bosons is N = 10.The calculations were carried out on a grid with 128 grid points in the interval [ − , ω i = E i − ε , are shown in Table C.1. With respectto ε , one observes that M = 4 and M = 6 lead to highly accurate results. With respect to the excitedstates, it can be observed that the BdG approach, i.e., for M = 1, misses several excitations in the low-energyspectrum. In contrast to that, the many-body calculations yield all states. The numerical accuracy clearlyimproves upon increasing the number of self-consistent orbitals. For M = 6, even the higher c.m. excitedstates ω , ω , and ω are obtained to very high accuracy. Together with the benchmark from Table 2, itcan be deduced that the implementation of LR-MCTDHB is capable of accurately describing the low-energyspectrum of both repulsive and attractive bosons.Table C.2 shows the numerical results obtained for the anisotropic HIM in 2D with trap frequenciesΩ x = 1 . y = √ . 9. The remaining system parameters are the same as in the isotropic case from Table3. One first observes that the anisotropy lifts all previous degeneracies of the isotropic case. Apart from that,the numerical results are qualitatively similar. In the BdG case, only the first c.m. excitation is obtainednumerically exact, all other c.m. excited states are inaccessible. The pure relative excitations are also foundto a good accuracy. For M = 2, one gets access to higher c.m. excited states in the x -direction, together withcombined excitations of that type with excitations of relative coordinates. As in the isotropic case, all excitedstates obtained in the y -coordinate deteriorate in accuracy compared to the BdG case, see, e.g., ω and ω .The reason is the same as for the isotropic case, namely that the additional orbital resembles a p x -orbital andleads to preferred direction in the description. Including a third orbital that resembles a p y -orbital solves thisproblem. It additionally improves the overall accuracy of the low-energy excitations, and moreover ensuresthat all low-lying excited states are obtained. M = 1 M = 4 M = 6 ( m x , n x ) Exact ε , 0) 9.038150 ω , 0) 1.000000 ω n/a 2.000222 2.000000 (2 , 0) 2.000000 ω n/a 3.000432 3.000000 (3 , 0) 3.000000 ω , 2) 3.794733 ω n/a 4.011839 4.000012 (4 , 0) 4.000000 ω n/a 4.794870 4.794733 (1 , 2) 4.794733 ω n/a 5.022646 5.000028 (5 , 0) 5.000000 ω , 3) 5.692100 ω n/a 5.797912 5.794737 (2 , 2) 5.794733 ω n/a 6.142955 6.000676 (6 , 0) 6.000000 Table C.1: Benchmark of LR-MCTDHB to the attractive 1D HIM. Shown are the ground-state energy ε and the energies ω i of the first ten excitations for N = 10 bosons and different numbers of orbitals M . The trapping frequency is Ω = 1 . λ = 0 . 13, yielding a ground-state depletion of f = 0 . = 1 M = 2 M = 3 ( m x , m y n x , n y ) Exact ε , , , 0) 110.003452 ω , , , 0) 1.000000 ω , , , 0) 1.378405 ω , , , 0) 1.788854 ω n/a 2.000438 2.000477 (2 , , , 0) 2.000000 ω , , , 1) 2.198268 ω n/a 1.940430 2.378704 (1 , , , 0) 2.378405 ω , , , 2) 2.607681 ω , , , 0) 2.683282 ω n/a n/a 2.757078 (0 , , , 0) 2.756810 ω n/a 2.789191 2.789227 (1 , , , 0) 2.788854 ω n/a 3.000758 3.000816 (3 , , , 0) 3.000000 Table C.2: Benchmark of LR-MCTDHB to the anisotropic 2D HIM with N = 100 bosons. The trapping frequencies are Ω x = 1 . y = √ . x - and y -directions, and the two-body interaction strength id λ = − . M = 1 M = 3 ( n, m, l z ) Exact ε ω ω ω ω ω n/a 1.800439 (0,2,2) 1.800000 ω ω n/a 2.000657 (0,2,0) 2.000000 ω n/a 2.200438 (0,2,-2) 2.200000 ω ω n/a 2.487088 (2,1,3) 2.488854 ω ω n/a 2.689270 (2,1,1) 2.688854 ω n/a 2.691535 (2,1,1) 2.688854 ω n/a 2.700528 (0,3,3) 2.700000 ω Table C.3: Benchmark of LR-MCTDHB to the isotropic 2D HIM in the rotating frame. The angular velocity around the z -axisis Ω rot = 0 . 1. Results are presented for N = 100 bosons and different numbers of orbitals M . The trapping frequencies areΩ x = Ω y = 1 . 0, whereas the strength of the repulsive interaction is λ = − . ε , and of the first 15 excited states, ω i = E i − ε . Underlined digits denote deviations from the exact values from Eqs. (160),(161) and (C.1). The quantum numbers n and m refer to the relative and c.m. coordinates, whereas l z refers to the angularmomentum in the z -direction. All quantities are dimensionless. See text for more details. Finally, excitations in the rotating frame of reference are discussed. The energies of the ground and excitedstates, given by Eqs. (160) and (161), become shifted by the term E rot = − Ω rot l z (C.1)where Ω rot denotes the angular velocity around the z -axis and l z denotes the angular momentum of the89xited state in the z -direction. In Table 3, all obtained excitations are labeled by corresponding quantumnumbers for relative and c.m. excitations, n and m , as well as by l z . The increased numerical accuracy ofLR-MCTDHB(3) in comparison to the BdG results becomes clearly obvious. Moreover, all states are obtainedfor M = 3, whereas the BdG spectrum misses several states. As a conclusion, the utilized implementation ofLR-MCTDHB is also capable of accurately describing the low-energy excitation spectra in the rotating frame. Appendix D Many-particle variance: A sensitive quantity to many-body effects In this Appendix, a brief introduction of the many-body variance, which is employed in the main text tostudy many-body excitations involved in the dynamics of trapped BECs, is given. The many-body variancehas been of interest in recent studies on correlations between bosons in the ground state [242] and for out-of-equilibrium BECs [243, 244], in particular in a 1D bosonic Josephson junction [245] or in anharmonic andanisotropic trapping potentials [246, 81]. It has also been found that the variance is a sensitive measure forthe numerical convergence of MCTDHB [80]. The focus is laid on the position variance in the x -direction inthe following.The position operator of N bosons is given byˆ X = N (cid:88) i =1 ˆ x i (D.1)where the operator ˆ x i denotes the position operator of the i -th particle. To compute its variance per particlein the state Ψ( r , ..., r N ), which reads1 N ∆ X = 1 N (cid:16) (cid:104) Ψ | ˆ X | Ψ (cid:105) − (cid:104) Ψ | ˆ X | Ψ (cid:105) (cid:17) , (D.2)one also needs to evaluate the expectation value of the square of the position operator. The latter is given byˆ X = N (cid:88) i =1 ˆ x i + (cid:88) i Bibliography [1] S. N. Bose, Plancks Gesetz und Lichtquantenhypothese. , Z. Phys. , 178 (1924).[2] A. Einstein, Quantentheorie des einatomigen idealen Gases. , Sitzungsber. Preuss. Akad. Wiss., Bericht , 261 (1924).[3] A. Einstein, Quantentheorie des einatomigen idealen Gases. Zweite Abhandlung. , Sitzungsber. Preuss.Akad. Wiss., Bericht , 3 (1925).[4] M. H. Anderson, J. R. Ensher, M. R. Matthews, C. E. Wieman, and E. Cornell, Observation of Bose-Einstein condensation in a dilute atomic vapor , Science , 198 (1995).[5] K. B. Davis, M.-O. Mewes, M. R. Andrews, N. J. van Druten, D. S. Durfee, D. M. Kurn, and W. Ketterle, Bose-Einstein condensation in a gas of sodium atoms , Phys. Rev. Lett. , 3969 (1995).[6] C. C. Bradley, C. A. Sackett, J. J. Tollett, and R. G. Hulet, Evidence of Bose-Einstein Condensation inan Atomic Gas with Attractive Interactions , Phys. Rev. Lett. , 1687 (1995).[7] J. Kasprzak, M. Richard, S. Kundermann, A. Baas, P. Jeambrun, J. M. J. Keeling, F. M. Marchetti,M. H. Szyma´nska, R. Andr´e, J. L. Staehli, V. Savona, P. B. Littlewood, B. Deveaud, and L. S. Dang, Bose-Einstein condensation of exciton polaritons , Nature , 409 (2006).[8] H. Deng, H. Haug, and Y. Yamamoto, Exciton-polariton Bose-Einstein condensation , Rev. Mod. Phys. , 1489 (2010).[9] S. O. Demokritov, V. E. Demidov, O. Dzyapko, G. A. Melkov, A. A. Serga, B. Hillebrands, and A. N.Slavin, Bose-Einstein condensation of quasi-equilibrium magnons at room temperature under pumping ,Nature volume , 430 (2006).[10] J. Klaers, J. Schmitt, F. Vewinger, and M. Weitz, Bose-Einstein condensation of photons in an opticalmicrocavity , Nature , 545 (2010).[11] P. Courteille, R. S. Freeland, D. J. Heinzen, F. A. van Abeelen, and B. J. Verhaar, Observation of aFeshbach Resonance in Cold Atom Scattering , Phys. Rev. Lett. , 69 (1998).[12] S. Inouye, M. R. Andrews, J. Stenger, H.-J. Miesner, D. M. Stamper-Kurn, and W. Ketterle, Observationof Feshbach resonances in a Bose-Einstein condensate , Nature , 151 (1998).[13] C. Chin , R. Grimm, P. Julienne , and E. Tiesinga, Feshbach resonances in ultracold gases , Rev. Mod.Phys. , 1225 (2010).[14] A. G¨orlitz, J. M. Vogels, A. E. Leanhardt, C. Raman, T. L. Gustavson, J. R. Abo-Shaeer, A. P. Chikkatur,S. Gupta, S. Inouye, T. Rosenband, and W. Ketterle, Realization of Bose-Einstein Condensates in LowerDimensions , Phys. Rev. Lett. , 130402 (2001).9315] O. Lahav, A. Itah, A. Blumkin, C. Gordon, S. Rinott, A. Zayats, and J. Steinhauer, Realization of aSonic Black Hole Analog in a Bose-Einstein Condensate , Phys. Rev. Lett. , 240401 (2010).[16] J. Steinhauer, Observation of self-amplifying Hawking radiation in an analogue black-hole laser , Nat.Phys. , 864 (2014).[17] J. Steinhauer, Observation of quantum Hawking radiation and its entanglement in an analogue black hole ,Nat. Phys. , 959 (2016).[18] E. P. Gross, Structure of a quantized vortex in boson systems , Nuovo Cimento , 454 (1961).[19] L. P. Pitaevskii, Vortex lines in an imperfect Bose gas , Zh. Eksp. Teor. Fiz. , 646 (1961) [Sov. Phys.JETP , 451 (1961)].[20] L. P. Pitaevskii and S. Stringari, Bose-Einstein Condensation , Oxford University Press, New York, 2003.[21] C. J. Pethick and H. Smith, Bose-Einstein Condensation in Dilute Gases , 2nd edition, Cambridge Uni-versity Press, New York, 2008.[22] A. J. Leggett, Bose-Einstein condensation in the alkali gases: Some fundamental concepts , Rev. Mod.Phys. , 307 (2003).[23] N. N. Bogoliubov, On the theory of superfluidity , J. Phys. (USSR), 23 (1947).[24] P.-G. de Gennes, Superconductivity of metals and alloys , Benjamin, New York, 1966.[25] G. C. Katsimiga, P. G. Kevrekidis, B. Prinari, G. Biondini, and P. Schmelcher, Dark-bright soliton pairs:Bifurcations and collisions , Phys. Rev A , 043623 (2018).[26] L. Jia, A.-B. Wang, and S. Yi, Low-lying excitations of vortex lattices in condensates with anisotropicdipole-dipole interaction , Phys. Rev. A , 043614 (2018).[27] B. Galilo, D. K. K. Lee, and R. Barnett, Topological Edge-State Manifestation of Interacting 2D Con-densed Boson-Lattice Systems in a Harmonic Trap , Phys. Rev. Lett. , 203204 (2017).[28] J.-X. Zhu, Bogoliubov-de Gennes Method and Its Applications , Springer, Heidelberg, 2016.[29] A. Spuntarelli, P. Pieri, G. C. Strinati, Solution of the Bogoliubov-de Gennes equations at zero temperaturethroughout the BCS-BEC crossover: Josephson and related effects , Phys. Rep. , 111 (2010).[30] F. Dalfovo, S. Giorgini, L. P. Pitaevskii, and S. Stringari, Theory of Bose-Einstein condensation in trappedgases , Rev. Mod. Phys. , 463 (1999).[31] O. E. Alon, A. I. Streltsov, and L. S. Cederbaum, Unified view on multiconfigurational time-propagationfor systems consisting of identical particles , J. Chem. Phys. , 154103 (2007).[32] O. E. Alon, A. I. Streltsov, and L. S. Cederbaum, Multiconfigurational time-dependent Hartree methodfor bosons: Many-body dynamics of bosonic systems , Phys. Rev. A , 033613 (2008).[33] A. U. J. Lode, K. Sakmann, O. E. Alon, L. S. Cederbaum, and A. I. Streltsov, Numerically exact quantumdynamics of bosons with time-dependent interactions of harmonic type , Phys. Rev. A , 063606 (2012).[34] J. Grond, A. I. Streltsov, A. U. J. Lode, K. Sakmann, L. S. Cederbaum, and O. E. Alon, Excitationspectra of many-body systems by linear response: General theory and applications to trapped condensates ,Phys. Rev. A , 023606 (2013). 9435] O. E. Alon, A. I. Streltsov, and L. S. Cederbaum, Unified view on linear-response of interacting identicaland distinguishable particles from multiconfigurational time-dependent Hartree methods , J. Chem. Phys , 034108 (2014).[36] Markus Fierz, ¨Uber die relativistische Theorie kr¨aftefreier Teilchen mit beliebigem Spin , Helvetica PhysicaActa , 3 (1939).[37] Wolfgang Pauli, The Connection Between Spin and Statistics , Phys. Rev. , 716 (1940).[38] P. Kramer and M. Sarenco, Geometry of the Time-Dependent Variational Principle , Springer, Heidelberg,1981.[39] P. A. M. Dirac, Note on Exchange Phenomena in the Thomas Atom , Proc. Cambridge Philos. Soc. ,376 (1930).[40] J. Frenkel, Wave Mechanics: advanced general theory , Oxford, Clarendon Press (1934).[41] A. D. McLachlan, A variational solution of the time-dependent Schr¨odinger equation , Mol. Phys. , 39(1964).[42] J. Kucar, H.-D. Meyer, and L. S. Cederbaum, Time-dependent rotated Hartree approach , Chem. Phys.Lett. , 525 (1987).[43] O. Penrose and L. Onsager, Bose-Einstein Condensation and Liquid Helium , Phys. Rev. , 576 (1956).[44] P. Nozi`eres, D. Saint James, Particle vs. pair condensation in attractive Bose liquids , J. Phys. (France) , 1133 (1982).[45] A. J. Leggett, Quantum Liquids: Bose condensation and Cooper pairing in condensed matter systems ,Oxford University Press, Oxford (2006).[46] K. Sakmann, Numerically exact dynamics of the interacting many-body Schr¨odinger equation for Bose-Einstein condensates: comparison to Bose-Hubbard and Gross-Pitaevskii theory , PhD thesis, HeidelbergUniversity, 2010.[47] H.-D. Meyer, U. Manthe, and L. S. Cederbaum, The Multi-Configurational Time-Dependent HartreeApproach , Chem. Phys. Lett. , 73 (1990).[48] U. Manthe, H.-D. Meyer, and L. S. Cederbaum, Wave-packet dynamics within the multiconfigurationHartree framework: General aspects and application to NOCl , J. Chem. Phys. , 3199 (1992).[49] H.–D. Meyer, F. Gatti, and G. A. Worth, Multidimensional Quantum Dynamics: MCTDH Theory andApplications , Wiley-VCH, 2009.[50] O. E. Alon, A. I. Streltsov, and L. S. Cederbaum, Many-body theory for systems with particle conversion:Extending the multiconfigurational time-dependent Hartree method , Phys. Rev. A , 022503 (2009).[51] O. E. Alon, A. I. Streltsov, and L. S. Cederbaum, Multiconfigurational time-dependent Hartree methodfor mixtures consisting of two types of identical particles , Phys. Rev. A. , 062501 (2007).[52] A. U. J. Lode and C. Bruder, Dynamics of Hubbard Hamiltonians with the multiconfigurational time-dependent Hartree method for indistinguishable particles , Phys. Rev. A , 013616 (2016).[53] A. U. J. Lode, Multiconfigurational time-dependent Hartree method for bosons with internal degrees offreedom: Theory and composite fragmentation of multicomponent Bose-Einstein condensates , Phys. Rev.A , 063601 (2016). 9554] S. R. White, Density matrix formulation for quantum renormalization groups , Phys. Rev. Lett. , 2363(1992).[55] K. A. Hallberg, New trends in density matrix renormalization , Adv. Phys. , 477 (2006).[56] E. M. Stoudenmire and S. R. White, Studying Two-Dimensional Systems with the Density Matrix Renor-malization Group , Annu. Rev. Cond. Matt. Phys. , 111 (2012).[57] A. I. Streltsov, O. E. Alon, L. S. Cederbaum, General variational many-body theory with complete self-consistency for trapped bosonic systems , Phys. Rev. A , 063626 (2006).[58] Y. Saad, Iterative Methods for Sparse Linear Systems ∼ Numerical Methods for Eigenvalue Problems , De Gruyter,Berlin/Boston, 2012.[64] Y. Saad, Numerical Methods for Large Eigenvalue Problems , Halstead Press, New York, 1992.[65] A. I. Streltsov, L. S. Cederbaum, O. E. Alon, K. Sakmann, A. U. J. Lode, J. Grond, O. I. Streltsova, S.Klaiman, and R. Beinke, The Multiconfigurational Time-Dependent Hartree for Bosons Package , version3.x, http://mctdhb.org, Heidelberg/Kassel (2006-present).[66] A. I. Streltsov, O. E. Alon, and L. S. Cederbaum, General mapping for bosonic and fermionic operatorsin Fock space , Phys. Rev. A 81, 022124 (2010).[67] A. I. Streltsov and O. I. Streltsova, MCTDHB-Lab Non-equilibrium quantum dynamics of ultra-coldbosonic mixtures: the multi-layermulti-configuration time-dependent Hartree method for bosons , New J.Phys. , 063018 (2013).[70] L. Cao, S. Kr¨onke, O. Vendrell, and P. Schmelcher, The multi-layer multi-configuration time-dependentHartree method for bosons: Theory, implementation, and applications , J. Chem. Phys. , 134103(2013).[71] L. Cao, V. Bolsinger, S. I. Mistakidis, G. M. Koutentakis, S. Kr¨onke, J. M. Schurer, and P. Schmelcher, Aunified ab initio approach to the correlated quantum dynamics of ultracold fermionic and bosonic mixtures ,J. Chem. Phys. , 044106 (2017).[72] L. Cohen and C. Lee, Exact reduced density matrices for a model problem , J. Math. Phys. , 3105(1985).[73] J. Yan, Harmonic Interaction Model and Its Applications in Bose-Einstein Condensation , J. Stat. Phys. , 314 (2003).[74] A. I. Streltsov, O. E. Alon, and L. S. Cederbaum, Role of Excited States in the Splitting of a TrappedInteracting Bose-Einstein Condensate by a Time-Dependent Barrier , Phys. Rev. Lett. , 030402 (2007).9675] K. Sakmann, A. I. Streltsov, O. E. Alon, and L. S. Cederbaum, Reduced density matrices and coherenceof trapped interacting bosons , Phys. Rev. A , 023615 (2008).[76] S. Klaiman, A. U. J. Lode, A. I. Streltsov, L. S. Cederbaum, and O. E. Alon, Breaking the resilience ofa two-dimensional Bose-Einstein condensate to fragmentation , Phys. Rev. A , 043620 (2014).[77] O. E. Alon, A. I. Streltsov, and L. S. Cederbaum, Build-up of coherence between initially-independentsubsystems: The case of Bose-Einstein condensates , Phys. Lett. A , 9 (2009).[78] A. I. Streltsov, O. E. Alon, and L. S. Cederbaum, Efficient generation and properties of mesoscopicquantum superposition states in an attractive Bose-Einstein condensate threaded by a potential barrier ,J. Phys. B: At. Mol. Opt. Phys. , 091004 (2009).[79] A. I. Streltsov, O. E. Alon, and L. S. Cederbaum, Formation and Dynamics of Many-Boson FragmentedStates in One-Dimensional Attractive Ultracold Gases , Phys. Rev. Lett. , 130401 (2008).[80] J. G. Cosme, C. Weiss, and J. Brand, Center-of-mass motion as a sensitive convergence test for variationalmultimode quantum dynamics , Phys. Rev. A , 043603 (2016).[81] O. E. Alon and L. S. Cederbaum, Attractive Bose-Einstein condensates in anharmonic traps: Accuratenumerical treatment and the intriguing physics of the variance , Chem. Phys. , 287 (2018).[82] C. Weiss and Y. Castin, Creation and Detection of a Mesoscopic Gas in a Nonlocal Quantum Superpo-sition , Phys Rev. Lett. , 010403 (2009).[83] A. I. Streltsov, O. E. Alon, and L. S. Cederbaum, Scattering of an attractive Bose-Einstein condensatefrom a barrier: Formation of quantum superposition states , Phys. Rev. A , 043616 (2009).[84] A. I. Streltsov, O. E. Alon, and L. S. Cederbaum, Swift Loss of Coherence of Soliton Trains in AttractiveBose-Einstein Condensates , Phys. Rev. Lett. , 240401 (2011).[85] J. Grond, J. Schmiedmayer, and U. Hohenester, Optimizing number squeezing when splitting a mesoscopiccondensate , Phys. Rev. A , 021603(R) (2009).[86] J. Grond, G. von Winckel, J. Schmiedmayer, and U. Hohenester, Optimal control of number squeezing intrapped Bose-Einstein condensates , Phys. Rev. A , 053625 (2009).[87] J. Grond, J. Schmiedmayer, and U. Hohenester, Shaking the condensates: Optimal number squeezing inthe dynamic splitting of a Bose-Einstein condensate , Physica E , 432 (2010).[88] G. J¨ager, T. Berrada, J. Schmiedmayer, T. Schumm, and U. Hohenester, Parametric-squeezing amplifi-cation of Bose-Einstein condensates , Phys. Rev. A , 053632 (2015).[89] J. Grond, U. Hohenester, I. Mazets, and J. Schmiedmayer, Atom interferometry with trapped Bose-Einstein condensates: impact of atom-atom interactions , New J. Phys. , 065036 (2010).[90] J. Grond, T. Betz, U. Hohenester, N. J. Mauser, J. Schmiedmayer, and T. Schumm, The Shapiro effectin atomchip-based bosonic Josephson junctions , New J. Phys. Quantumspeed limit and optimal control of many-boson dynamics , Phys. Rev. A , 062110 (2015).[92] J. Grond, U. Hohenester, J. Schmiedmayer, and A. Smerzi, Mach-Zehnder interferometry with interactingtrapped Bose-Einstein condensates , Phys. Rev. A , 023619 (2011).[93] I. Bˇrezinov´a, A. U. J. Lode, A. I. Streltsov, O. E. Alon, L. S. Cederbaum, and J. Burgd¨orfer, Wave chaosas signature for depletion of a Bose-Einstein condensate , Phys. Rev. A , 013630 (2012).9794] I. Bˇrezinov´a, J. Burgd¨orfer, A. U. J. Lode, A. I. Streltsov, L. S. Cederbaum, O. E. Alon, L. A. Collins,and B. I. Schneider, Elastic scattering of a Bose-Einstein condensate at a potential landscape , J. Phys.:Conf. Ser. , 012032 (2014).[95] K. Sakmann and M. Kasevich, Single-shot simulations of dynamic quantum many-body systems , Nat.Phys. , 451 (2016).[96] R. Schmitz, S. Kr¨onke, L. Cao, and P. Schmelcher, Quantum breathing dynamics of ultracold bosons inone-dimensional harmonic traps: Unraveling the pathway from few- to many-body systems , Phys. Rev. A , 043601 (2013).[97] S. Kr¨onke and P. Schmelcher, Many-body processes in black and gray matter-wave solitons , Phys. Rev. A , 053614 (2015).[98] G. C. Katsimiga, S. I. Mistakidis, G. M. Koutentakis, P. G. Kevrekidis, and P. Schmelcher, Many-bodyquantum dynamics in the decay of bent dark solitons of Bose-Einstein condensates , New J. Phys. ,123012 (2017).[99] S. Kr¨onke, J. Kn¨orzer, and P. Schmelcher, Correlated quantum dynamics of a single atom collisionallycoupled to an ultracold finite bosonic ensemble , New J. Phys. , 053001 (2015).[100] J. M. Schurer, A. Negretti, and P. Schmelcher, Unraveling the Structure of Ultracold MesoscopicCollinear Molecular Ions , Phys. Rev. Lett. , 063001 (2017).[101] Q.-S. Tan, H.-Y. Lu, and S. Yi, Spin squeezing of a dipolar Bose gas in a double-well potential , Phys.Rev. A , 013606 (2016).[102] B. Chatterjee and A. U. J. Lode, Order parameter and detection for a finite ensemble of crystallizedone-dimensional dipolar bosons in optical lattices , Phys. Rev A , 053624 (2018).[103] B. Chatterjee, M. C. Tsatsos, and A. U. J. Lode, Correlations of strongly interacting one-dimensionalultracold dipolar few-boson system in optical lattices , arXiv:1806.05048 [cond-mat.quant-gas] (2018), ac-cepted at New J. Phys.[104] A. U. J. Lode and C. Bruder, Fragmented Superradiance of a Bose-Einstein Condensate in an OpticalCavity , Phys. Rev. Lett. , 013603 (2017).[105] P. Molignini, L. Papariello, A. U. J. Lode, and R. Chitra, Superlattice switching from parametric insta-bilities in a driven-dissipative Bose-Einstein condensate in a cavity , Phys. Rev A , 053620 (2018).[106] A. U. J. Lode, F. S. Diorico, R. Wu, P. Molignini, L. Papariello, R. Lin, C. L´evˆeque, L. Exl, M. C.Tsatsos, R. Chitra, and N. J. Mauser, Many-body physics in two-component Bose-Einstein condensatesin a cavity: fragmented superradiance and polarization , New J. Phys. , 055006 (2018).[107] M. Albiez, R. Gati, J. F¨olling, S. Hunsmann, M. Cristiani, and M. K. Oberthaler, Direct Observationof Tunneling and Nonlinear Self-Trapping in a Single Bosonic Josephson Junction , Phys. Rev. Lett. ,010402 (2005).[108] S. Levy, E. Lahoud, I. Shomroni, and J. Steinhauer, The a.c. and d.c. Josephson effects in a Bose-Einstein condensate , Nature (London) , 579 (2007).[109] G. J. Milburn J. Corney, E. M. Wright, and D. F. Walls, Quantum dynamics of an atomic Bose-Einsteincondensate in a double-well potential , Phys. Rev. A , 4318 (1997).[110] A. Smerzi, S. Fantoni, S. Giovanazzi, and S. R. Shenoy, Quantum Coherent Atomic Tunneling betweenTwo Trapped Bose-Einstein Condensates , Phys. Rev. Lett. Coherent oscillations between two weakly coupledBose-Einstein condensates: Josephson effects, π oscillations, and macroscopic quantum self-trapping ,Phys. Rev. A , 620 (1999).[112] E. A. Ostrovskaya, Y. S. Kivshar, M. Lisak, B. Hall, F. Cattani, and D. Anderson, Coupled-mode theoryfor Bose-Einstein condensates , Phys. Rev. A , 031601(R).[113] D. Ananikian and T. Bergeman, Gross-Pitaevskii equation for Bose particles in a double-well potential:Two-mode models and beyond , Phys. Rev. A , 013604 (2006).[114] X. Y. Jia, W. D. Li, and J. Q. Liang, Nonlinear correction to the boson Josephson-junction model , Phys.Rev. A , 023613 (2008).[115] Y. Zhou, H. Zhai, R. L¨u, Z. Xu, and L. Chang, Quantum coherence of double-well Bose-Einsteincondensates: An SU(2) coherent-state path-integral approach , Phys. Rev. A , 043606 (2003).[116] G. Ferrini, A. Minguzzi, and F. W. J. Hekking, Number squeezing, quantum fluctuations, and oscillationsin mesoscopic Bose Josephson junctions , Phys. Rev. A , 023606 (2008).[117] V. S. Shchesnovich and M. Trippenbach, Fock-space WKB method for the boson Josephson model de-scribing a Bose-Einstein condensate trapped in a double-well potential , Phys. Rev. A A bosonic Josephson junction , J. Phys. B , R61 (2007).[119] J. M. Schurer, A. Negretti, and P. Schmelcher, Capture dynamics of ultracold atoms in the presence ofan impurity ion , New J. Phys. Ground-state properties of ultracold trapped bosons withan immersed ionic impurity , Phys. Rev. A , 033601 (2014).[121] J. M. Schurer, R. Gerritsma, P. Schmelcher, and A. Negretti, Impact of many-body correlations on thedynamics of an ion-controlled bosonic Josephson junction , Phys. Rev. A , 063602 (2016).[122] M. Heimsoth, D. Hochstuhl, C. E. Creffield, L. D. Carr, and F. Sols, Effective Josephson dynamics inresonantly driven Bose-Einstein condensates , New J. Phys. , 103006 (2013).[123] K. Sakmann, A. I. Streltsov, O. E. Alon, and L. S. Cederbaum, Exact Quantum Dynamics of a BosonicJosephson Junction , Phys. Rev. Lett. , 220601 (2009).[124] K. Sakmann, A. I. Streltsov, O. E. Alon, and L. S. Cederbaum, Quantum dynamics of attractive versusrepulsive bosonic Josephson junctions: Bose-Hubbard and full-Hamiltonian results , Phys. Rev. A ,013620 (2010).[125] K. Sakmann, Many-Body Schr¨odinger Dynamics of Bose-Einstein Condensates , Springer, Heidelberg,2011.[126] K. Sakmann, A. I. Streltsov, O. E. Alon, and L. S. Cederbaum, Universality of fragmentation in theSchr¨odinger dynamics of bosonic Josephson junctions , Phys. Rev. A , 023602 (2014).[127] J. G. Cosme, M. F. Andersen, and J. Brand, Interaction blockade for bosons in an asymmetric doublewell , Phys. Rev. A , 013616 (2017).[128] R. Beinke, S. Klaiman, L. S. Cederbaum, A. I. Streltsov, and O. E. Alon, Many-body tunneling dynamicsof Bose-Einstein condensates and vortex states in two spatial dimensions , Phys. Rev. A , 043627 (2015).[129] D. P. Arovas and A. Auerbach, Quantum tunneling of vortices in two-dimensional superfluids , Phys.Rev. B , 094508 (2008). 99130] A. M. Martin, R. G. Scott, and T. M. Fromhold, Transmission and reflection of Bose-Einstein conden-sates incident on a Gaussian tunnel barrier , Phys. Rev. A , 065602 (2007).[131] O. Fialko, A. S. Bradley, and J. Brand, Quantum tunneling of a vortex between two pinning potentials ,Phys. Rev. Lett. , 015301 (2012).[132] P. Ao and D. J. Thouless, Tunneling of a quantized vortex: Roles of pinning and dissipation , Phys. Rev.Lett. , 132 (1994).[133] A. Auerbach, D. P. Arovas, and S. Ghosh, Quantum tunneling of vortices in two-dimensional conden-sates , Phys. Rev. B , 064511 (2006).[134] J. R. Salgueiro, M. Zacar´es, H. Michinel, and A. Ferrando, Vortex replication in Bose-Einstein conden-sates trapped in double-well potentials , Phys. Rev. A , 033625 (2009).[135] S. Klaiman and O. E. Alon, Spatially-partitioned many-body vortices , J. Phys.: Conf. Ser. , 012015(2016).[136] J. Christensson, C. Forss´en, S. ˚Aberg, and S. M. Reimann, Effective-interaction approach to the many-boson problem , Phys. Rev. A , 012707 (2009).[137] R. A. Doganov, S. Klaiman, O. E. Alon, A. I. Streltsov, and L. S. Cederbaum, Two trapped particlesinteracting by a finite-range two-body potential in two spatial dimensions , Phys. Rev. A , 033631 (2013).[138] U. R. Fischer, A. U. J. Lode, and B. Chatterjee, Condensate fragmentation as a sensitive measure ofthe quantum many-body behavior of bosons with long-range interactions , Phys. Rev. A , 063621 (2015).[139] K. Sakmann and J. Schmiedmayer, Conserving symmetries in Bose-Einstein condensate dynamics re-quires many-body theory , arXiv:1802.03746v2 [cond-mat.quant-gas] (2018).[140] M. Razavy, Quantum Theory of Tunneling , World Scientific, Singapore, 2003.[141] H. A. Kramers, Wave mechanics and half-integral quantization , Z. Phys. A , 828 (1926).[142] R. W. Gurney and E. U. Condon, Wave mechanics and radioactive disintegration , Naure (London) ,439 (1928).[143] R. W. Gurney and E. U. Condon, Quantum mechanics and radioactive disintegration , Phys. Rev. ,127 (1929).[144] N. Teeny, E. Yakaboylu, H. Bauke, and C. H. Keitel, Ionization Time and Exit Momentum in Strong-Field Tunnel Ionization , Phys. Rev. Lett. , 063003 (2016).[145] N. Teeny, C. H. Keitel, and H. Bauke, Virtual-detector approach to tunnel ionization and tunnelingtimes , Phys. Rev. A , 022104 (2016).[146] L. D. Carr, M. J. Holland, and B. A. Malomed, Macroscopic quantum tunneling of Bose-Einsteincondensates in a finite potential well , J. Phys. B , 3217 (2005).[147] N. Moiseyev and L. S. Cederbaum, Resonance solutions of the nonlinear Schr¨odinger equation: Tunnel-ing lifetime and fragmentation of trapped condensates , Phys. Rev. A , 033605 (2005).[148] S. Kim and J. Brand, Decay modes of two repulsively interacting bosons , J. Phys. B , 195301 (2011).[149] S. Hunn, K. Zimmermann, M. Hiller, and A. Buchleitner, Tunneling decay of two interacting bosons inan asymmetric double-well potential: A spectral approach , Phys. Rev. A , 043626 (2013).100150] A. U. J. Lode, A. I. Streltsov, O. E. Alon, H.-D. Meyer, and L. S. Cederbaum, Exact decay and tunnellingdynamics of interacting few-boson systems , J. Phys. B: At. Mol. Opt. Phys. , 044018 (2009).[151] A. U. J. Lode, A. I. Streltsov, O. E. Alon, H.-D. Meyer, and L. S. Cederbaum, Corrigendum: Exactdecay and tunnelling dynamics of interacting few-boson systems , J. Phys. B: At. Mol. Opt. Phys. ,029802 (2010).[152] Axel U. J. Lode, A. I. Streltsov, K. Sakmann, O. E. Alon, and L. S. Cederbaum, How an interactingmany-body system tunnels through a potential barrier to open space , Proc. Natl. Acad. Sci. USA ,13521 (2012).[153] Axel U. J. Lode, S. Klaiman, O. E. Alon, A. I. Streltsov, L. S. Cederbaum, Controlling the velocitiesand the number of emitted particles in the tunneling to open space dynamics , Phys. Rev. A , 053620(2014).[154] Axel U. J. Lode, Tunneling Dynamics in Open Ultracold Bosonic Systems , Springer, Heidelberg, 2015.[155] K. Kasamatsu, M. Tsubota, and M. Ueda, Vortex Phase Diagram in Rotating Two-Component Bose-Einstein Condensates , Phys. Rev. Lett. , 150406 (2003).[156] J. Lovegrove, M. O. Borgh, and J. Ruostekoski, Energetic Stability of Coreless Vortices in Spin-1 Bose-Einstein Condensates with Conserved Magnetization , Phys. Rev. Lett. , 075301 (2014).[157] E. A. Cornell and C. E. Wieman, Nobel Lecture: Bose-Einstein condensation in a dilute gas, the first70 years and some recent experiments , Rev. Mod. Phys. , 875 (2002).[158] W. Ketterle, Nobel lecture: When atoms behave as waves: Bose-Einstein condensation and the atomlaser , Rev. Mod. Phys. , 1131 (2002).[159] J. R. Abo-Shaeer, C. Rama, J. M. Vogels, and W. Ketterle, Observation of vortex lattices in Bose-Einstein condensates , Science , 476 (2001).[160] D. Dagnino, N. Barber´an, M. Lewenstein, and J. Dalibard, Vortex nucleation as a case study of symmetrybreaking in quantum systems , Nat. Phys. , 431 (2009).[161] D. A. Butts and D. S. Rokhsar, Predicted signatures of rotating Bose-Einstein condensates , Nature ,327 (1999).[162] K. W. Madison, F. Chevy, W. Wohlleben, and J. Dalibard, Vortex Formation in a Stirred Bose-EinsteinCondensate , Phys. Rev. Lett. , 806 (2000).[163] F. Chevy, K. W. Madison, and J. Dalibard, Measurement of the Angular Momentum of a RotatingBose-Einstein Condensate , Phys. Rev. Lett. , 2223 (2000).[164] K. W. Madison, F. Chevy, V. Bretin, and J. Dalibard, Stationary States of a Rotating Bose-EinsteinCondensate: Routes to Vortex Nucleation , Phys. Rev. Lett. , 4443 (2001).[165] S. E. Weiner, M. C. Tsatsos, L. S. Cederbaum, and A. U. J. Lode, Phantom vortices: hidden angularmomentum in ultracold dilute Bose-Einstein condensates , Sci. Rep. , 40122 (2017). License for use ofcontent: https://creativecommons.org/licenses/by/4.0/legalcode.[166] M. C. Tsatsos and A. U. J. Lode, Resonances and Dynamical Fragmentation in a Stirred Bose-EinsteinCondensate , J. Low Temp. Phys. , 171 (2015).[167] M. Theisen and A. I. Streltsov, Many-body excitations and deexcitations in trapped ultracold bosonicclouds , Phys. Rev. A , 053622 (2016). 101168] R. Beinke, S. Klaiman, L. S. Cederbaum, A. I. Streltsov, and O. E. Alon, Many-body effects in theexcitation spectrum of weakly interacting Bose-Einstein condensates in one-dimensional optical lattices ,Phys. Rev. A , 063602 (2017).[169] S. Burger, F. S. Cataliotti, C. Fort, F. Minardi, M. Inguscio, M. L. Chiofalo, and M. P. Tosi, Superfluidand Dissipative Dynamics of a Bose-Einstein Condensate in a Periodic Optical Potential , Phys. Rev.Lett. , 4447 (2001).[170] I. Bloch, Ultracold quantum gases in optical lattices , Nature Physics , 23 (2005).[171] K. Suthar, A. Roy, and D. Angom, Fluctuation-driven topological transition of binary condensates inoptical lattices , Phys. Rev. A , 043615 (2015).[172] O. Morsch and M. Oberthaler, Dynamics of Bose-Einstein condensates in optical lattices , Rev. Mod.Phys. , 179 (2006).[173] A. K. Tuchman, C. Orzel, A. Polkovnikov, and M. A. Kasevich, Nonequilibrium coherence dynamics ofa soft boson lattice , Phys. Rev. A , 051601(R) (2006).[174] A. K. Tuchman, W. Li, S. Dettmer, and M. A. Kasevich, Localization and anomalous transport in a 1Dsoft boson optical lattice , New J. Phys. , 311 (2006).[175] K. Xu, Y. Liu, D. E. Miller, J. K. Chin, W. Setiawan, and W. Ketterle, Observation of Strong QuantumDepletion in a Gaseous Bose-Einstein Condensate , Phys. Rev. Lett. , 180405 (2006).[176] S. I. Mistakidis, L. Cao, and P. Schmelcher, Negative-quench-induced excitation dynamics for ultracoldbosons in one-dimensional lattices , Phys. Rev. A , 033611 (2015).[177] G. M. Koutentakis, S. I. Mistakidis, and P. Schmelcher, Quench-induced resonant tunneling mechanismsof bosons in an optical lattice with harmonic confinement , Phys. Rev. A , 013617 (2017).[178] J. Neuhaus-Steinmetz, S. I. Mistakidis, and P. Schmelcher, Quantum dynamical response of ultracoldfew-boson ensembles in finite optical lattices to multiple interaction quenches , Phys. Rev. A , 053610(2017).[179] S. I. Mistakidis, T. Wulf, A. Negretti, and P. Schmelcher, Resonant quantum dynamics of few ultracoldbosons in periodically driven finite lattices , J. Phys. B: At. Mol. Opt. Phys. , 244004 (2015).[180] S. I. Mistakidis and P. Schmelcher, Mode coupling of interaction quenched ultracold few-boson ensemblesin periodically driven lattices , Phys. Rev. A , 013625 (2017).[181] A. I. Streltsov, K. Sakmann, O. E. Alon, and L. S. Cederbaum, Accurate multi-boson long-time dynamicsin triple-well periodic traps , Phys. Rev. A , 043604 (2011).[182] R. Roy, A. Gammal, M. C. Tsatsos, B. Chatterjee, B. Chakrabarti, and A. U. J. Lode, Phases, many-body entropy measures, and coherence of interacting bosons in optical lattices , Phys. Rev. A , 043625(2018).[183] S. Dutta, M. C. Tsatsos, S. Basu, and A. U. J. Lode, Management of the Correlations of UltracoldBosons in Triple Wells , arXiv:1802.02407 [cond-mat.quant-gas] (2018).[184] L. Cao, S. Kr¨onke, J. Stockhofe, J. Simonet, K. Sengstock, D.-S. L¨uhmann, and P. Schmelcher, Beyond-mean-field study of a binary bosonic mixture in a state-dependent honeycomb lattice , Phys. Rev. A ,043639 (2015).[185] P. Grech and R. Seiringer, The Excitation Spectrum for Weakly Interacting Bosons in a Trap , Commun.Math. Phys. , 559 (2013). 102186] R. Seiringer, The Excitation Spectrum for weakly interacting Bosons , Commun. Math. Phys. , 565(2011).[187] R. Beinke, L. S. Cederbaum, and O. E. Alon, Enhanced many-body effects in the excitation spectrum ofa weakly-interacting rotating Bose-Einstein condensate , arXiv:1803.04762 [cond-mat.quant-gas] (2018).[188] M. R. Matthews, B. P. Anderson, P. C. Haljan, D. S. Hall, C. E. Wieman, and E. A. Cornell, Vorticesin a Bose-Einstein Condensate , Phys. Rev. Lett. , 2498 (1999).[189] N. Regnault and T. Jolicoeur, Quantum Hall Fractions in Rotating Bose-Einstein Condensates , Phys.Rev. Lett. , 030402 (2003).[190] S. Stock, B. Battelier, V. Bretin, Z. Hadzibabic, and J. Dalibard, Bose-Einstein condensates in fastrotation , Laser Phys. Lett. , 275 (2005).[191] C.-C. Chang, N. Regnault, T. Jolicoeur, and J. K. Jain, Composite fermionization of bosons in rapidlyrotating atomic traps , Phys. Rev. A From rotating atomic rings to quantum Hall states , Sci. Rep. , 43 (2011).[193] I. Coddington, P. Engels, V. Schweikhard, and E. A. Cornell, Observation of Tkachenko Oscillations inRapidly Rotating Bose-Einstein Condensates , Phys. Rev. Lett. , 100402 (2003).[194] V. Schweikhard, I. Coddington, P. Engels, V. P. Mogendorff, and E. A. Cornell, Rapidly RotatingBose-Einstein Condensates in and near the Lowest Landau Level , Phys. Rev. Lett. , 040404 (2004).[195] A. Aftalion, X. Blanc, and J. Dalibard, Vortex patterns in a fast rotating Bose-Einstein condensate ,Phys. Rev. A , 023611 (2005).[196] T. Mizushima, M. Ichioka, and K. Machida, Beliaev Damping and Kelvin Mode Spectroscopy of a Bose-Einstein Condensate in the Presence of a Vortex Line , Phys. Rev. Lett. , 180401 (2003).[197] T. P. Simula and K. Machida, Kelvin-Tkachenko waves of few-vortex arrays in trapped Bose-Einsteincondensates , Phys. Rev. A , 063627 (2010).[198] L. O. Baksmaty, S. J. Woo, S. Choi, and N. P. Bigelow, Tkachenko Waves in Rapidly Rotating Bose-Einstein Condensates , Phys. Rev. Lett. , 160405 (2004).[199] F. Chevy, Low-lying twisting and acoustic modes of a rotating Bose-Einstein condensate , Phys. Rev. A , 041604(R) (2006).[200] A. Collin, Collective modes and the broken symmetry of a rotating attractive Bose gas in an anharmonictrap , Phys. Rev. A , 013611 (2006).[201] F. Ancilotto and F. Toigo, Dipolar Bose gas in a highly anharmonic trap , Phys. Rev. A , 023617(2014).[202] S. Viefers, T. H. Hansson, and S. M. Reimann, Bose condensates at high angular momenta , Phys. Rev.A , 053604 (2000).[203] S. M. Reimann, M. Koskinen, Y. Yu, and M. Manninen, Rotating quantum liquids crystallize , New J.Phys. , 59 (2006).[204] J. C. Cremon, G. M. Kavoulakis, B. R. Mottelson, and S. M. Reimann, Vortices in Bose-Einsteincondensates: Finite-size effects and the thermodynamic limit , Phys. Rev. A , 053615 (2013).103205] J. C. Cremon, A. D. Jackson, E. ¨O. Karabulut, G. M. Kavoulakis, B. R. Mottelson, and S. M. Reimann, Rotating Bose-Einstein condensates: Closing the gap between exact and mean-field solutions , Phys. Rev.A , 033632 (2015).[206] P. T. Nam and R. Seiringer, Collective Excitations of Bose Gases in the Mean-Field Regime , Arch.Rational Mech. Anal. , 381 (2015).[207] A. L. Fetter, Rotating trapped Bose-Einstein condensates , Rev. Mod. Phys. , 647 (2009).[208] O. E. Alon, Solvable model of a generic trapped mixture of interacting bosons: reduced density matricesand proof of Bose-Einstein condensation , J. Phys. A , 295002 (2017).[209] S. I. Mistakidis, G. C. Katsimiga, P. G. Kevrekidis, and P. Schmelcher, Correlation effects in thequench-induced phase separation dynamics of a two species ultracold quantum gas , New J. Phys. ,043052 (2018).[210] M. Pyzh, S. Kr¨onke, C. Weitenberg, and P. Schmelcher, Spectral properties and breathing dynamics ofa few-body Bose-Bose mixture in a 1D harmonic trap , New J. Phys. , 015006 (2018).[211] V. Cikojevi´c, L. Vranjeˇs Marki´c, and J. Boronat, Harmonically trapped Bose-Bose mixtures: a quantumMonte Carlo study , New J. Phys. Effective interactions in a quantum Bose-Bose mixture ,Phys. Rev. A , 053617 (2018).[213] P. Cheiney, C. R. Cabrera, J. Sanz, B. Naylor, L. Tanzi, and L. Tarruell, Bright Soliton to QuantumDroplet Transition in a Mixture of Bose-Einstein Condensates , Phys. Rev. Lett. , 135301 (2018).[214] Y. Asano, M. Narushima, S. Watabe, and T. Nikuni, Collective excitations in Bose-Fermi mixtures ,arXiv:1807.06747 [cond-mat.quant-gas] (2018).[215] J. Chen, J. M. Schurer, and P. Schmelcher, Bunching-antibunching crossover in harmonically trappedfew-body Bose-Fermi mixtures , Phys. Rev. A , 023602 (2018).[216] J. Chen, J. M. Schurer, and P. Schmelcher, Entanglement Induced Interactions in Binary Mixtures , Phys.Rev. Lett. , 043401 (2018).[217] P. Siegl, S. I. Mistakidis, and P. Schmelcher, Many-body expansion dynamics of a Bose-Fermi mixtureconfined in an optical lattice , Phys. Rev. A , 053626 (2018).[218] R. S. Lous, I. Fritsche, M. Jag, F. Lehmann, E. Kirilov, B. Huang, and R. Grimm, Probing the Interfaceof a Phase-Separated State in a Repulsive Bose-Fermi Mixture , Phys. Rev. Lett. , 243403 (2018).[219] S. H¨afner, J. Ulmanis, E. D. Kuhnle, Y. Wang, C. H. Greene, and M. Weidem¨uller, Role of the in-traspecies scattering length in the Efimov scenario with large mass difference , Phys. Rev. A , 062708(2017).[220] O. E. Alon, A. I. Streltsov, and L. S. Cederbaum, Unified view on multiconfigurational time propagationfor systems consisting of identical particles , J. Chem. Phys. , 154103 (2007).[221] E. Fasshauer and A. U. J. Lode, Multiconfigurational time-dependent Hartree method for fermions:Implementation, exactness, and few-fermion tunneling to open space , Phys. Rev. A , 033635 (2016).[222] C. L´evˆeque and L. B. Madsen, Time-dependent restricted-active-space self-consistent-field theory forbosonic many-body systems , New. J. Phys. , 043007 (2013).104223] C. L´evˆeque and L. B. Madsen, Multispecies time-dependent restricted-active-space self-consistent-fieldtheory for ultracold atomic and molecular gases , J. Phys. B , 155302 (2018).[224] E. Polizzi, Density-Matrix-Based Algorithms for Solving Eigenvalue Problems , Phys. Rev. B. Vol. Proof of Bose-Einstein Condensation for Dilute Trapped Gases , Phys. Rev.Lett. , 17 (2002).[227] S. Klaiman and L. S. Cederbaum, Overlap of exact and Gross-Pitaevskii wave functions in Bose-Einsteincondensates of dilute gases , Phys. Rev. A , 063648 (2016).[228] L. S. Cederbaum, Exact many-body wave function and properties of trapped bosons in the infinite-particlelimit , Phys. Rev. A , 013615 (2017).[229] P. A. Ruprecht, M. Edwards, K. Burnett, and C. W. Clark, Probing the linear and nonlinear excitationsof Bose-condensed neutral atoms in a trap , Phys. Rev. A , 4178 (1996).[230] Y. Castin, and R. Dum, Instability and Depletion of an Excited Bose-Einstein Condensate in a Trap ,Phys. Rev. Lett. , 3553 (1997).[231] C. W. Gardiner, Particle-number-conserving Bogoliubov method which demonstrates the validity of thetime-dependent Gross-Pitaevskii equation for a highly condensed Bose gas , Phys. Rev A , 1414 (1997).[232] Y. Castin, and R. Dum, Low-temperature Bose-Einstein condensates in time-dependent traps: Beyondthe U(1) symmetry-breaking approach , Phys. Rev. A , 3008 (1998).[233] O. E. Alon, A. I. Streltsov, and L.S. Cederbaum, Time-dependent multi-orbital mean-field for fragmentedBose-Einstein condensates , Phys. Lett. A , 453 (2007).[234] O. E. Alon, A. I. Streltsov, and L. S. Cederbaum, Multiorbital mean-field approach for bosons, spinorbosons, and Bose-Bose and Bose-Fermi mixtures in real-space optical lattices , Phys. Rev. A , 013611(2007).[235] L. S. Cederbaum and A. I. Streltsov, Best mean-field for condensates , Phys. Lett. A , 564 (2003).[236] L. S. Cederbaum and A. I. Streltsov, Self-consistent fragmented excited states of trapped condensates ,Phys. Rev. A , 023610 (2004).[237] O. E. Alon and L. S. Cederbaum, Pathway from Condensation via Fragmentation to Fermionization ofCold Bosonic Systems , Phys. Rev. Lett. , 140402 (2005).[238] O. E. Alon, A. I. Streltsov, and L. S. Cederbaum, Zoo of Quantum Phases and Excitations of ColdBosonic Atoms in Optical Lattices , Phys. Rev. Lett. , 030405 (2005).[239] L. S. Cederbaum, A. I. Streltsov, and O. E. Alon, Fragmented Metastable States Exist in an AttractiveBose-Einstein Condensate for Atom Numbers Well above the Critical Number of the Gross-PitaevskiiTheory , Phys. Rev. Lett. , 040402 (2008).[240] J. Grond, A. I. Streltsov, and L. S. Cederbaum, Excitation spectra of fragmented condensates by linear-response: General theory and application to a condensate in a double-well , Phys. Rev. A , 063607(2012).[241] O. E. Alon, Many-body excitation spectra of trapped bosons with general interaction by linear response ,J. Phys.: Conf. Ser. , 012039 (2015). 105242] S. Klaiman and O. E. Alon, Variance as a sensitive probe of correlations , Phys. Rev. A , 063613(2015).[243] S. Klaiman, A. I. Streltsov, and O. E. Alon, Uncertainty product of an out-of-equilibrium many-particlesystem , Phys. Rev. A , 023605 (2016).[244] S. Klaiman, A. I. Streltsov, and O. E. Alon, Uncertainty product of an out-of-equilibrium Bose-Einsteincondensate , J. Phys.: Conf. Ser. Impact of the range of the interaction on the quantum dynamics of abosonic Josephson junction , Chem. Phys. , 72 (2018).[246] S. Klaiman, R. Beinke, L. S. Cederbaum, A. I. Streltsov, and O. E. Alon, Variance of an anisotropicBose-Einstein condensate , Chem. Phys.509