Many-body excitations and de-excitations in trapped ultracold bosonic clouds
MMany-body excitations and de-excitations in trapped ultracold bosonic clouds
Marcus Theisen ∗ and Alexej I. Streltsov † Theoretische Chemie, Physikalisch-Chemisches Institut, Universität Heidelberg,Im Neuenheimer Feld 229, D-69120 Heidelberg, Germany (Dated: September 26, 2018)We employ the MultiConfiguraional Time-Dependent Hartree for Bosons (MCTDHB) method tostudy excited states of interacting Bose-Einstein condensates confined by harmonic and double-welltrap potentials. Two approaches to access excitations, a static and a dynamic one, have been studiedand contrasted. In static simulations the low-lying excitations have been computed by utilizing theLR-MCTDHB method - a linear response theory constructed on-top of a static MCTDHB solution.Complimentary, we propose two dynamic protocols that address excitations by propagating theMCTDHB wave-function. In particular, we investigate dipole-like oscillations induced by shiftingthe origin of the confining potential and breathing-like excitations by quenching frequency of aparabolic part of the trap. To contrast static predictions and dynamic results we have computedtime-evolutions and their Fourier transforms of several local and non-local observables. Namely, westudy evolution of the (cid:104) x ( t ) (cid:105) , its variance Var( x ( t )) , and of a local density computed at a selectedposition. We found out that the variance is the most sensitive and informative quantity - along withexcitations it contains information about the de-excitations even in a linear regime of the induceddynamics. The dynamic protocols are found to access the many-body excitations predicted by thestatic LR-MCTDHB approach. ∗ [email protected] † [email protected] a r X i v : . [ c ond - m a t . qu a n t - g a s ] A ug I. INTRODUCTION
Many-body approaches are crucial for a modern understanding and physical description of ultracold quantum gases.Particularly, systems of strongly interacting ultracold atoms demand for theories that predict dynamic behavior andexcitation spectra [1–4]. A system of major interest is the Bose-Einstein condensate (BEC). With its experimentalrealization [5, 6] a well-controllable mesoscopic quantum object is at hand offering ways to study many-body systemsand excitations therein empirically. Experimentally, owing to their controllability, excitations of BECs have beenwidely explored [7–10]. At the same time tremendous theoretical efforts seek to account for the observed correlationphenomena and their proper many-body descriptions [1–4]. Since analytic solutions to many-body problems are rare,numerical simulations are inevitable.The most popular mean-field theory to treat BECs in dilute, ultracold vapors is based upon the Gross-Pitaevskii(GP) equation [4]. Within the GP approach all constituting bosons are assumed to occupy a single one-particlestate, describing thereby a simple, i.e., fully condensed and coherent condensate [11]. Clearly, whenever the BEC ispartially condensed, depleted, or multi-fold fragmented [12] a single wave-function approach is insufficient. On physicalside all these phenomena originate from strong interparticle interactions and multi-well topologies of the confiningtraps. Experimentally, connections between strong interactions and quantum depletion have been established in [13],the fragmentation phenomena in double-well traps have been observed in [14, 15]. First experimental realization ofthe Mott insulator states in multi-well traps formed by optical lattices has been reported in [16]. On theoreticalside to grasp depletion and fragmentation phenomena one has to go beyond the one-mode GP ansatz and use morecomplicated many-body theories where the bosonic wave-function is constructed with multiple (one-particle) orbitalmodes [17–21]. In this paper to describe condensation, depletion and fragmentation phenomena in static setups andtime-dependent processes on the same ground we use the Multi-Configurational Time-Dependent Hartree method forBosons (MCTDHB) [22, 23], which is available within the MCTDHB-Laboratory package [24].The standard static, i.e., time-independent approach to calculate excitation spectra of BECs at the one-orbitalmean-field level is to apply the Linear Response (LR) theory to the GP ground state. The underlying equationsare often called Bogoliubov-de Gennes (BdG ≡ LR-GP) equations [25]. Recently, the same linear-response idea wasapplied to the multi-orbital best mean-field [26] and many-body MCTDHB [27] ground states allowing thereby toaccess excitations of the multi-fold fragmented BECs at the mean-field and many-body levels. The invented [27],developed [24], generalized [28] and benchmarked [29] LR-MCTDHB theory has allowed us to discover low-lyingexcitations in trapped systems which are not described within the BdG theory. As one might expect such states arepresent in fragmented systems [26] and, surprisingly, also in condensed systems [27–29] where BdG was believed togovern the physics.Naturally, questions of interest concern the origin of many-body excited states, their properties and possible classifi-cation schemes. Here, we attempt to generalize characterization of the excitations obtained in the non-interacting caseto the interacting systems. We also distinguish two kinds of many-body excitations - one branch describes collectiveexcitations while the second one represents excitations involving multiple particles. These excitations will be definedand described in some detail below.Another issue tackled in this paper is a dynamic control. The long-term perspective is to find dynamic protocolsthat allow for a controlled population of the desired excitations. Within the framework of the MCTDHB method thisproblem can, in principle, be solved by applying the optimal-control theory as proposed in [30–33] and, more recently,by merging the general optimal-control Chopped RAndom Basis (CRAB) algorithms [34, 35] with the MCTDHBmethod [36]. Here, however, we present and consider two simple protocols that excite BECs via sudden modificationof the system’s Hamiltonian. Such manipulations manifest in oscillation of the many-body wave-function, and thusthe particle density, in time. In an experimental setup changes on the Hamiltonian can be realized by introducingcontrol fields, altering the confining potential or toggling the inter-particle interaction using Feshbach resonances [37].An open question is how to detect and verify many-body excitations confidently. A first step towards the answer isto find out which observables provide more detailed information about the many-body excitations. To this end, wepropose to use the positional variance as a sensitive probe for many-body excitations and de-excitations.The structure of this paper is as follows. In Sec. II we introduce the systems of interest by giving the Hamiltonianand specifying the external trapping potential. Sec. III contains numerical results of BECs described at mean-field(GP) and many-body (MCTDHB) levels. Two methods to analyze excitation spectra are discussed and compared:A static approach using linear response theory (Sec. III A) and a dynamic approach using wave-packet propagation(Sec. III B). Finally, in Sec. IV we summarize the results and give future prospects.
II. SYSTEM AND HAMILTONIAN
The MCTDHB( M ) algorithm [22, 23] effectively solves the time-dependent many-body Schrödinger equation ˆ H Ψ = i (cid:126) ∂∂t Ψ , (1)taking the sum over symmetrized Hartree products (or permanents ) as the ansatz for the many-body wave-function | Ψ( t ) (cid:105) = (cid:88) (cid:126)n C (cid:126)n ( t ) | (cid:126)n ; t (cid:105) . (2)Here, the summation runs over all possible configurations of the state vector (cid:126)n = ( n , n , . . . , n M ) that preserve thetotal number of bosons N = (cid:80) Mi n i and M denotes the number of one-particle wave-functions φ i ( r , t ) (or orbitals )used to construct the respective permanents.The shape of the orbitals φ i ( r , t ) and the expansion coefficients C (cid:126)n ( t ) , which account for normalization, arevariational time-dependent parameters of the MCTDHB method. We would like to stress that the single orbitalMCTDHB( M = 1 ) approach is fully equivalent to the Gross-Pitaevskii theory. Using the MCTDHB( M ) method with M ≥ orbitals implies an above mean-field, many-body treatment.Given the MCTDHB wave-function at every time-slice t we can construct and diagonalize the reduced one-bodydensity matrix ρ ( r , r (cid:48) ; t ) . It can be expressed in terms of the natural orbitals and occupations, i.e., its eigenstates φ NOi ( r , t ) and eigenvalues n i ( t ) : ρ ( r , r (cid:48) ; t ) ≡ N (cid:90) Ψ ∗ ( r (cid:48) , r , . . . , r N , t ) Ψ( r , r , . . . , r N , t ) d r · · · d r N = M (cid:88) i,j ρ ij ( t ) φ ∗ j ( r (cid:48) , t ) φ i ( r , t )= M (cid:88) i n i ( t ) φ ∗ NOi ( r (cid:48) , t ) φ NOi ( r , t ) . (3)The many-body state is called condensed if only one natural eigenfunction – condensate orbital is ≈ occupied.The BEC may be called depleted when the non-condensed fraction, i.e., the total occupation of other than thecondensate orbital ( (cid:80) Mi> n i ) is of the order of . When several eigenfunctions have macroscopic occupation thesystem is called fragmented. For completeness we also present the expression for the two-body density, which we willuse later: ρ ( r , r , r (cid:48) , r (cid:48) ; t ) ≡ N ( N − (cid:90) Ψ ∗ ( r (cid:48) , r (cid:48) , . . . , r N , t ) Ψ( r , r , . . . , r N , t ) d r · · · d r N = M (cid:88) i,j,k,l ρ ijkl ( t ) φ ∗ j ( r (cid:48) , t ) φ ∗ l ( r (cid:48) , t ) φ i ( r , t ) φ k ( r , t ) . (4)For N interacting particles the Hamiltonian reads ˆ H = N (cid:88) j =1 ˆ h ( r j ) + N (cid:88) j 26 Hz and T = 1 . 37 ms , respectively.In this paper we focus on one-dimensional (1D) systems composed of N = 10 bosons with contact interaction W ( R ) = δ ( R ) . The bosons are trapped by a harmonic potential with a Gaussian barrier: V ( x ) = ax + b exp( − cx ) . (7)Let us discuss some properties of this trapping potential. Obviously, for b = 0 it reduces to a harmonic potentialwith frequency ω H = √ a . Furthermore, for general a, b, c ∈ R + and k = bc/a > the potential is a double-well withminima at x , = ±√ c − ln k . b=0n =99.9% t o t a l e n e r g y b=5n =95.3%spatial coordinate b=10n =60.1% FIG. 1. The three regimes studied: Shapes of the trap potential (red lines), its eigenvalues (black lines) and wave-packetdensity (grey) are illustrated. The potentials (7) are normalized as to vanish in the minima. For the double-wells the harmonicapproximation at the potential minimum is shown (green potentials, blue eigenvalues). Occupation numbers of the first naturalorbital n and densities of the ground state are calculated for N = 10 particles with Λ = λ ( N − 1) = 1 at the MCTDHB( M = 2 )level. Typically, we set a = 1 / ( ω H = 1 ) and c = 1 and distinguish between three major regimes according to thebarrier height and the population of the most occupied (first) natural orbital: i) a fully condensed BEC in a harmonicpotential ( b = 0 ), ii) a depleted BEC in a shallow double-well trap ( b = 5 ) and iii) an almost fully fragmented BECin a deep double-well potential ( b = 10 ). The regimes are sketched in Fig. 1.Since the chosen potentials are symmetric to the origin, the Hamiltonian Eq. (5) is reflection invariant. This allowsa separation of the underlying Hilbert space into symmetric and an anti-symmetric manifolds. We call states livingin the symmetric and anti-symmetric subspaces gerade ( g ) and ungerade ( u ), respectively. III. RESULTSA. Static picture In this section we analyze the excitation spectra of bosonic systems by virtue of the LR theory. The procedureinvolves two steps [27–29]. First, we compute the ground state E GS for a given system using MCTDHB( M ) method.Second, we compute the excitation spectrum by applying the LR-MCTDHB( M ) atop this static solution. The numberof orbitals M determines the quality of the ground state. Hence, the more orbitals M used to describe the groundstate, the greater the variety of excited states. In particular, we distinguish between excitations that appear due todeformation of the involved orbitals φ i ( orbital-like ) and excitations that originate to the redistribution of the particlesbetween given orbitals via configuration interaction ( CI-like ) C (cid:126)n , see Eq. (2) and Refs. [27–29] for more details. Foreach excitation the response amplitude to a given perturbation ( ∝ x or ∝ x ) is calculated. Its magnitude determinesthe excited state’s intensity in the respective excitation spectrum. Physically it expresses the probability to excite thestate by that specific perturbation. For a detailed description of the LR-MCTDHB theory and its implementationsee [27–29].We exploit the separability of the Hilbert space by calculating the response of the system to small perturbationsof ungerade ( ∝ x ) and gerade ( ∝ x ) symmetry. This allows us to uniquely identify each excitation from the groundstate by its energy ∆ E = E − E GS and its symmetry (either u or g ).In this work we consider the gradual transformation of a harmonic potential into a deep double-well potential byintroducing a Gaussian barrier. For weakly interacting bosons the excitation spectrum is shown in Fig. 2. The standardLR-GP ≡ BdG ≡ LR-MCTDHB( M = 1 ) results are indicated by black circles and the many-body LR-MCTDHB( M =2 ) results by colored lines, where green and red colors mark gerade and ungerade states, respectively. Δ E b 0 1 2 3 4 0 2 4 6 8 101u1g' 0th1st Δ E b FIG. 2. Weak interaction: Excitation spectrum as a function of barrier height b calculated for N = 10 bosons and Λ = λ ( N − 1) = 0 . . ∆ E is the excitation energy with respect to the ground state. Mean-field, LR-GP ≡ BdG, results aregiven by open black circles and many-body, LR-MCTDHB(2), predictions by colored lines. Green lines indicate gerade and redlines ungerade symmetry. States are labeled by their nodal structure and symmetry; primes are used to mark pure many-bodyexcitations. See text for discussion. From Fig. 2 it is clear that there are states which are reproduced at both the mean-field LR-GP ≡ LR-MCTDHB( )and the many-body LR-MCTDHB( ) levels. However, one can also see states that are only present in the many-body treatment. Following Refs. [27–29] we classify the obtained solutions according to this property further into mean-field-like and many-body excitations.Let us discuss the differences between mean-field-like and pure many-body excitations. For a non-interactingsystem the accessible states can be visualized in terms of particle-hole excitations involving eigenstates (orbitals) ofthe respective trap potential. Within LR-GP (BdG) all bosons reside in the ground state and only h − p excitationsare allowed. They describe the excitation of a single particle ( p ) to a higher-order one-particle state leaving a hole ( h )in the ground state. The order of a one-particle function is given by its nodal structure, where an even (odd) numberof nodes corresponds to an excitation of gerade (ungerade) symmetry. At the many-body level excitations may havea more complicated structure involving multiple particle excitations at the same time.For the harmonic regime ( b = 0 in Fig. 2), excitations have been studied elsewhere [27, 29] by decoupling relativeand center-of-mass (CoM) motion. Essentially, there are mean-field-like and many-body states of equal symmetrythat are energetically close. As suggested there, we label these excitations according to their nodal structure andthe symmetry they exhibit. For example, u is the lowest-in-energy state with a one-node structure. It describesthe ungerade CoM excitation in a harmonic oscillator at ∆ E = ω H = 1 . The higher-order h − p states g, u, . . . can correspond to the excitations of the CoM and relative motion. Pure many-body states are indicated by a prime,e.g., g (cid:48) . As shown in [27], the lowest-in-energy many-body excitations correspond to higher-order CoM modes. It isworthwhile to mention that in harmonically trap systems with a non-contact, e.g., parabolic interparticle interaction,the lowest-in-energy many-body states can correspond to the higher-order excitations of the relative motion [29].Next, we would like to discuss the almost fully fragmented regime in a deep double-well trap ( b = 10 in Fig. 2). Oneclearly observes the formation of band-like structures for the asymptotic limit where the energy difference betweeneach pair (green and red curves) vanishes. In this limit of high barriers the entire excitation spectrum becomestwo-fold degenerate and the ground state of the many-body system is two-fold fragmented.In the deep double-well regime the LR-GP theory predicts only one excited state in the lowest band depicted inFig. 2 by open black circles. This state corresponds to the h − p excitation from the ground state GP orbital of thegerade symmetry to its complimentary ungerade one. Similar 1h-1p excitations to higher quasi-degenerate excitedorbitals form the mean-field description of the higher bands. By contrasting these mean-field predictions with the LR-MCTDHB( M = 2 ) results plotted in Fig. 2 by green and red lines we see that in the many-body description additionalstates contribute to the bands. In the many-body picture the zeroth band is composed of ten states, equal to the totalnumber of particles: N = 10 . To understand this, let us once again return to the non-interacting case and recall thattwo lowest-in-energy eigenstates of the deep double-well trap are almost degenerate two-hump symmetric ( gerade )and an anti-symmetric ( ungerade ) orbitals. In the non-interacting ground state all N bosons occupy the geradefunction. The first single-particle excitation to the other orbital is described by u , changing, thereby, the overallsymmetry. This h − p excitation is the same as within the BdG theory. Next excitation corresponds to the situationin which two particles are excited simultaneously from the ground state orbital to its anti-symmetric counterpart. Asa consequence the overall gerade symmetry of this many-body excited state marked as g (cid:48) is re-established. A threeparticle excitation process is described by u (cid:48) and so on. This branch of many-body excitations thereby correspondsto the successive population of the respective anti-symmetric orbital. Formally, all these excitations can be describedby different redistributions (configurations) of the particles among the two orbitals. Therefore, these excitations areCI-like.In a presence of interactions the many-body states of the lowest (0th) band become split. However, this branch ofexcitations is hard to access, due to the involvement of multi-particle processes. This physical conclusion is confirmedby the values of the corresponding linear response amplitudes, see Tables I, II and III and their discussions.In contrast to the lowest band, bands of higher energy at the LR-MCTDHB(2) level of description contain onlyfour states. Closer examination of the first band shows that two states are GP-like and the other two are pure many-body excitations. Depending on the orbital in which the h − p excitation is realized, the GP-like state is eitherof gerade or ungerade symmetry. The corresponding many-body excited states have a h − p − p structure: Inaddition to the GP-like excitation a second particle is excited to a higher-in-energy orbital. At the LR-MCTDHB(2)level of description more complex excitation structures like h − p − p , h − p − p are not reachable in higherbands. This absence of the higher-order excitations reflects the fact that the MCTDHB(2) is an approximation tothe exact solution and, therefore, higher-level theories are required for their proper descriptions, see e.g., discussionof the LR-MCTDHB(4) results.For dynamic studies the clustering and degeneracy of states in bands complicates identification and attributionof the distinct excitations. One way to overcome this issue is to investigate BECs in a shallow double-well ( b = 5 in Fig. 2), where the bands are not yet formed. Another possibility to split the excited states is to increase theinterparticle interaction strength. In Fig. 3 we show the same study as before but for stronger interparticle interactionstrength Λ = λ ( N − 1) = 1 . For convenience, we keep the introduced labeling.Let us compare and contrast the excitation spectra at weak (Fig. 2) and intermediate (Fig. 3) interparticle inter-actions. At the mean-field level (BdG - black circles) the overall spectra are quite similar. The mean-field excitationenergies are slightly shifted due to interparticle interaction. All these mean-field excitations are reproduced at themany-body level (LR-MCTDHB( ) - color lines) in both cases. The main difference to the weak interaction, however,is a splitting of states in the zeroth band in the deep double-well regime. At intermediate interactions the ultracoldsystems become less coherent i.e., more correlated and a single wave-function description is definitely insufficient.By contrasting LR-MCTDHB( ) and LR-MCTDHB( ) (black dots) results depicted in Fig. 3 we see that, asexpected, at higher energies the four-orbital theory introduces additional many-body states and shifts positions ofsome excited states. The observed shifts suggest that an even higher-level MCTDHB theory with more than M = 4 Δ E b 0 1 2 3 4 0 2 4 6 8 10 Δ E b 0 1 2 3 4 0 2 4 6 8 10 Δ E b FIG. 3. Intermediate interaction: Excitation spectrum as a function of barrier height b for N = 10 bosons and Λ = λ ( N − 1) = 1 . ∆ E is the excitation energy with respect to the ground state. Open black circles depict the mean-fieldLR-GP ≡ BdG results. The two-orbital LR-MCTDHB( M = 2 ) many-body results are plotted by colored lines, four-orbitalLR-MCTDHB( M = 4 ) many-body results – by black dots. Green lines indicate gerade and red lines ungerade symmetry. Thelow-energy excitations are seen to converge. See text for the discussion. orbitals is needed to converge the respective many-body states. However, for the excitation energies smaller than ∆ E (cid:46) and not to small barriers ( b (cid:38) ) the deviations between two- and four-orbital LR-MCTDHB results arenegligible. We conclude that the usage of the two-orbitals LR-MCTDHB( ) theory is sufficient to describe the low-lying excitations for all the here studied systems and in all the studied regimes.We would like to mention that with stronger interaction strength avoided crossings become visible. They reflect thefact that states of same symmetry do not cross but rather exchange properties. An example is the avoided crossingat b ≈ and ∆ E ≈ in Fig. 3.To conclude the static part, the LR-MCTDHB( M ) theory has been used to compute and study systematicallyexcitation spectra of the interacting bosonic systems trapped in different external potential. By introducing a Gaussianbarrier b the harmonic potential has been gradually transformed into a deep double-well potential. The computedexcitation energies ∆ E and intensities I for the three important regimes ( b = 0 , , ) are listed in Tables I, II andIII. Next, we plan to contrast these static results with the dynamic studies in the following sections. B. Dynamics In this section we solve the many-body Schrödinger equation, Eq. (1), using MCTDHB( M ) wave-packet propagation.We use the following protocol to access excited states: Firstly, we obtain a given system in its ground state. Secondly,we manipulate the Hamiltonian Eq. (5) and propagate the above obtained state in an altered environment for twobasic scenarios, namely shift and quench . In order to activate ungerade excitations, which correspond to the dipole -likemovement of the BEC in a harmonic trap we propose to shift the origin x → x + x shift of the external potential inEq. (7). The gerade ( breathing -like) excitations of the BEC can be activated by quenching the frequency of the outerharmonic potential a -parameter in Eq. (7). These sudden manipulations of the Hamiltonian induce a non-equilibriumevolution of the system which manifests itself as oscillation patterns of the one-particle density and in the temporal FIG. 4. Propagation: Right panel: Breathings in a harmonic potential induced by a quench of a from . to . . Leftpanel: Dipole-like oscillation in a deep double-well after shifting the potential by x shift = 0 . . Calculations are done atMCTDHB( M = 2 ) level for N = 10 particles and Λ = λ ( N − 1) = 1 . evolution of other observables. We shall see that the proper choice of observables is very crucial as each one is sensitiveto different aspects of the excitations and, therefore, provides different information about them.Let us start our discussion on suitable observables with the particle density which is computed as the diagonal partof the reduced one-body density Eq. (3): ρ ( x, t ) ≡ ρ ( x, x ; t ) = M (cid:88) i =1 n i ( t ) φ ∗ NOi ( x, t ) φ NOi ( x, t ) . (8)In Fig. 4 the typical evolution of ρ ( x, t ) in the shift and quench scenarios is shown. Here, time is plotted againstspace ( x coordinate) in a Minkovskii-like manner and the particle density ρ ( x, t ) is indicated by a color map: Yellowand black colors mark regions of high and low density, respectively. The left panel shows the dipole-like movement ofa BEC after shifting the origin of the double-well potential. The right panel shows typical evolution of the system inthe harmonic trap activated by the quench scenario. The visible symmetric deformation of the density is known as abreathing mechanism [38, 39].In addition to the particle density, we also analyze evolution of the expectation values of several one- and two-bodyoperators. The first one is the position operator defined for N -particle case as ˆ x = (cid:80) Ni x i . Its expectation valuereads: (cid:104) x (cid:105) = M (cid:88) i,j ρ ij ( t ) x ij , (9)where ρ ij are elements of the reduced density matrix, see Eq. (3) and x ij = (cid:82) φ ∗ j ( x, t ) x φ i ( x, t )d x . The second one isthe square of the position operator, which contains contributions from one- and two-body operators ˆ x = ( (cid:80) Ni x i ) = (cid:80) Ni x i + (cid:80) Ni,j x i x j . So, computation of its expectation value is a more involved task [40]: (cid:104) x (cid:105) = M (cid:88) i,j ρ ij ( t ) x bij + M (cid:88) i,j,k,l ρ ijkl ( t ) x bijkl , (10)where x bij = (cid:82) φ ∗ j ( x, t ) x φ i ( x, t ) d x and x bijkl = (cid:82) φ ∗ j ( x , t ) φ ∗ l ( x , t ) x x φ i ( x , t ) φ k ( x , t ) d x d x . Finally, we cancompute the corresponding variance Var( x ) = (cid:104) x (cid:105) − (cid:104) x (cid:105) . (11)All these expectation values are integral quantities and, therefore, contain global information about the system’sevolution: (cid:104) x (cid:105) gives the mean position of the cloud and its variance Var( x ) - the deviation from the mean.In all the defined above expectation values time dependency is assumed implicitly. Let us now derive formal explicitexpressions. The time evolution of any initial state | Ψ( t = 0) (cid:105) can be re-expressed in terms of the exact eigenstates | Ψ n (cid:105) of the Hamiltonian | Ψ( t ) (cid:105) = (cid:88) n a n e − iE n t | Ψ n ( t ) (cid:105) , (12)where a n = (cid:104) Ψ( t = 0) | Ψ n (cid:105) . The evolution of the expectation value of any operator ˆ O , thus, reads (cid:104) Ψ( t ) | ˆ O | Ψ( t ) (cid:105) = (cid:88) n,m a n a m e − i ( E n − E m ) t (cid:104) Ψ m | ˆ O | Ψ n (cid:105) . (13)It is important to stress that along with the differences in energy between the ground and excited states – excitationenergies, it also contains contributions from the differences in energy between other involved excited states. Fromnow on for the sake of simplicity we refer to the contributions from these higher-order processes which do not involvethe ground state as de-excitations and mark them correspondingly, see for example g (cid:48) ↔ g de-excitation in Table I.Their probabilities depend on the initial wave-function via expansion coefficients a n in Eq. (12) and on the particularintegrals (cid:104) Ψ m | ˆ O | Ψ n (cid:105) . A naive expectation is that contributions to the evolution from the de-excitations should benegligible if the activated dynamics keeps the system in a linear response regime. Formally, such a regime can bereached if only a is dominant in the expansion of the initial wave-packet Eq. (12). In the present study we wouldlike to verify it and to find out which expectation values – of the one- or two-body operators are more sensitive tothe dynamics and, therefore, contain more information about the excitations and de-excitations.Complimentary to the integral observables, we would also like to monitor the density at some fixed-point position x ρ ( x = x , t ) , (14)which by definition is a local observable. We expect that the integral quantities can be used to characterize the BECas an entity and the local particle density ρ ( x = x , t ) to grasp possible local features. From this perspective theintegral observable can be regarded as a weighted average of the local quantities computed at all possible positions.Let us mention some obvious implications of our dynamic protocols to the above defined quantities of the interest.Obviously, (cid:104) x (cid:105) = 0 holds for the quench scenario in a symmetric trap potential and, therefore, Var( x ) = (cid:104) x (cid:105) . Inother words, quench scenario can address only gerade excitations. In the shift scenario, in general Var( x ) (cid:54) = 0 and theungerade excitations manifest themselves in oscillations of (cid:104) x (cid:105) . Due to the above symmetry arguments, we plan touse the mean position (cid:104) x (cid:105) to study ungerade excitations in the shift scenario and the variance Var( x ) to study geradeexcitations in the quench scenario.As mentioned above, in the dynamic scenarios excitations manifest themselves in temporal oscillations of the many-body wave-function and observables. Hence, we use the MCTDHB wave-function computed at each propagationtime-step to calculate the corresponding local ρ ( x = x ) and integral (cid:104) x (cid:105) , Var( x ) observables. Next, by applying asubsequent Discrete Fourier Transformation (FT) to each evolving quantity we extract the dominating oscillationfrequencies, the corresponding excitation energies ω ≡ ∆ E and their intensities, see appendix for more computationaldetails. One of the primary goals of the present study is to contrast the static LR predictions obtained in the previoussection with these dynamic FT-extracted excitations and intensities.A delicate issue is the length of propagation. Naturally, a long propagation time results in a large number ofoscillation circles that allows for an accurate extraction of frequencies. The numerical task is to assure a small error inpropagation especially for large quenches and displacements. For our purpose, we shift the trap origin to x shift = 0 . and quench the trap frequency a from . to . in Eq. (7). We assume that these values are in a linear regime suchthat excitations are predictable by the LR theory. See appendix for more information on quantitative characterizationof the linear and non-linear regimes in the studied dynamic scenarios.In the following we study the shift and quench scenarios for the three regimes introduced in Sec. I.0 1. Harmonic Trap We start by analysing the dynamic scenarios in the harmonic potential. In Table I the static LR predictions anddynamic results are contrasted directly. Excitation energies ∆ E and corresponding intensities I (in brackets) arelisted for the different theories and dynamic scenarios. We label the excitations as described in Sec. III A. BdG LR-MCTDHB(2) "shift" (GP) "shift" (MB) "quench" (GP) "quench" (MB)label ∆ E ( I ) ∆ E ( I ) ∆ E ( I ) ∆ E ( I ) ∆ E ( I ) ∆ E ( I )1 u ac ) 1.000 (1.00 ac ) - - g b ,0.64 c ) 1.919 (0.40 b ,0.79 c ) g (cid:48) - 2.055 (0.13) - - - 2.050 (0.53 b ) u (cid:48) - 2.792 (0.02) - - - - u g c ) 3.847 (0.21 c ) g (cid:48) ↔ g - 0.135 - - - 0.130 (0.04 b ) a observed in (cid:104) x (cid:105) b observed in Var( x ) c observed in ρ ( x = 1) TABLE I. Excited states in the harmonic potential: Static and dynamic results for excitation energies ∆ E ≤ . and normalized intensities I (in brackets) are presented at mean-field (GP) and many-body (MB) levels. The labeling is asdescribed in Sec. III A. For the statics the intensities I are normalized for ungerade ( u ) and gerade ( g ) manifolds separately.In the dynamics normalization is done by summing over the intensities of all constituting peaks. States with intensity smallerthan one percent are omitted. We attribute the lowest-in-energy dynamic process with ∆ E = 0 . to a higher-order g (cid:48) ↔ g excitation (de-excitation) because only these two lowest static LR states g (cid:48) , g have appropriate energies ∆ E g (cid:48) − ∆ E g = 0 . and symmetry. Footnotes in the dynamic picture indicate the quantity in which the excitation is observed. The dynamics arecomputed for x shift = 0 . and a quench a = 0 . → . using N = 10 and Λ = λ ( N − 1) = 1 . Before going into a detailed analysis note that all the dynamic excitations and de-excitations extracted from theevolution of the respective observables are predicted by the static LR-MCTDHB theory. We show that de-excitationenergies can be obtained within the static LR-MCTDHB approach by taking differences between appropriate excitedstates. The frequencies obtained within the static and dynamic simulations match nicely implying that it is indeedpossible to excite corresponding states dynamically by the suggested protocols. For a quantitative comparison of thestatic and dynamic results the FT-computed intensities have been normalized over all constituting states of equalsymmetry. It is seen that dynamic results and LR predictions share a general tendency: Excitations that have largeresponse amplitudes in the LR predictions tend to be easily excited by the dynamic scenarios.The two left panels of Fig. 5 present a detailed analysis of the shift scenario in a harmonic trap. They depict theFourier spectra of the mean position (cid:104) x (cid:105) and a local density at fixed position ρ ( x = 1) . Black lines plot mean-fieldMCTDHB( M = 1 ) results, red dashed lines show many-body MCTDHB( M = 2 ) results. We directly see that (cid:104) x (cid:105) and ρ ( x = 1) oscillate with a single frequency ω = 1 at both mean-field and many-body levels. By comparison withthe LR results (Table I), we identify this frequency as the lowest-in-energy ungerade CoM excitation, labeled u .This result is independent of the magnitude of the applied shift and of the number M of orbitals used in the staticand dynamic MCTDHB( M ) simulations confirming thereby the separability of CoM and relative motion in harmonictraps.The right panels of Fig. 5 show the analysis of the quench scenario in the harmonic trap. Now, the local density ρ ( x = 1) and the variance Var( x ) which is a global quantity contain different amount of information. Let us startby analyzing the variance. The Fourier spectrum of Var( x ) at the mean-field GP level (black line) reveals a singlepeak attributed to the lowest gerade excited state g . Within the many-body MCTDHB(2) treatment (red dashedlines) this picture changes drastically. In addition to the g excitation a second peak corresponding to the puremany-body excitation g (cid:48) is visible. It is the most intense peak in the spectrum of the variance (see Table I) statingthat a many-body treatment is by all means inevitable for an accurate description of the dynamics. The lowest-in-energy dynamic process seen as a small peak at low frequency with ∆ E ≈ . can be attributed to a higher-order g (cid:48) ↔ g excitation (de-excitation) because only these two lowest-in-energy states g (cid:48) , g in the static LR picturehave appropriate symmetry and energetics E g (cid:48) ↔ g ≡ E g (cid:48) − E g = ∆ E g (cid:48) − ∆ E g = 0 . . The occurrence of thisde-excitation might indicate that we are beyond the linear response regime.The local density ρ ( x = 1) shows a different picture than the integral observable Var( x ) . On the one hand the1 ω i n t e n s it y [ a r b . un it s ] shift quench ρ (x=1) 1 2 3 4frequency ω i n t e n s it y [ a r b . un it s ] shift quench ρ (x=1) 0 0.5 frequency ω i n t e n s it y [ a r b . un it s ] shift quench Harmonic trap: The left frame shows the frequency analysis of an oscillating BEC after a sudden shift x shift = 0 . of the confining harmonic potential. The right frames shows the study of the quench scenario a = 0 . to a = 0 . . The GP (blacklines) and MCTDHB(2) (red dashed lines) results are shown. The intensities are normalized by the sum of all constitutingpeaks. many-body excitation g (cid:48) is absent and on the other hand the higher-order mean-field excitation g becomes visible.The former becomes clear when thinking of g (cid:48) as an excitation involving several particles. Since ρ ( x = 1) has alocal character it simply cannot see this global excitation. However, this locality is the reason for the pronouncedappearance of g corresponding to a higher mean-field excitation. First, notice that the density ρ ( x = x ) is high inregions close to the potential minimum x = 0 and low in outer regions. Second, recall that integral operators e.g., Var( x ) yield the globally most intense excitations by integration (weighted sum) over the local densities ρ ( x = x ) taken at all possible positions. Since g is visible in ρ ( x = 1) but not in Var( x ) , it gives negligible contribution to the2global dynamics but important to the local dynamics. From this observation we conclude that local quantities takenin the regimes of low densities can be more sensitive to higher excitations than the global ones.The dynamic studies in the harmonic potential reveal that many-body excitations play a crucial role even in thecase of a fully condensed system. Many-body aspects of BECs, e.g., g (cid:48) excitation can be observed in the proposedquench scenario by monitoring the evolution of the global operator Var( x ) . We conclude that the variance beingvery informative measure of the many-body correlations [40, 41] also constitutes a sensitive probe for many-bodyexcitations. The fact that de-excitations are visible implies that we are out of the linear response regime, see appendix.We also note that information about higher order mean-field excitations become visible in the local density ρ ( x = x ) taken at regions of comparable low density. In other words, local quantities can probe excitations that are negligibleon a global scale. 2. Shallow Double-Well As for the harmonic potential we analyze the integral (cid:104) x (cid:105) , Var( x ) and local ρ ( x = 1) observables for the shiftand quench scenarios now in the shallow double-well. Figure 6 depicts the discrete FT-spectra extracted from theevolution of these quantities. The most striking observation is that a plethora of low-lying excitations become visible.Here we would like to stress that the ground state of this system is slightly depleted – fragmentation ratio is ≈ . BdG LR-MCTDHB(2) "shift" (GP) "shift" (MB) "quench" (GP) "quench" (MB)label ∆ E ( I ) ∆ E ( I ) ∆ E ( I ) ∆ E ( I ) ∆ E ( I ) ∆ E ( I )0 u a ,0.08 c ) 0.164 (0.25 a ,0.08 c ) - - g (cid:48) - 0.308 (0.02) - - - 0.305 (0.29 b ) u (cid:48) - 0.456 (0.03) - - - - g b ,0.89 c ) 1.654 (0.32 b ,0.91 c ) u (cid:48) - 1.799 (0.15) - 1.782 (0.15 a ,0.12 c ) - - u a ,0.74 c ) 1.983 (0.58 a ,0.69 c ) - - g (cid:48) - 2.134 (0.09) - - - 2.129 (0.13 b ) g b ,0.05 c ) 3.063 (0.05 b ,0.03 c ) u a ,0.14 c ) 3.771 (0.02 a ,0.11 c ) - - g (cid:48) ↔ g - 0.480 - - - 0.476 (0.08 b ) g ↔ g (cid:48) - 1.346 - - - 1.350 (0.04 b ) g ↔ g b ,0.02 c ) 1.408 (0.01 b ) u ↔ u c ) - - - a observed in (cid:104) x (cid:105) b observed in Var( x ) c observed in ρ ( x = 1) TABLE II. Excited states in the shallow double-well: The structure of the table is as in Table I. In addition low-lyingexcitations with intensities smaller than one percent are neglected. These are, in particular, higher-order excitations of thezeroth band, compare to Fig. 2. As seen in Table II every excitation energy obtained in the dynamic study can be attributed and identified with thecorresponding static LR result. Here we adopt the labeling from the deep double-well as described in Sec. III A. Itis worthwhile to mention that the static LR amplitudes computed as responses of the system to small perturbationsof ungerade ( ∝ x ) and gerade ( ∝ x ) symmetry in the double-well traps can give only a rough estimate for theexcitability of the corresponding states by the dynamic protocols. For a more quantitative characterizations a propermodification of the perturbation operators is required, which is out of the scope of the present study.Before going into details let us discuss several important aspects of the shift scenario in the double-well traps. Firstof all there is no separation of the center-of-mass and relative motions, all these states are coupled and the shiftscenario can excite them. The individual wells constituting the double-well are asymmetric, contrast in the Fig. 1the red and green curves depicting the double-well traps and the harmonic approximations to each of these wellscorrespondingly. So, the shift of the origin of the double-well can lead to the excitations of all local modes of theindividual wells.3 ω i n t e n s it y [ a r b . un it s ] shift quench ρ (x=1) 1 2 3 4frequency ω i n t e n s it y [ a r b . un it s ] shift quench ρ (x=1) 0 0.5 frequency ω i n t e n s it y [ a r b . un it s ] shift quench Shallow double-well trap: The description is as for Fig. 5. The left panels of Fig. 6 show the studies in the shift scenario. As can be seen, the many-body spectra include theGP results and introduce the additional many-body excitation with pronounced intensity attributed in the Table II as u (cid:48) . Surprisingly, this state is seen in the evolution of both the global operator (cid:104) x (cid:105) and local quantity ρ ( x = 1) . Thefact that u (cid:48) is also visible in local ρ ( x = 1) indicates that relative and CoM motion are indeed coupled as opposedto the harmonic case.The right panels of Fig. 6 show the FT analysis of the Var( x ) and ρ ( x = 1) in the quench scenario. Concerningthe variance, the many-body treatment in addition to the mean-field GP peaks reveals two quite intense many-bodyexcitations g (cid:48) and g (cid:48) , see Table II. The most-intense mean-field g peak and these two many-body excitations havecomparable intensities stressing the importance of a beyond mean-field treatment in shallow double-wells. We alsofind that multiple de-excitations contribute to the variance spectrum (see Table II). This can be explained by the fact4that the applied quench and shift manipulations pump more energy into the final system in the double-well studiescomparing with that in the harmonic case, see the appendix for more information. More pumped energy means thata larger number of the excited states are populated increasing thereby the probability of the de-excitations seen, e.g.,in the evolution of the integral two-body observables.The situation is different for the local quantity. In the quench scenario the local density ρ ( x = 1) does not provideinformation about the many-body excitations. As one can see in the right lower panel of Fig. 6 the mean-field GPexcitation g is the main and only excitation available. Since we measure ρ ( x = 1) close to the density minimum, higherorder excitations do not contribute significantly. However, taking the density at regions away from the minimum,the g excitation starts to contribute to the local dynamics. Concluding, the local density can be used for accessingmany-body excitations but a proper choice of the local position is very crucial. 3. Deep Double-Well BdG LR-MCTDHB(2) "shift" (GP) "shift" (MB) "quench" (GP) "quench" (MB)label ∆ E ( I ) ∆ E ( I ) ∆ E ( I ) ∆ E ( I ) ∆ E ( I ) ∆ E ( I )0 u a ) - - - g (cid:48) - 0.064 - - - 0.066 (0.36 b ) g b ,0.93 c ) 2.092 (0.33 b ,0.89 c ) u a ,0.66 c ) 2.095 (0.93 a ,0.70 c ) - - u (cid:48) - 2.158 (0.08) - 2.163 (0.02 a ,0.02 c ) - - g (cid:48) - 2.161 (0.09) - - - 2.161 (0.12 b ) g c ) - 3.914 (0.07 c ) 3.934 (0.07 b ,0.10 c ) u a ,0.23 c ) 3.943 (0.03 a ,0.15 c ) - - u (cid:48) - 4.029 (0.05) - 4.032 (0.02 a ,0.07 c ) - - g (cid:48) - 4.039 - - - - g (cid:48) ↔ g - 0.069 - - - g ↔ g b ) u ↔ u c ) 1.849 (0.02 c ) - - g ↔ g (cid:48) - 2.028 - - - 2.032 (0.02 b ) a observed in (cid:104) x (cid:105) b observed in Var( x ) c observed in ρ ( x = 1) TABLE III. Excited states in the deep double-well: The structure of the table is as in Table I. Finally let us study the dynamics in the deep double-well potential depicted in the left panel of Fig. 1. For thischoice of trapping potential the ground state is almost fully fragmented ( ≈ ) and, therefore, not described withinthe standard GP mean-field theory. We recall that the BdG approach is capable of predicting the linear responsesonly form a fully condensed state, so its applicability in the present case of the deep double-well and fragmentedground state is questionable, nevertheless, we provide the LR-GP result for the sake of completeness. Indeed, themean-field and many-body approaches already vary significantly in the description of the energy of the first excitedstate u . Table III shows that the corresponding excitation energy ∆ E within static LR-MCTDHB(2) is almost twotimes the energy of the lowest-in-energy excitation obtained by the mean-field BdG theory.In Fig. 7 we depict the discrete FT-results obtained from the evolutions of the local and integral observables afterapplying the quench and shift scenarios to the ground state of the deep double-well. From Table III it becomes clearthat the overall situation is similar to the shallow double-well. Most of the dynamic excitations can be identified withdistinct LR states. The main difference, however, is that excited states form bands, which result in more discretespectra.By analyzing the evolution of the (cid:104) x (cid:105) observable in the shift scenario we see that the u mean-field-like excitationcontributes most to the dynamics as expected from the corresponding LR intensities. This is true for both the mean-field and many-body descriptions. The predicted many-body excitations u (cid:48) and u (cid:48) have quite small LR intensitiesand, therefore, do not contribute to the dynamics. In the evolution of the local density ρ ( x = 1) the higher ordermean-field-like excitation u is more pronounced. Still, in the shift scenario the integral and local observables give5 ω i n t e n s it y [ a r b . un it s ] shift quench ρ (x=1) 1 2 3 4frequency ω i n t e n s it y [ a r b . un it s ] shift quench ρ (x=1) 0 0.5 frequency ω i n t e n s it y [ a r b . un it s ] shift quench Deep double-well trap: The description is as for Fig. 5. us the same picture of the excitations. It is important to note that the lowest-in-energy excitation u is not visiblealthough predicted by the LR intensities. A possible explanation could be in a small energy difference between thisand the ground state. One has to monitor the evolution of the system for a longer propagation time in order to reacha higher resolution and detect contributions from this transition.The right panels of Fig. 7 show the excitation spectra obtained in the quench scenario of the deep double-well.The overall situation is very similar to that observed in the shallow double-well. Namely, the evolution of thevariance Var( x ) contains contributions from the excitations to the first ( g , g (cid:48) ) and second ( g ) bands. However, theidentification of the most intense peak at almost zero energy is not trivial, because energetically it can be attributedto the many-body excitation g (cid:48) as well as to the inter-band de-excitation g (cid:48) ↔ g . According to the LR resultspresented in the Table III, the appearance of the g (cid:48) excitation is rather unfavorable, because its LR intensity is6negligible. Also unlikely is that the de-excitation g (cid:48) ↔ g constitutes the peak because its intensity in the comparableshallow double-well regime is small. However, by monitoring the many-body evolution for a longer propagation timewith the subsequent DFT analysis one should be able to split this lowest-in-energy peak into the contributing statesand give further insight into the identification problem. The local observable ρ ( x = 1) , similarly to the situationdescribed in the shallow double-well, can access the mean-field-like g excitations to the second band, but not thelow-lying many-body excitations.The dynamic studies on the shallow and deep double-well potentials reveal one common feature of the local density ρ ( x = x ) . This local quantity can be quite informative for detection of the many-body excitations in the shift scenariobut not in the quench scenario. Hence, the ungerade many-body excitations can in principle be measured by usinglocal observables, however, to account for gerade many-body ones, it is desirable to use the many-body operators,such as Var( x ) . IV. SUMMARY In this paper we have investigated excitations and de-excitations of 1D BECs in the crossover from a fully condensedto a fully fragmented regimes. The primary physical goal was to identify the nature and possible classes of the low-lying excitations and to compare and contrast predictions of the mean-field based (GP) and many-body (MCTDHB)theories. On the methodological side we have employed two different techniques to access the excited states: A staticapproach based on the linear response theories and a dynamic approach which utilizes the wave-packet propagation.In the first section to compute the excitation spectra we have applied the static linear response mean-field BdG ≡ LR-GP and many-body LR-MCTDHB theories. By comparing and contrasting the obtained results we have found thatthe many-body theory reproduces the excitations predicted at the mean-field level and, in addition, introduces newclass of excitations which have pure many-body origin. In particular, the low-lying many-body excitations havebeen observed in depleted and fragmented BECs confined in double-well potentials and, surprisingly, in condensedultracold bosons trapped in harmonic traps. In order to quantify the probability of these states to be excited byapplied perturbations we have computed their linear response amplitudes. For interacting systems the many-bodyexcitations have non-vanishing response amplitudes implying their relevance for quantum dynamics. The comparablevalues of the response amplitudes obtained for the many-body and mean-field excited states imply that the many-bodytreatment is unavoidable for a proper description of these systems. For further characterization of the excited stateswe have used their symmetry and nodal structure.In the dynamic studies we have proposed two simple protocols that address excitations manifesting themselves intemporal oscillations of the many-body wave-function and observables. We have investigated the dipole-like ungeradeoscillations by shifting the origin of the confining potential and the breathing-like gerade excitations by quenchingthe frequency of the parabolic part of the trap. The MCTDHB wave-function computed at each propagation time-step has been used to calculate the corresponding local ρ ( x = x ) and integral (cid:104) x (cid:105) , Var( x ) observables. By applyinga subsequent Fourier transformation to each evolving quantity we have extracted the excitation energies and theirintensities and compared them with the corresponding static linear-response predictions. The main methodologicalconclusion is that all the dynamic excitations obtained within the studied protocols are predicted by and contained inthe static LR-MCTDHB theory. Moreover, the intensities of the FT-extracted excited states and the correspondingstatic LR predictions share a general tendency: Excitations that have large response amplitudes in the LR predictionstend to be easily excited by the dynamic scenarios. This holds true for mean-field and many-body excitations.We have demonstrated that in the shift scenario the one-body observable (cid:104) x (cid:105) can be used to access the mean-field-likeand many-body excitations. The excitation probabilities of the mean-field-like states are, however, higher than thatof the many-body states. The two-body observable Var( x ) used in the quench scenario is more informative quantitybecause along with the mean-field-like and many-body excitations it contains information about some de-excitations.We have shown that these de-excitation energies can be obtained within the static LR-MCTDHB approach by takingdifferences between appropriate excited states. The intensities of some pure many-body states are found to be largeand comparable with mean-field-like states, implying that they can be directly detected. We would like to mentionthat the local density taken at different regions of the trap can also be used to access excited states. In particular,we have shown that the local densities taken at the regions of a comparatively low particle density, i.e., far from thepotential minima can be used to access some higher order excitations which have not been seen in the global pictureas given by (cid:104) x (cid:105) and Var( x ) . However, these predictions are found to be very sensitive to the selected position and thedynamical scenario used.Finally, in the present dynamic studies the quench and shift parameters of the Hamiltonian have been chosensuch that the pumped energy was less than a few percents. Nonetheless, we have observed small de-excitationscontributing to the evolution of the local and global observables indicating the beyond the linear regime was reached.As a consequence, in real experiments implementing simple dynamical scenarios it can be difficult to keep the system7in the linear regime and, therefore, a detailed knowledge of both excitations and de-excitations is crucial for a properinterpretation of the experimental data.Alternatively, one can think about more sophisticated dynamical scenarios and control protocols where a majortask would be to populate a desired excited state exclusively. This kind of population control required for quantumcomputing, see e.g. [42], can, in principle, be reached by merging the optimal control algorithms with appropriatereal-space many-body theories [30–33, 36]. Acknowledgments We are indebted to O.E. Alon and S. Klaiman for fruitful scientific discussions. Computation time on the bwGRiD[43], and HybriLIT [44] clusters as well as financial support by the DFG are greatly acknowledged. Appendix: Linear Regime In this appendix we specify, quantify and determine the regimes of the dynamic changes of the Hamiltonian. Inthe present studies all the dynamic manipulations on the system have been implemented by a sudden transition froman initial Hamiltonian H i to the final one H f involving, thereby, finite shifts and quenches. This observation shouldbe contrasted with the fact that the derivations of the LR methods are based on infinitesimal perturbations of thesystem. Moreover, the term linear in this respect specifies that only small first-order (linear) corrections induced bythe applied perturbations to the many-body wave-function are taken into account and considered.Hence we can formulate the problem as follow: To what extend in a given dynamic scenario can we assure a linearityas required in the context of the LR approaches?In order to answer this question let us study the energy which is pumped into a system by the applied dynamicprotocols. Let E GS be the ground state energy of the final system described by the H f and E tot be the total energy ofthe system after the applied dynamical scenario, i.e., after the sudden transition from H i → H f . The energy pumpedinto the system is then given by E pump = E tot − E GS . This amount of energy is redistributed between the excitedstates. Figure 8 shows the pumped energy with respect to the ground state energy when the potential, given byEq. (7), is shifted (left) or quenched (right).Obviously, we have observed the non-linear behavior of the energy as a function of the presented control parametersin both cases. So, for the studied systems the linear dependency of the pumped energy on the control parameteris unreachable in principle and, therefore, it cannot be considered as a criterion of the linear behavior. Instead, wepropose to rely on the magnitude of the pumped energy. For our studies we choose shifts or quenches that pumpa maximal amount of of the ground state energy into the systems. This choice is practical, because it offers acompromise between two extremes. On the one hand the magnitude of the pumped energy is large enough to beobserved experimentally, while on the other hand non-linear effects like de-excitations are small but visible. Appendix: Computational Details The derivation of the LR theory atop the MCTDHB wave-function is reported in Refs. [27–29]. The final result forthe resulting LR-MCTDHB theory boils down to the diagonalization of the non-hermitian LR matrix, i.e., takes onthe form of the eigenvalue equation: L u k v k C ku C kv = ω k u k v k C ku C kv . (A.1)The linear-response matrix L of the many-boson MCTDHB wave-function Ψ is more involved than the commonly-employed Bogoliubov-de Gennes (BdG) linear-response matrix. Physically, the response amplitudes of all modes, u k and v k , and of all expansion coefficients, C ku and C kv , combine to give the many-body excitation spectrum ω k = E k − E GS . Here E GS and E k are the energies of the MCTDHB ground state and excited states respectively. Wehave successfully managed to explicitly construct L and obtained the many-body excitation spectrum for bosonsinteracting by contact [27] and general [29] interparticle interaction potentials.8 E p u m p / E G S shiftwidth 0.3 0.5 0.7quench in a FIG. 8. Pumped energy: Energy pumped to the system by the proposed dynamic protocols which are based on the shiftand quench of the trap potential, see Eq. 7 and its discussion. In the left panel, the pumped energy is plotted against theapplied shift x shift . The right panel shows the energy induced to the system by quenching the trap frequency a . Vertical linesand arrows denote the shift and quench parameters used in the paper. Red lines indicate results for the harmonic potential,green lines for the shallow double-well and blue lines for the deep double-well. Calculations are done at MCTDHB( M = 4 )level for N = 10 particles and Λ = λ ( N − 1) = 1 . To quantify the intensity of the response we compute the response weights: I k ≡ γ k = (cid:104) u k | f + ( ρ ) / | φ (cid:105) + (cid:104) v k | f − ( ρ ) / | φ , ∗ (cid:105) + (cid:16)(cid:82) d r φ , ∗ i f + φ j (cid:17) (cid:104) C ku | ˆ a † i a j | C (cid:105) + (cid:16)(cid:82) d r φ , ∗ i f − , ∗ φ j (cid:17) (cid:104) C kv | (cid:16) ˆ a † j a i (cid:17) ∗ | C , ∗ (cid:105) . (A.2)Here all quantities with 0-superscript relate to the static MCTDHB solution and f ± define the driving amplitudes ofthe applied perturbations, see Ref. [27] for details.To identify the symmetry and nodal structure of the excited state we have computed and analyzed the real oscillatoryparts of the corresponding response densities ∆ ρ = Re[∆ ρ ko ( r ) + ∆ ρ kc ( r )] . The contributing orbital and CI-terms aregiven as following ∆ ρ ko ( r ) = M (cid:88) i,j =1 ρ ij φ i ( r ) (cid:8) ˜ u kj ( r ) + ˜ v kj ( r ) (cid:9) , (A.3)where ˜ u kj = ( ρ ) − / ji u ki , and ˜ v kj = ( ρ ) − / ji v ki . Further, we defined ˆ ρ = (cid:80) ij ˆ a † i ˆ a j φ , ∗ i ( r ) φ j ( r ) . ∆ ρ kc ( r ) = (cid:104) C | ˆ ρ | C ku (cid:105) + (cid:104) C k, ∗ v | ˆ ρ | C (cid:105) . (A.4)In our time dependent studies for the temporal evolution of the quantity f ( t ) given in time-slices f (0) , f (1) , f (2) , .., f ( N − we compute its Discrete Fourier Transform (DFT) which is the equivalent of the continuous Fourier Transform fora time-evolving signal given only at N instants (time-points) separated by equidistant sample times ∆ t = τ (i.e. a9finite sequence of data): F ( iω ) = (cid:90) ( N − τ f ( t ) e − iωt dt → N − (cid:88) k =0 f ( k ) e − ikωτ , (A.5)or, more precisely with ω = 0 , πNτ × , πNτ × , · · · , πNτ × ( N − F ( n ) = N − (cid:88) k =0 f ( k ) e − ik πnN . (A.6)We have normalized the obtained F ( n ) amplitudes such that their sum would be equal to the unity.In all our numerical simulations with the MCTDHB-Laboratory package [24] we have use the Sin-DVR, Exp-DVRor FFT-grid in a box ( − 10 : 10) with N g = 512 grid points. All the time propagations have been done till T = 500 . [1] F. Dalfovo, S. Giorgini, L. P. Pitaevskii, and S. Stringari, Rev. Mod. Phys. , 463 (1999).[2] A. J. Leggett, Rev. Mod. Phys. , 307 (2001).[3] I. Bloch, J. Dalibard, and W. Zwerger, Rev. Mod. Phys. , 885 (2008).[4] C. Pethick and H. Smith, Bose-Einstein Condensation in Dilute Gases , 2nd ed. (Cambridge University Press, New York,2008).[5] M. H. Anderson, J. R. Ensher, M. R. Matthews, C. E. Wieman, and E. A. Cornell, Science , 198 (1995).[6] K. B. Davis, M. O. Mewes, M. R. Andrews, N. J. van Druten, D. S. Durfee, D. M. Kurn, and W. Ketterle, Phys. Rev.Lett. , 3969 (1995).[7] R. Folman, P. Krüger, J. Schmiedmayer, J. Denschlag, and C. Henkel, Adv. At. Mol. Opt. Phys. , 263 (2002).[8] Quantum Gases: Finite Temperature and Non-Equilibrium Dynamics (Vol. 1 Cold Atoms Series), edited by N. P. Proukakis,S. A. Gardiner, M. J. Davis, and M. H. Szymanska (Imperial College Press, London, 2013).[9] R. Ozeri, N. Katz, J. Steinhauer, and N. Davidson, Rev. Mod. Phys. , 187 (2005).[10] I. Shammass, S. Rinott, A. Berkovitz, R. Schley, and J. Steinhauer, Phys. Rev. Lett. , 195301 (2012).[11] O. Penrose and L. Onsager, Phys. Rev. , 576 (1956).[12] P. Nozières and D. Saint James, J. Phys. France , 1133 (1982).[13] K. Xu, Y. Liu, D. E. Miller, J. K. Chin, W. Setiawan, and W. Ketterle, Phys. Rev. Lett. , 180405 (2006).[14] C. Orzel, A. K. Tuchman, M. L. Fenselau, M. Yasuda, and M. A. Kasevich, Science , 2386 (2001).[15] S. Hofferberth, I. Lesanovsky, B. Fischer, J. Verdu and J. Schmiedmayer, Nature Physics , 710 (2006).[16] M. Greiner, O. Mandel, T. Esslinger, T. W. Hänsch, and I. Bloch, Nature (London) , 39 (2002).[17] C. Menotti, J. R. Anglin, J. I. Cirac, and P. Zoller, Phys. Rev. A , 023601 (2001).[18] R. W. Spekkens and J. E. Sipe, Phys. Rev. A , 3868 (1999).[19] A. I. Streltsov, O. E. Alon, and L. S. Cederbaum, Phys. Rev. A , 063626 (2006).[20] E. J. Mueller, T. -L. Ho, M. Ueda, and G. Baym, Phys. Rev. A , 033612 (2006).[21] M.-K. Kang and U. R. Fischer, Phys. Rev. Lett. , 140404 (2014).[22] A. I. Streltsov, O. E. Alon, and L. S. Cederbaum, Phys. Rev. Lett. , 030402 (2007).[23] O. E. Alon, A. I. Streltsov, and L. S. Cederbaum, Phys. Rev. A , 23 (1947), in english.[26] J. Grond, A. I. Streltsov, L. S. Cederbaum, and O. E. Alon, Phys. Rev. A , 063607 (2012).[27] J. Grond, A. I. Streltsov, A. U. J. Lode, K. Sakmann, L. S. Cederbaum, and O. E. Alon, Phys. Rev. A , 023606 (2013).[28] O. E. Alon, A. I. Streltsov, and L. S. Cederbaum, J. Chem. Phys. , 034108 (2014).[29] O. E. Alon, J. Phys.: Conf. Ser. , 012039 (2015).[30] J. Grond, J. Schmiedmayer, and U. Hohenester, Phys. Rev. A , 021603(R) (2009).[31] J. Grond, G. von Winckel, J. Schmiedmayer, and U. Hohenester, Phys. Rev. A , 053625 (2009).[32] R. Bücker, T. Berrada, S. van Frank, J.-F. Schaff, T. Schumm, J. Schmiedmayer, G. Jäger, J. Grond, U. Hohenester, J.Phys. B: At. Mol. Opt. Phys. , 104012 (2013).[33] G. Jäger, T. Berrada, J. Schmiedmayer, T. Schumm, and U. Hohenester, Phys. Rev. A , 053632 (2015).[34] P. Doria, T. Calarco, and S. Montangero, Phys. Rev. Lett. , 190501 (2011).[35] T. Caneva, T. Calarco, R. Fazio, G. E. Santoro, and S. Montangero, Phys. Rev. A , 012312 (2011).[36] I. Brouzos, A. I. Streltsov, A. Negretti, R. S. Said, T. Caneva, S. Montangero, and T. Calarco, Phys. Rev.A , 062110(2015).[37] C. Chin, R. Grimm, P. Julienne, and E. Tiesinga, Rev. Mod. Phys. , 1225 (2010).[38] M. Edwards, P. A. Ruprecht, K. Burnett, R. J. Dodd, and C. W. Clark, Phys. Rev. Lett. , 1671 (1996). [39] C. R. McDonald, G. Orlando, J. W. Abraham, D. Hochstuhl, M. Bonitz, and T. Brabec, Phys. Rev. Lett. , 256801(2013).[40] S. Klaiman and O. E. Alon, Phys. Rev. A , 063613 (2015).[41] S. Klaiman, A. I. Streltsov, and O. E. Alon, Phys. Rev. A , 023605 (2016).[42] M. G. Bason, M. Viteau, N. Malossi, P. Huillery, E. Arimondo, D. Ciampini, R. Fazio, V. Giovannetti, R. Mannella, andO. Morsch, Nat. Phys.8