Capture Dynamics of Ultracold Atoms in the Presence of an Impurity Ion

We explore the quantum dynamics of a one-dimensional trapped ultracold ensemble of bosonic atoms triggered by the sudden creation of a single ion. The numerical simulations are performed by means of the ab initio multiconfiguration time-dependent Hartree method for bosons which takes into account all correlations. The dynamics is analyzed via a cluster expansion approach, adapted to bosonic systems of fixed particle number, which provides a comprehensive understanding of the occurring many-body processes. After a transient during which the atomic ensemble separates into fractions which are unbound and bound with respect to the ion, we observe an oscillation in the atomic density which we attribute to the additional length and energy scale induced by the attractive long-range atom-ion interaction. This oscillation is shown to be the main source of spatial coherence and population transfer between the bound and the unbound atomic fraction. Moreover, the dynamics exhibits collapse and revival behavior caused by the dynamical build-up of two-particle correlations demonstrating that a beyond mean-field description is indispensable.


I. INTRODUCTION
The theoretical description of degenerate atomic quantum gases has advanced significantly in the past two decades. It relies particularly on the separation of length scales of the external trapping potential, the interparticle distances and the range of atomic interactions. In the ultracold regime, this allows to approximate the interaction by a contact potential (or simply pseudopotential) [1], which has proven tremendously powerful when applied to bosonic gases. The latter holds especially for the mean-field approximation leading to the well-known Gross-Pitaevsikii (GP) equation [2,3] which describes the condensed state of a weakly interacting bosonic gas yielding a variety of phenomena such as collective excitations, solitons, and vortices [4,5]. For state-of-theart methods, such as the density-matrix renormalization group [6] and the multiconfiguration time-dependent Hartree method for bosons (MCTDHB) [7] as well as its multi-layer extension (ML-MCTDHB) [8,9], the pseudopotential represents an essential simplification for the quantum dynamical description of many-boson systems in order to investigate the physics beyond the mean-field approximation. Even though the short-range interaction has proven to lead to exotic states of quantum matter like supersolidity [10] or crystalline phases [11] and therefore represents an important case, not all ultracold atomic systems exhibit these short-range interactions. An example for the latter are chromium atoms which possess long-range magnetic dipolar interactions [12,13].
Hybrid atom-ion systems represent a specific class of systems with a new scale of interactions due to the inter- * jschurer@physnet.uni-hamburg.de play between the charge of an ion and the induced dipole moment of a neutral atom [14]. They have recently become available experimentally [15][16][17][18] and attracted increasing interest. Most of the current experiments, however, are based on the Paul trap scheme, whose drawback is the so-called micromotion which so far prevents from reaching the ultracold regime as shown in theoretical classical and quantum analyses [19,20]. Corresponding studies showed that a large ion-atom mass ratio might help circumventing this limitation [19][20][21], although it is not yet clear whether the s-wave regime can be reached. Still, a recent detailed study has shown that specifically for the atom-ion pair Li-Yb + the ultracold regime can be reached experimentally [22]. Given these findings, it is desirable to extend the theoretical understanding of ultracold neutral quantum many-body systems by the presence of ions. This attractive interaction induces, additional to the confinement, a further length and energy scale and it is therefore expected to lead to intriguing effects such as the formation of molecular ions [23] and ioninduced density bubbles in the atomic cloud [24]. Apart from this fundamental point of view, these hybrid systems show versatile applications in quantum information processing as, for example, the controlled creation of entanglement [21,25], the realization of quantum gates [26], and the simulations of solid-state systems [27].
In a previous study, we have investigated in detail the ground-state properties of an atom-ion hybrid system, consisting of a single static ion in the center of a bosonic atomic cloud. The dependence of relevant observables on the atom number and the interaction strength has been analyzed and we showed that the presence of an ion strongly affects them [28]. For weakly interacting atoms, we found that the ion impedes the transition to the Thomas-Fermi regime while for strong atom-atom interaction it modifies the fragmentation behavior depend-ing on the atom number parity. In the present work, we explore the impact of the second length and energy scale generated by the atom-ion interaction onto the dynamics of the atomic cloud. In particular, we envision the scenario in which we create a single ionic impurity in an ensemble of N interacting atomic bosons. If the atomic cloud, initially prepared in the ground state of a harmonic trap, and the ion are suddenly brought into contact, this process resembles, to some extent, the sudden ionization of a single impurity atom within the atomic cloud. Here, however, we still neglect the ionic motion which is justified in case the ion is tightly trapped. We will see in the following that the second length and energy scale induces a coherent oscillation between states bound in the atom-ion potential and states of the harmonic trap. This oscillation occurs in addition to the usual harmonic excitation and reveals a collapse and revival behavior caused by the dynamical build-up of correlations.
This work is organized as follows. In Sec. II, we define the setup and explain the model of our hybrid atom-ion system. In addition, we introduce a cluster-expansion scheme which enables us to distinguish the single-particle dominated physics from the processes induced by the build-up of (quantum) correlations. In Sec. III, we explore the dynamical evolution of the hybrid system and identify the most important excited modes by means of our cluster-expansion approach. Sec. IV contains a brief convergence analysis which shows the necessity to use an advanced method like MCTDHB. Finally, we summarize our findings in Sec. V together with the conclusions and an outlook on future investigations.

II. MODEL AND THEORETICAL APPROACH
In the following, we introduce our model and define the process to initiate the dynamics. Further, we briefly outline the MCTDHB used for the simulations, define important quantities for the analysis, and introduce our cluster expansion approach which we exploit in order to analyze the many-body wave function.

A. Model of the System
We consider N interacting bosonic atoms in a onedimensional (1D) harmonic trap at zero temperature initially prepared in the ground state of the system which we compute by imaginary time propagation [29] of an initial guess wave function. The short-range intra-atomic interaction is modeled by a contact pseudopotential. Into this atomic cloud, we immerse a single trapped impurity atom which does not interact with the other atoms. This could be achieved by tuning the inter-atomic interaction to zero by exploiting a Feshbach resonance [30]. At time t 0 = 0 this impurity atom is ionized by a laser pulse such that a single ion is created within the atomic cloud. The ionizing laser pulse is assumed to be far detuned from any possible resonance of the atomic cloud. Moreover, we assume pulse durations such that the ionization process takes place on much faster time scales than any possible response of the atomic cloud. This allows us to treat the ionization as an effectively instantaneous process. If the atomic impurity is trapped in a sufficiently deep and tight trap, the proposed scenario can be indeed achieved experimentally, as recently reported for optical trapping of a single ion [31]. Thus, at time t 0 = 0, the ion is assumed to be statically trapped at z I = 0 in the atomic cloud. The motional excitations following this ionization process result in a rich dynamics which we analyze in Sec. III.
The interaction between the ion at z I and an atom at z A behaves in 1D at large distances as −αe 2 /(2(z A −z I ) 4 ) up to a minimal cutoff distance R 1D [32], where α is the polarizability of the atoms and e the elementary charge. For numerical many-body simulations with MCTDHB it is more convenient, however, to define a model potential as [28]: Here z denotes the relative coordinate z A −z I . The model parameters v 0 ,γ, and ω are determined by the shortrange quantum defect parameters. The above potential asymptotically approach the −1/z 4 behavior at large distances, whereas at short distances it has a barrier which is designed such that the quantum defect theory results for the atom-ion scattering [33] can be reproduced. Additionally to the harmonic confinement, this interaction introduces a second length R * = αe 2 m/ 2 and energy scale E * = 2 /(2mR * 2 ) to the system, where m is the mass of a single neutral atom. In summary, the system can be described by the following many-body Hamiltonian (in E * and R * units) with the harmonic trap of frequency ω 0 and characteristic length l = /(mω 0 )/R * , and the intra-atomic contact interaction strength g. The ionic part is switched on at time t 0 by a step-function θ(t−t 0 ). Thereby, we term the Hamiltonian of a single boson asĥ i . In the following, we denote the stationary eigenenergies ofĥ i (t > t 0 ) for fixed t with ǫ j and the corresponding single-particle eigen- Hereafter, for the sake of numerical convenience, we set l = 0.5R * , which would corresponds to ω 0 = 2π · 3.3 kHz and l = 188 nm for 87 Rb atoms. We will see that this choice does not affect qualitatively the observed dynamics. Furthermore, we choose a weak intra-atomic interaction strength g = 2E * R * and small atomic ensembles The effective potential (gray shaded area in background) consisting of the harmonic trap and the atom-ion interaction potential for the 87 Rb atom. Further, we show the single-particle energy levels ǫj indicated by straight horizontal lines. In addition, the energetically lowest single-particle state φ 0 1 (z) (dark gray area) and the trap state φ 0 5 (z) (light gray area) are sketched with arbitrary but equal scaling. Possible types of dynamics occurring in the system are indicated by arrows.
consisting of N = 2 up to N = 10 neutral atoms. Besides, we fix the model parameters to ω = 80(R * ) −4 , v 0 = 3ω and γ = 4 √ 10ω. We refer here to Ref. [28] for a detailed discussion of the chosen parameters as well as for the experimental conditions needed for the quasi-1D regime. Our above choice leads to two bound states for the atoms in the atom-ion potential which are localized on both sides of the ion but vanish at z I . Even though we neglect the motion of the ion, we term the two states below E = 0 bound states while we refer to the remaining states as trap states (with E > ω 0 /2). In Fig. 1, we show the total potential for the atoms together with the energies ǫ j as well as the lowest single-particle state φ 0 1 (z), bound in the atom-ion potential, and the state φ 0 5 (z). The arrows indicate the possible processes that can occur: within the trap (green arrow), in the atomion potential (white arrow), and between those two scales (oblique magenta arrows).

B. Theoretical Approach
We explore the quantum dynamics of the many-body system described by the wave function |Ψ by means of the numerically exact ab initio method MCTDHB. Its main idea is that by using m time-dependent variationally optimized single-particle basis functions, the number of basis functions can be kept rather small. More precisely, the many-body wave function |Ψ for N bosons is expanded in bosonic number states |n(t) in order to take into account the indistinguishability of the bosons. Note that in a number state |n(t) , each boson occupies one of the m time-dependent singleparticle functions (SPFs) |Ψ j (t) and that the vector n = (n 1 , · · · , n m ) contains the occupation numbers n j of every SPF. Besides, the sum in Eq. (3) goes over all possible n with j n j = N , which is denoted by the symbol n|N . With this ansatz for the many-body wave function, the temporal evolution of the wave function |Ψ(t) is obtained by means of the Dirac-Frenkel variational principle [34,35] which guarantees a variational optimal many-body solution. We would like to emphasize that not only the coefficients A n (t) but also the SPFs |Ψ j (t) are adapted in time to the many-body dynamics in order to allow for the largest possible overlap between the ansatz (3) and the true many-body wave function. We refer for a detailed description of the method to Refs. [7][8][9].
The analysis of the full many-body wave function |Ψ is generally a complicated task due to the underlying high dimensionality resulting from the N degrees of freedom. In order to analyze the time-dependent many-body wave function in detail, we inspect two different quantities: The one-and two-particle reduced density matrices which allow for the investigation of spatial coherence and correlations [36], and clusters enabling us to analyze coherence and correlations in terms of any single-particle basis. We transfer and adapt the notion of clusters from Ref. [37] to bosonic systems of few particles.

Density Matrices
The reduced one-and two-particle density matrices are defined via the expectation values of the field operatorŝ Ψ † (x) andΨ † (x) as ρ 1 (x, y, t) = Ψ † (x, t)Ψ † (y, t) and ρ 2 (x, y, y ′ , x ′ , t) = Ψ † (x, t)Ψ † (y, t)Ψ † (y ′ , t)Ψ † (x ′ , t) , respectively. Their spectral decomposition can be written in terms of the natural populations λ j (t) and the natural orbitals Φ j (x, t) and the natural populations γ j (t) and natural geminals Φ j (x, x ′ , t) respectively. This natural representation of the reduced density matrices has several advantages. For example, the natural populations can be used to identify the degree of fragmentation of the system [38] and to judge the convergence of our numerical simulations [39]. Further, the natural orbitals build always a very suitable singleparticle basis set for the description of the dynamical system, hence in many cases only a few basis functions are needed to represent the (many-body) wave function in this basis. Despite these advantages, the analysis and a deep physical understanding of the quantum dynamics in this basis is very difficult, since the natural orbitals are time-dependent eigenfunctions of the density matrix in its spatial representation. Given this, we use instead the so-called clusters for the detailed analysis of the manybody dynamics. We define in the following the notion of clusters in the framework of the cluster-expansion approach.

Definition of Clusters and Correlations
The cluster-expansion approach is a very powerful technique to describe the quantum dynamics of interacting many-body systems because it allows for a systematic truncation of the Bogolyubov-Born-Green-Kirkwood-Yvon (BBGKY) hierarchy [40]. By expanding M -particle expectation values into cumulants or correlated clusters, a consistent theory up to the single-, two-, or even M -particle level can be developed [41]. In this spirit, we shall use the cluster-expansion approach in order to separate the many-body dynamics, obtained by means of MCTDHB, into single-and two-particle contributions which will enable us to get more physical insight. In particular, it will be useful for the identification and classification of the most relevant excitations of the system that are created during the dynamical evolution (see Sec. III).
To this end, we will first briefly review the clusterexpansion approach, which will be helpful for a better understanding of our modifications to the traditional approach. Indeed, we shall adapt the traditional clusterexpansion to bosonic systems with a fixed particle number. This constraint is particularly important for atomic ensembles of only a few particles.
To begin with, let us consider an arbitrary orthonormal time-independent single-particle basis φ j (x) and the annihilation and creation operatorsâ † j andâ † j for the corresponding single-particle state, respectively. Then, the expectation value of any observable can be expressed in terms of the M -particle expectation values or M -particle clusters with the index sets j = {j 1 , ..., j M } and k = {k 1 , ..., k M } and M ≤ N . These are nothing else but the matrix elements of the M -particle reduced density matrix in an arbitrary basis and can be, most conveniently for MCT-DHB, obtained from the spectral representation of the M -particle reduced density matrix, as shown for M = 1 and M = 2 in App. A . Now, the single-particle properties of the system are given by the single-particle clusters â † iâ † j , also called singlets. We distinguish between occupations f j = â † jâ † j , which describe the population of the state φ j (x), and coherences p ij = â † iâ † j (i = j) which can be understood as a transition amplitude between the i-th and the j-th state.
Starting from the singlets, the cluster-expansion is recursively build-up based on the consistent factorization of an M -particle cluster into independent particles (singlets), correlated pairs, correlated three particle clusters, up to correlated M -particle clusters [37,42], that is: Here the terms [ M ] S represent the single-particle contributions, while the terms ∆ M C contain the correlated part of the M -particle cluster. Note that we omitted the indices for brevity such that all cluster products, denoted in square brackets, include a sum over all unique permutations. For instance, for the two-particle Given this, the correlated part of the twoparticle cluster is given by ∆ 2 C := 2 − [ 2 ] S . Now the correlated clusters ∆ M C with M > 2 can be determined recursively which completes the formulation of the cluster expansion.
At this point a remark is in order: For bosonic systems with a fixed number of particles (i.e., in a number-conserving theory), even in a GP type meanfield state (i.e., with only one single-particle orbital), the term ∆ M C is non-zero because of the bosonic symmetry. This implies that the systematic truncation of the BBGKY hierarchy, for which the correlated parts of any M -particle cluster beyond a certain size (M > M T ) have to be neglected, cannot be performed in this case. Note that this issue does not arise neither for fermions nor for bosons with particle number fluctuations. Thus, in order to circumvent this problem, we shall introduce a slightly different definition of the correlated parts ∆ M C . Our strategy will be to define the correlated parts of any Mparticle cluster in such a way that all the (M −1)-particle contributions are indeed removed. As a result of such a strategy, for example, the correlated parts automatically vanish in any mean-field state.
To this end, let us note that an M -particle cluster contains all the information about the M ′ -particle clusters with M ′ < M . This can be easily seen by using the recursive relation between the density matrices for M ≥ 2, which leads to In order to identify the correlated part of an M -particle cluster, we decompose the M -particle cluster into two parts: one consisting of all contributions from clusters with M ′ < M , denoted by M j,k <M , and one which contains only the M -particle contributions, that is, There are several ways one could perform such a decomposition. For instance, in the traditional cluster-expansion approach, as we have discussed above, one would choose the following definition for the single-particle part of the two-particle cluster: Here, however, we shall define the term M j,k <M by requiring that it has to fulfill the condition This expression is assumed to hold for any many-body quantum state and implies that q ∆ M {j,q},{q,k} = 0 [43]. Besides, if we consider, for instance, the case M = 2, then we see that the right-hand side of Eq. (13) accounts only for single-particle contributions of two-particle clusters.
In order to find a definition for M j,k <M such that Eq. (13) is fulfilled, let us first investigate, as an example, the case where |Ψ is a general mean-field state, which is defined as a single permanent. More precisely, given some single-particle orbitals χ j (x) with the associated creation operatorsα † j , a mean-field state is defined as a single permanent like |MF k = N iα † ki |vac with |vac being the vacuum and k = {k 1 , ..., k N } (for the commonly known GP state k i = k 1 ∀i). For such a state the two-particle clusters can be written as with the bosonic correlations [44] given by Here n j is the number of bosons in the single-particle state χ j (x) defined via n j = i δ ki,j (δ i,j is the Kronecker-Delta). It is easy to verify that for such a mean-field state Eq.
With this choice, the term ∆ 2 {k,q},{q ′ ,k ′ } is zero in a mean-field state per definition. In this spirit, the terms ∆ M j,k can be understood as the M -particle correlations.
Now for a general many-body quantum state we replace in Eq. (15) the mean-field occupations n j and the mean-field basis functions χ j (x) by the natural populations λ j and the natural orbitals Φ j (x), respectively, yielding the following analogue expression: Note that in general λ j is not an integer. Hence, we define the single-particle contributions of a two-particle cluster as whereas the two-particle correlations, also called doublets, are defined as Here some considerations are in order. Our choice for the non-correlated part of the two-particle cluster [Eq. (17)] indeed fulfills the condition (13) which justifies the above replacement of the mean-field occupations and basis functions by the natural populations and natural orbitals, respectively. Although we focused here on the twoparticle clusters, we note that expressions like Eq. (18) can be obtained for any M -particle clusters, too. However, for the present study the singlets and doublets are sufficient to understand the dynamics of the system. Further, we would like to stress that the decomposition into singlets and M -particle correlations does not correspond to a separation into mean-field and beyond mean-field contributions. But even if the condition (13) only separates the clusters into singlets, doublets, three-particle correlations, etc., it enables us to identify genuine correlations of any many-body quantum state which are beyond mean-field. Indeed, in the limit of a mean-field state Eq. (16) boils down to Eq. (15), and therefore all correlated parts ∆ 2 {k,q},{q ′ ,k ′ } = 0 vanish. This would not be possible with the traditional cluster-expansion approach.
Finally, we would like to highlight that one could instead search for the best mean-field state [45] and then separate the M -particle clusters into mean-field part and contributions beyond that. It turns out, however, that such a choice does not satisfy Eq. (13). This shows that a mean-field approximation of a real many-body state does not necessarily result in the exact single-particle properties of the system.

Singlet Dynamics
Even though we are able to derive the time evolution of every M -particle cluster from our MCTDHB solutions, it is worth to investigate the equations of motion of the clusters since they make it possible to understand the coupling between the singlets themselves and between singlets and doublets. The dynamics of the M -particle clusters can be derived from the Heisenberg equation of motion and the above outlined definitions. For the analysis of the dynamics in Sec. III, we focus onto the equations of motion of the singlets. One can easily show that the time evolution of the coherences is given by while the one of the populations is governed by Here we used the definition of Eq. (18) for the twoparticle cluster that appears in the corresponding Heisenberg equation, because of the atom-atom interaction and the asterisk for the complex conjugation. Besides, in the above equations, we have introduced the renormalized single-particle energiesǫ i = ǫ i + Σ ii , the singlet couplings Σ ij = 2 kq V kijq â † kâ † q , the coupling to the bosonic correlations and to the doublets and the interaction matrix elements Note that the equations of motion for the singlets form a system of non-linear differential equations, since the mean-field couplings are also defined by the singlets themselves. Furthermore, they are coupled to the doublets (via Γ ij ) which can be interpreted as a source term (or inhomogeneity) for the singlets. On the other hand, the dynamics of the doublets depends on the singlets as well as on the three-particle correlations ∆ 3 , which renders the coupling time-dependent. We would like to emphasize that the coupling to ∆ 3 is a manifestation of the BBGKY-hierarchy. For the truncation of the hierarchy at the singlet level, one has to ensure that the contribution of the doublets to the singlet dynamics remains very small during the time interval of interest, which would imply that one can set Γ ij = 0. In addition, we remind here again that in order to obtain a closed singlet theory, the coupling to the bosonic correlations has to be included which is the price to pay for our choice of factorization. In the following, we show how a consistent and closed singlet theory can be accomplished. At first, by means of Eqs. (13) and (17), we obtain an exact relation among the bosonic correlations and the singlets of the system: With this expression for the coupling to the bosonic correlations, the equations of motion of the singlets can be completely decoupled from higher than single-particle clusters such that a consistent and closed singlet theory is obtained. In the following, we will see the power of such a singlet theory in the analysis of the complicated many-body dynamics, especially in combination with the MCTDHB method.

III. DYNAMICAL EVOLUTION
We investigate now the dynamics induced by the instantaneous creation of an ion in an atomic cloud. First, we analyze the one-body density (matrix) and the components of the energy as well as the many-body excitation spectrum. The observations are then discussed and analyzed in detail in terms of our singlet-doublet theory developed above.

A. Observations
In Fig. 2, we show exemplarily the temporal evolution of the one-body density ρ(z, t) = ρ(z, z, t) of the Time evolution of the one-particle density ρ(z, t) and the components of the total energy per particle for a system with N = 2. The green, cyan, and magenta lines represent the trapping, kinetic, and ionic energy per particle, respectively.
atomic cloud for g = 2E * R * and N = 2. In order to give a feeling of the actual time scales involved, we note that for 87 Rb we have /E * ≈ 0.39 ms, while for 7 Li it corresponds to /E * ≈ 0.001 ms. For short times, the suddenly created ion captures very quickly most of the atomic cloud in its bound states. The remaining atomic density fraction is emitted as a beam into the outer region of the harmonic trap. Consequently, this fraction is decelerated (t < 0.2 /E * ) and back reflected (0.2 /E * < t < 0.4 /E * ) by the harmonic confinement. Subsequently, this sequence repeats with approximatively constant frequency. Furthermore, a second faster oscillation in the density fraction captured within the ionic potential is visible (see the holes in the density plot at z ≈ ±R * /2). Additionally, we show in Fig. 2 the components of the energy per particle. The trapping energy (green line) perfectly oscillates with a single frequency which coincides with the oscillation frequency of the outer density fraction. In contrast, the ionic energy (magenta line), which represents the expectation value of the ionic potential (1), oscillates with a higher frequency matching the inner density oscillation. On top of this oscillation, we observe a short pulse when the outer fraction "crashes" into the inner density part. Note that these events are not visible in the trapping energy since the harmonic trap energy is negligible in the vicinity of the trap center. The kinetic energy (cyan line) can be understood as the negativ sum of trapping and ionic energy, since the interaction energy (not shown) is comparably small and the total energy is conserved during the dynamics. Hereafter, we term the oscillation of the outer fraction of the atomic cloud harmonic and the inner fraction ionic oscillation. Note that the above observations are qualitatively independent of the atom number N such that Fig. 2 is representative also for larger N . After multiple oscillation periods (see Fig. 3), we observe that the ionic energy (magenta line), and therefore the ionic oscillation, exhibits a clear collapse and revival behavior on this long time scale. On the other hand, the trapping energy (green line), and thus the harmonic oscillation, becomes only slightly damped. While for N = 2 these two aspects of the long-time behavior can be barely seen, it becomes strongly pronounced for larger particle numbers.
The collapse and revival behavior can also be observed in the long time dynamics of the atomic density. Figure  4 shows the disappearance of the ionic oscillation during the time interval [1.5, 3.5] /E * and its recurrence around t ≈ 4.0 /E * . This effect seems to be directly connected to the loss and regain of spatial coherence between the inner and outer density fraction which can be observed in the snapshots of the one-particle density matrix at times t = 0.54 /E * (left panel), t = 2.20 /E * (middle panel), and t = 4.61 /E * (right panel) in Fig. 4. Below, we will understand the relation between the ionic oscillation and the spatial coherence in detail through the singletdoublet analysis.
Let us now discuss the frequencies that are involved in the dynamics. To this aim, we investigate the fidelity defined by the overlap of the many-body wavefunction at time t 0 = 0 with the one at time t [46]: Since this fidelity F (t) can be understood as the expectation value of the time-evolution operator, thus an Nbody operator, its Fourier transform contains information of all involved excited eigenstates of the interacting N atom system. Due to the discrete nature of the spectrum and because in our numerical simulations the propagation time T is finite, it turns out to be more efficient for the computation of the Fourier transform to use the compressed sensing (CS) method [47,48]. With compressed sensing, one can indeed obtain a resolution in frequency space better than ∆ω = 2π T [49]. To this end, we have used the matlab package "SPGL1" for compressed sensing from Ref. [50]. The algorithms used in this package can be found in Refs. [51,52]. In Fig. 5, we show the Fourier transform F (ω) of the fidelity (red continuous line). Several prominent resonances become apparent. We observe one dominant mode at a frequency ω ≈ 34 E * / . This corresponds to the ionic oscillation frequency. In contrast, the harmonic frequency can not directly be found in the spectrum. Furthermore, we see that modes with high frequency are excited and seem to have an equidistant spacing. In addition to this, there is a low energy mode around ω ≈ 15 E * / , whose origin we will explain in the subsequent sections.

B. Singlet Dynamics
Let us start the analysis of the many-body spectrum in terms of the singlets. Therefore, we choose from here on the non-interacting single-particle functions φ 0 j (x) as the basis φ j (x) in which the clusters are expressed. In general, the identification of the modes corresponding to the observed resonances is a very complicated task. Nevertheless, we were able to identify the most important modes by linearizing the equation of motion for the singlets [see Eqs. (19) and (20), and App. B]. The obtained energies together with the "exact" spectrum belonging to the Fourier transform of the fidelity [see Eq. (26)] are shown in Fig. 5 (blue circles). One observes a very good agreement of the obtained peak positions with the Fourier spectrum of F (t) (red continuous line). On the other hand, the relative heights of the peaks can be only obtained approximatively by our singlet theory (see App. B) showing qualitative agreement with the Fourier spectrum. Further, we can identify which singlet has the dominant contribution to a certain resonance (see App. B for details). This dominant singlet is indicated in Fig.  5 at the corresponding peak.
We find that the most dominant frequency corresponds to p 13 = â † 1â † 3 excitations. These are coherences between the number states |N, 0, 0, ... and |N −1, 0, 1, 0, ... (in the basis of the non-interacting single-particle states φ 0 j (z)) which oscillates with the frequency ω I = (ǫ 3 − ǫ 1 )/ = 38.7 E * / for g = 0. Thus, the ionic oscillation corresponds to an oscillation between a bound state and  F (t). The blue circles mark the singleparticle energies (see text). We indicate the corresponding dominant singlet p kq at each peak by the label pi,j for the sake of better readability. The dashed vertical lines illustrate the non-interacting limit, thus correspond to g = 0. Note that p1,5 is hard to identify due to its small amplitude. the first trap state (see also the oblique magenta arrows in Fig. 1). Hence, it connects the inner part of the atomic cloud with the outer fraction establishing spatial coherence between them as we have seen in Fig. 4. Moreover, this mode induces, as we will see in the following, population transfer between the inner and the outer atomic fraction. On the other hand, the harmonic oscillation can not be attributed to a single resonance, because no resonance shows up at its frequency. Nevertheless, we can understand its origin. We note that singlets â † 1â † i with i = 5, 7, 9, ... to very high i are excited and therefore present in F (ω). They correspond to oscillations between the numberstate |N, 0, 0, ... and the states with |N − 1, 0, ..., 0, 1, 0, ... , where a single particle is excited to the ith state and thus they oscillate, for g = 0, with frequencies ω i = (ǫ i − ǫ 1 )/ with i = 5, 7, 9, .. (note that excitations into states with i even are symmetry forbidden). These frequencies are shown in Figs. 5 as dashed lines. Since the difference between two neighboring resonances is approximatively constant (ω i+2 −ω i ) ≈ 2ω 0 due to the equidistant spacing of the energy levels in the harmonic trap, all of these singlets are in phase again after a time T ≈ π/ω 0 ≈ 0.4 /E * . Although this approximation does not work perfectly at all times, because of the presence of the ionic potential, the harmonic oscillation is visible for the entire simulation time due to the contribution of high energy modes which are only slightly affected by the presence of the ion. Nevertheless, the damping of the harmonic oscillation, visible in the trap energy (see also Fig. 3, green line), can be attributed to the induced slight dephasing.
In the inset of Fig. 5, the resonance positions in dependence of the particle number N are shown. We observe that the peaks which can be associated to a coherence p 1j shift to smaller energies for growing N . Our singlet theory is perfectly able to capture this behavior. Therefore, the shift of resonances can be understood as a meanfield phenomenon which effectively renormalizes the singlet resonances. Importantly, we numerically find that the property of equidistant spacing between the singleparticle modes seems to be nearly untouched by the interaction and, as a consequence, the harmonic oscillation is essentially unaffected.
The singlet theory explains most of the modes found in the Fourier spectrum F (ω). The low energy resonance at ω ≈ 15.5 E * / , however, does not appear in the singlet spectrum. Moreover, a closer inspection of the inset of Fig. 5 shows that this mode has the tendency of an increasing energy as the number of bosons increases, which is the opposite behavior of the other resonances. We show in the following that this resonance can be attributed to correlations which are dynamically build-up during the temporal evolution. To do so we have first to understand the coupling between the singlets and the doublets.
Towards that end, we proceed further with the analysis of the dynamical evolution of the singlets. In Fig.  6, the coherences p 13 , p 1T = k>3 p 1k , and the occupations f 1 , f 2 , f T = k>3 f k are shown. As stated be-forehand, the coherence p 13 oscillates with the frequency of the ionic oscillation. From the time evolution of the ionic energy (see Fig. 3), one expects in addition the collapse and revival behavior. Here we see that for N = 2 (lower left panel), the absolute value of p 13 is nearly constant, while collapse and revival become clearly visible for N = 10 (lower right panel). The coherence between the first bound and the trap states, p 1T , shows the aforementioned rephasing behavior in the distinct peaks with T ≈ π/ω 0 ≈ 0.4 /E * periodicity. The damping is also visible here for N = 10, but, importantly, the peaks are still very pronounced enabling the harmonic oscillation to be unaffected. Turning our attention to the dynamics of the occupations (Fig. 6, upper panels), we can identify that the population of the bound states, thus the ion population f I = f 1 + f 2 , is about 80 − 90%, and therefore only about 10% − 20% of the particles are emitted into the trap by the "ionization process". Furthermore, a population transfer between the two bound states occurs (upper left panel). For larger N (upper right panel), also transfer to and from the trap population (green line) becomes visible with a rate equal to the frequency of the ionic oscillation.
An inspection of Eq. (20) reveals that the population transfer between the bound state φ 0 1 and the trap states is mediated by coherences, primarily by p 13 , because it has the highest contribution in the Fourier spectrum F (ω) (see also Fig. 5), explaining its oscillation with the frequency of the ionic oscillation. In contrast, the dynamical evolution of f 2 has to be induced by the correlations Γ 22 , since all coherences in question vanish. Hence, the population transfer from the state φ 0 1 to the state φ 0 2 would not take place in a mean-field scenario, thus it is a clear signature of a genuine many-body effect. We also note that the frequency of this population transfer is ω ≈ 8.2 E * / for N = 2 which matches very well to the position of the additional low energy peak seen in the inset of Fig. 5. Further, the larger N is, the faster the transfer process between the first bound and the trapped states happens which fits to the shift of the low energy resonance to larger energies visible in Fig. 5. We therefore show, in the next section, which doublets play an important role for this oscillation and to which process/mode it corresponds. Further, we will explain the origin of the collapse and revival of the ionic oscillation.

C. Doublet Dynamics
In order to understand the observed collapse and revival phenomena and the incoherent population transfer, we need to investigate the dynamical evolution of the doublets. For the sake of clarity, we classify the most important doublets (see definition in Sec. II B 2) into incoherent doublets g kkkk , g kqkq , and coherent doublets g kkqq , g kqkk . This separation is justified by the fact that the incoherent doublets are always real-valued quantities, and therefore do not contribute to the dynamics of the singlet occupations [see Eq. (20)]. Moreover, the doublets have the properties g kqk ′ q ′ = g * k ′ q ′ kq and g kqk ′ q ′ = g kqq ′ k ′ = g qkk ′ q ′ = g qkq ′ k ′ . Therefore, we can identify, noting that only correlations involving states up to k = 3 are contributing to the system dynamics, the most important classes of correlations: incoherent {g 1111 , g 2222 , g 3333 , g 1212 , g 1313 } and coherent {g 1122 , g 1133 , g 2233 , g 1211 , g 2122 , g 1311 , g 3133 , g 2322 , g 3233 }.
In Fig. 7, we show for N = 2 (left panels) and N = 10 (right panels) only those doublets of the incoherent (upper panels) and the coherent (lower panels) doublets which are considerably occupied. At first, we observe that, even though we start in an essentially uncorrelated state, very quickly a tremendous amount of correlations is created in the course of the temporal evolution. With the help of these correlations, we are now able to answer the remaining open questions concerning the dynamical evolution of the system.
In order to understand the population transfer between the states φ 0 1 and φ 0 2 , we have to identify the correlations g qkk ′ q ′ contributing to Γ 11 and Γ 22 . The coherent doublets g 1122 (g 2211 ) is the only correlation in question which can be understood by inspecting Eq. (22) and Fig. 7. Comparing Figs. 6 and 7 (left panels), we see that they are build up and decay with the frequency of the population transfer between the state φ 0 1 to the state φ 0 2 . Since g 1122 and g 2211 which drive Γ 11 and Γ 22 , respectively, are by π out of phase (g 1122 = g * 2211 ), they induce the population transfer between those two bound states. Thus, we can understand these coherent doublets as mediators for the population transfer. Further, we can now identify the low energy mode in the spectrum of Fig. 5 as a coherent oscillation between the numberstates |N, 0, ... and |N − 2, 2, 0, ... which corresponds to a two-particle excitation. In the noninteracting case, this mode would therefore oscillate with a frequency ω D = 2(ǫ 2 − ǫ 1 )/ ≈ 7.4 E * / indicated by a dashed dotted line in the inset of Fig. 5. Here we can understand the following: First, this oscillation can not be seen in single-particle quantities like the one-body reduced density matrix. Second, this is also the reason why we could not explain the low energy peak by our singlet theory. Third, any mean-field approach would not be able to capture the associated dynamics. For example, in a GP theory the population transfer from |N, 0, ... to |N − 2, 2, 0, ... would be not possible, because in a general GP state ( k c kâ † k ) N |vac contributions from odd states are symmetry forbidden as long as only GP states with parity symmetry are considered.
Finally, we would like to discuss possible reasons for the collapse and revival behavior of the ionic oscillation occurring for larger N . The impact of the doublets onto the singlet p 13 is contained in Γ 13 . Even if many doublets contribute to Γ 13 , we can identify the responsible doublet by inspecting Fig. 7 (right panels). We observe that exactly at the times when the coherent doublets g 1133 are present, the ionic oscillation is strongly suppressed whereas when g 1133 nearly vanishes the ionic oscillation reappears (compare Figs. 3, 4, and 6). Therefore, the doublet g 1133 is the main source of loss of the coherence p 13 . Consequently, this doublet is responsible for the loss of spatial coherence between the inner and the outer density fraction during the dynamics, as it can be seen in Fig. 4. Nevertheless, the build-up of g 1133 does not only act as a source of damping for the singlet coherences, since the coherences p 13 recur when the value of g 1133 is reduced again. Further, we note that g 1133 (g 3311 ) is also present in Γ 11 (Γ 33 ) such that also the collapse and revival in the population transfer can be attributed to this coherent doublet.

IV. CONVERGENCE ANALYSIS
In the following, we briefly discuss the quality of our data and to which extend advanced tools as MCTDHB in the description of the dynamical behavior of such quantum systems are indeed necessary. To this aim, we inspect the natural populations [see Eq. (4)] which provide an assessment of the degree of fragmentation of the system [38]. Further, one can use them to systematically judge the convergence of our result to the exact solution of the many-body quantum system [39]. In Fig.  8, the evolution of the natural populations exemplarily for N = 10 particles is shown. We see that even if the initial state is nearly condensed, thus describable in the GP framework, depletion from the condensate state is created very quickly during the time evolution. Therefore, a GP description would be inaccurate and a multiorbital description becomes unavoidable. Here we are using m = 6 single-particle orbitals. The contribution of the lowest natural orbital (smallest natural population) stays below 1% such that only natural orbitals with even smaller contribution are expected to be not taken into account by truncating the Hilbert space. One could further speculate that in such cases a multi-orbital mean-field description might be sufficient. However, this is not the case and we stress that the non-steady behavior of the natural populations implies the necessity to go beyond a mean-field description [53]. Therefore, the utilization of MCTDHB is essential to capture the full correlated quantum dynamics of the system.

V. CONCLUSIONS AND OUTLOOK
We have investigated the dynamics of a trapped cloud of N interacting bosons after the sudden creation of a static ion. Thereby, most of the atoms are quickly captured in the bound states of the atom-ion potential, while the remaining atomic fraction is emitted into the trap. We were able to understand the subsequent dynamics and the many-body spectrum in detail by means of an underlying singlet-doublet theory.
Our atom-ion hybrid system exhibits, apart from the harmonic oscillation, that is, a density oscillation within the external harmonic confinement, an additional density oscillation which we were able to attribute to a coherent oscillation between a bound state and a trap state. This coherent oscillation results in the spatial coherence between the inner and outer density fraction. It also allows for population transfer between bound and trap states such that the fraction of atoms in the bound states becomes time-dependent. In contrast to the harmonic oscillation, this ionic oscillation shows a particle number dependence. Furthermore, we showed that, even though the atoms are weakly interacting, strong correlations are build-up during the temporal evolution. These correlations lead to population transfer between the bound states and show up as a strong doublet resonance in the many-body spectrum. Besides this, the correlations induce a periodic suppression of coherence between the inner and the outer fraction which gives rise to collapse and revival of the ionic oscillation on longer time scales.
The impact of the ion onto the bosonic cloud dynamics can be summarized as follows: Apart from the accumulation of atoms on both sides of the ion while depleting the cloud at the ion position, the ion induces a coherent oscillation between the bound states and the trap states. The oscillation essentially arises because of the additional scales present in the system. Indeed, this can be easily seen in the non-interacting case (g = 0) and taking the limit of small trapping frequency ω 0 → 0. In this case, the frequency of the ionic oscillation converges to the non-interacting single-particle energy ǫ 1 of the lower bound state (see Sec. II A for definition) which, thus, sets effectively a lower bound to the frequency of the ionic oscillation. Hence, although we have chosen a tight trap, which is numerically convenient to resolve the dynamics in a finite system, the ionic oscillation would be present with a frequency on the same order of magnitude for a weak confinement, too.
The present work can be viewed as a further step towards the simulation of the dynamics of the hybrid atomion system in the ultracold regime. The fact that quasi one-dimensional quantum Bose gases are routinely realized in several laboratories and given the recent advances in optical trapping of ions, puts the experimental realization of the hybrid system considered here within reach. Furthermore, since the multilayer extension of our method (ML-MCTDHB) is especially designed for the simulation of mixtures, we plan to investigate the impact of the ionic motion onto the ground state and the dynamical evolution of the bosonic ensemble in the nearer future. This will be done by treating the ion quantum mechanically as well. Due to the attractive inter-particle interaction, these future studies can reveal the dynamical formation of molecular ions as predicted in Ref. [23]. Furthermore, it would be of interest to study the energy transfer between ionic and atomic degrees of freedom in order to understand, for example, sympathetic cooling of the ion in the atomic cloud in more detail.  In summary, we see that the linearized singlet equations provide us with the singlet resonances ν n and the singlet modes S R n . We can derive these from the MCT-DHB solution in the following way: Via Eq. (A5), we can extract the dynamics of the singlets which we can use to derive the matrix M , and therefore the solution S(t). By diagonalizing M , we obtain the singlet energies and the dominant singlet in the eigenmodes (largest entry in the vector) shown in Fig. 5. Further, we can extract the oscillator strength |c n | 2 by Eq. (B7) which we can use as a measure for the relative height of the peaks in the many-body Fourier spectrum of F (t) just by scaling them to the global maximum of the spectrum. Note that, Eq. (B7) should be evaluated at small t in order to allow the linear approximation.