Spectral Structure and Many-Body Dynamics of Ultracold Bosons in a Double-Well

Abstract

We examine the spectral structure and many-body dynamics of two and three repulsively interacting bosons trapped in a one-dimensional double-well, for variable barrier height, inter-particle interaction strength, and initial conditions. By exact diagonalization of the many-particle Hamiltonian, we specifically explore the dynamical behavior of the particles launched either at the single-particle ground state or saddle-point energy, in a time-independent potential. We complement these results by a characterization of the cross-over from diabatic to quasi-adiabatic evolution under finite-time switching of the potential barrier, via the associated time evolution of a single particle’s von Neumann entropy. This is achieved with the help of the multiconfigurational time-dependent Hartree method for indistinguishable particles (MCTDH-X)—which also allows us to extrapolate our results for increasing particle numbers.

1. Introduction

The detailed microscopic understanding of interacting many-particle quantum dynamics in state-of-the-art experiments with ultracold atoms [1,2,3,4,5,6,7,8,9,10] in well-characterized potential landscapes remains a challenging task for theory: While a large arsenal of advanced numerical techniques has been developed over the past two decades to efficiently simulate interacting many-particle dynamics [11,12,13,14,15], all of them must ultimately surrender when confronted with truly complex dynamics, i.e., under conditions where a generic initial state fully explores, on sufficiently long time scales, an exponentially large Hilbert space in the number of particles and/or degrees of freedom. By the very meaning of complexity, even the most efficient numerical methods can only be expected to yield reliable results when the dynamics can be restricted to finite sub-spaces of the exponentially large Hilbert spaces—either by reducing the time window over which the evolution is followed, or by choosing physical situations which a priori confine the many-particle state. This has been long understood in the light–matter interaction of atoms and molecules [16], as well as in quantum chaos [17], and meets revived interest given the experimental progress in the control of cold matter [18].

While it is, therefore, clear that the only promising route for an efficient characterization of large and complex quantum systems can be through effective descriptions—such as offered, e.g., by the theory of open quantum systems [19,20,21,22], modern semiclassics [23], or random matrix theory [24,25]—there is an intermediate range of system sizes where efficient numerical methods can (a) be gauged against each other, to benchmark their quantitative reliability, without any a priori restriction on the explored portion of Hilbert space, and (b) contribute to gauge effective theories against (numerically) exact solutions [26,27,28], at spectral densities where quantum granular effects induce possibly sizeable deviations [29] from effective theory predictions (which always rely on some level of coarse graining). In our view, it is this intermediate system sizes where efficient methods of numerical simulation develop their full potential, since they can inspire and ease the development, e.g., of powerful statistical methods and paradigms (such as scaling properties [18,26,30])—which then enable robust predictions in the realm of fully unfolding complexity.

In the present paper, we contribute to this line of research by exploring the spectral and dynamical properties of a few bosonic particles loaded into a symmetric double-well potential, with static or switchable barrier. Prima facie, this is a well-known and text-book-like example, yet with a panoply of experimental realizations, and of paradigmatic relevance as an incarnation of Josephson dynamics [4,31,32,33,34,35,36,37,38,39,40,41,42,43,44,45,46,47] or as the elementary building block of quantum dynamics in lattice-like structures [48], and quickly defines a formidable numerical challenge if only one admits excitations far beyond the immediate vicinity of the ground state energy, and seeks to accurately monitor the long-term dynamics of two or more particles. We will see how the spectral structure of the single-particle problem is amended by adding a second, identical particle, and how finite-strength interactions restructure the many-particle spectrum and eigenstates, throughout the excitation spectrum up to the vicinity of the potential barrier.

Ultimately, we have two main goals: First, we want to shed light onto the nature of the tunneling processes of two, repulsively interacting bosons launched either at the ground state or saddle-point energy. Secondly, we want to study the transition between diabatic and adiabatic switching of the potential barrier for different particle numbers. Importantly, we will find that tunneling is described by a second-order two-particle process and not by a direct first-order two-particle process. We underpin our studies by information from the respective regimes in the spectrum. These scenarios come with very different challenges for the numerical treatment, because the evaluation of the time evolution generated by a time-independent ordinary differential equation in case of a static potential significantly differs from the one described by a time-dependent ordinary differential equation in case of a switchable barrier. To achieve our above two main goals and as a central aspect of our present contribution, we employ a variety of different numerical approaches which, by accounting for the complete spectral structure of the double-well (rather than the lowest-lying band alone), go far beyond widely used single-band Bose–Hubbard models.

Here and in the following, we use the term “many-body/particle”, albeit the systems we consider are composed of a relatively small number of particles. Please note that our considerations are from first principles and start from the many-body Hamiltonian. Moreover, it has been shown theoretically [17,49,50] and experimentally [1,2] that the physics of interacting few-body systems can very quickly approach the many-body limit.

The spectral information thus generated allows us to decipher characteristic features of the many-particle dynamics, for distinct choices of the initial condition, and over a wide range of interaction strengths, for static as well as for diabatically or (quasi-)adiabatically ramped potential barriers. Finally, we illustrate, through an analysis of the von Neumann entropy of the (reduced) single-particle density matrix, how such transition from diabatic to (quasi-) adiabatic switching controls the effectively explored sub-volume of Hilbert space, and how robust coarse grained features of the resulting “phase diagram” emerge as the particle number is increased from two to ten. The latter case can only be treated with the help of the MCTDH-X [28,51,52] method which has been verified against exact [27,28] and experimental [53] results and is reviewed in Ref. [15]. Here, we push MCTDH-X to its limits in monitoring long-term dynamics of rather moderate, mesoscopic particle numbers, in the presence of strong, switching-induced excitations (“quenches”).

The paper is organized as follows: The theoretical framework, including a brief description of the numerical methods, is presented in Section 2. Section 3 is devoted to the discussion of the spectral and eigenstate structure of the problem at hand. First, Section 3.1 discusses how the energy spectrum depends on both the tunneling barrier height and the inter-particle interaction strength, for two and three particles. Next, in Section 3.2, we study few-body correlations encoded in the few-body eigenstates. This prepares our analysis of the dynamics in Section 4. In Section 4.1, we investigate the dynamics of two particles in a static double-well potential, initially prepared in two different states: A superposition of low-lying states, and a superposition of excited states with energies close to the saddle point. Finally, we consider the scenario of a time-dependent potential in Section 4.2: With the atoms initially prepared in the ground state of a harmonic trap, a central barrier is ramped up, and the thereby induced dynamics can be tuned from diabatic to (quasi-) adiabatic by appropriate control of the ramping time. Our results are summarized in Section 5.

2. Hamiltonian and Methods

2.1. Hamiltonian of Trapped Interacting Bosons

The Hamiltonian of N spinless, ultracold atoms with repulsive contact interaction and confined to a one-dimensional double-well potential reads in atomic units

$$\begin{array}{ccc}\hfill H& =& \sum _i^N\left(-\frac{1}{2}\frac{d^2}{d{x}_i^2}+V(x_i,t)\right)+\frac{\lambda }{2}\sum _{i\ne j}\delta (x_i-x_j),\hfill \end{array}$$

(1)

where

$$\begin{array}{c}\hfill V(x_i,t)=\frac{x_i^2}{2}+A(t)e^{-x_i^2/2}\end{array}$$

(2)

allows for a non-trivial time dependence of the potential barrier, through the time dependence of $A(t)$, $x_i$ denotes the position of the ith particle, and the repulsive interaction strength $\lambda >0$ is determined by the s-wave scattering length and the transverse confinement [54].

The minimum of $V(x_i,t)$ is located at $x=0$ if $A(t)<1$ (single-well), or at $x=\pm \sqrt{2ln(A(t))}$ if $A(t)\ge 1$ (double-well). Both static and time-dependent barriers will be considered. In the static case, the central barrier amplitude is constant, $A(t)=A_{\max }$, whereas in the time-dependent scenario, the amplitude is ramped up linearly according to

$$A(t)=A_{\max }\times \left\{\begin{array}{c}t/T_{\mathrm{ramp}},t<T_{\mathrm{ramp}},\hfill \\ 1,t\ge T_{\mathrm{ramp}}.\hfill \end{array}\right.$$

(3)

2.2. Numerical Methods and Observables

The spectral and dynamical properties of the Hamiltonian (1) are numerically investigated by using three approaches: the Fourier Grid Hamiltonian (FGH), the Bose–Hubbard (BH) representation of a continuous potential, and the multiconfigurational time-dependent Hartree method for indistinguishable particles (MCTDH-X); see Appendix A, Appendix B and Appendix C, respectively.

Each of these is suited for a specific task. We use FGH and BH which, ultimately, rely on different basis set representations of the Hamiltonian, to infer the spectrum of $N\le 2$ and $N=3$ interacting bosons, by direct diagonalization. FGH is also useful for the investigation of the quenched dynamics when a harmonic potential with $A_{\max }=0$ at $t=0$ is suddenly transformed into a static double-well with fixed barrier $A_{\max }=const.$ at $t>0$ (in other words, $T_{\mathrm{ramp}}\to 0$ in Equation (3)). For our study of the case of $N\le 10$ interacting bosons in a time-dependent double-well with $T_{\mathrm{ramp}}\ne 0$, we use the MCTDH-X method which enables accurate results for the dynamics, but cannot provide the complete spectral information as the BH/FGH methods. Since dynamical properties of interacting many-particle systems emerge, already at rather small particle numbers [17], the combination of all three approaches can be considered complementary.

FGH and BH yield the N-particle eigenenergies

$$\begin{array}{c}\hfill H|{\Psi }_n\rangle =E_n^{NP}(A_{max},\lambda )|{\Psi }_n\rangle ,\end{array}$$

(4)

with $|{\Psi }_n\rangle$ the N-particle eigenvector with quantum number n. All eigenstates are normalized to unity, throughout this paper. The quantity

$$\begin{array}{c}\hfill {\left|{\psi }_n(x_1,x_2,\dots ,x_N)\right|}^2=|\langle x_1,\dots ,x_n|{\Psi }_n\rangle {|}^2\end{array}$$

(5)

yields the associated probability density to find N bosons located at positions $x_1,x_2,\dots ,x_N$, respectively. Visualizations thereof reflect the correlations between the positions of the particles [37,38,40,41,43,44,53,55,56], which can be assessed, e.g., through their entanglement. A possible (though certainly non-exhaustive) quantifier of the non-separability of a general many-particle state $|\Psi (t)\rangle$ is given by the von Neumann entropy

$$S(t)=-\mathrm{Tr}\left[{\rho }^{1P}(t)ln{\rho }^{1P}(t)\right]$$

(6)

of the reduced single-particle density matrix [37,38,57,58,59], where ${\rho }^{1P}(x,x^{\prime },t)$ is defined as the trace over all degrees of freedom of all but one boson of the full density operator, i.e.,

$${\rho }^{1P}(t)={\mathrm{Tr}}_{2,\dots ,N}\left[|\Psi (t)\rangle \langle \Psi (t)|\right].$$

(7)

In particular, $S=0$ if the state is separable, while large values of S are a hallmark of a strongly entangled many-particle state [60,61,62,63,64]. We note that wide-spread mean-field approaches like the time-dependent Gross–Pitaevskii equation presuppose a separable many-body state with $S=0$; any $S\ne 0$ thus heralds the breakdown of such a mean-field description.

To characterize the dynamics of two bosons, we monitor the time evolution of the particles’ probabilities to reside both in the right (RR) or left (LL) well, or of each occupying one well (LR), given by [55]

$$\begin{array}{cc}\hfill P_{(LL)}(t)=& {\int }_{x_{\min }}^0\mathrm{d}{x}_1{\int }_{x_{\min }}^0\mathrm{d}{x}_2{|\psi (x_1,x_2;t)|}^2,\hfill \\ \hfill P_{(RR)}(t)=& {\int }_0^{x_{\max }}\mathrm{d}{x}_1{\int }_0^{x_{\max }}\mathrm{d}{x}_2{|\psi (x_1,x_2;t)|}^2,\hfill \\ \hfill P_{(LR)}(t)=& 2·{\int }_{x_{\min }}^0\mathrm{d}{x}_1{\int }_0^{x_{\max }}\mathrm{d}{x}_2{|\psi (x_1,x_2;t)|}^2,\hfill \end{array}$$

(8)

where we defined the three mutually distinct domains $(LL)=(x_1<0,x_2<0)$, $(RR)=(x_1\ge 0,x_2\ge 0)$, and $(LR)=(x_1<0,x_2\ge 0)\vee (x_1\ge 0,x_2<0)$. We also introduced the minimum ($x_{\min }$) and maximum ($x_{\max }$) values of the grid in configuration space employed in the numerical approaches. In addition, we evaluate the time-integrated probability current

$$J_{(RR\to LR)}(t)=2·{\int }_0^t\mathrm{d}{t}^{\prime }{\int }_0^{x_{\max }}\mathrm{d}{x}_2\mathrm{Im}{\left[{\psi }^{*}(x_1,x_2;t^{\prime })\frac{\partial }{\partial x_1}\psi (x_1,x_2;t^{\prime })\right]}_{x_1=0},$$

(9)

where the factor 2 accounts for the bosonic symmetry. $J_{(RR\to LR)}$ is derived [65] from the continuity equation and measures the probability flux within a time interval t from domain $(RR)$ to domain $(LR)$. This quantity is particularly important to distinguish first-order pairwise tunneling $(RR\to LL)$ from second-order pairwise tunneling $(RR\to LR\to LL)$. First-order pairwise tunneling was observed [65], e.g., for attractively interacting bosons in a double-well, where $J_{(RR\to LR)}(t)=0,\forall t,$ when the particles are initially prepared in one well.

3. Structure of Spectrum and Eigenstates

3.1. Few-Body Excitation Spectra

Since the dynamics of the system is ultimately encoded in its spectrum, we first discuss the parametric evolution of the eigenvalues (4) of $N=1,2$ and 3 bosons with both the central barrier height $A_{\max }$ and the interaction strength $\lambda$.

The single-particle spectrum is obtained by solving the time-independent Schrödinger equation

$$\begin{array}{c}\hfill \left(-\frac{1}{2}\frac{d^2}{d{x}^2}+\frac{x^2}{2}+A_{\max }{e}^{-x^2/2}\right){\psi }_n(x)=E_n^{1P}{\psi }_n(x).\end{array}$$

(10)

Figure 1 shows the evolution of the single-particle eigenenergies $E_n^{1P}$ as the central barrier height $A_{max}$ is continuously increased from a harmonic trap ($A_{max}=0$) to a deep double-well ($A_{max}=30$).

In the harmonic limit, the spectrum exhibits the well-known harmonic oscillator structure $E_n^{1P}(A_{max}=0)=n+1/2$. As the eigenenergies dive into the region below the barrier $A_{max}$ (indicated by the red diagonal in Figure 1a), the odd and even harmonic oscillator states become (nearly) degenerate. Sufficiently above $A_{\max }$, the energies are only weakly perturbed by the central barrier and we essentially recover the harmonic oscillator energy levels. In the limit $A_{\max }\to \infty$, the two wells decouple, leading to a fully degenerate harmonic oscillator spectrum.

From the structure of the single-particle spectrum, we can already anticipate that different dynamical behaviors can be expected for initial conditions with energies chosen below or above $A_{max}$, as will be elaborated upon, subsequently.

We now turn our attention to the spectrum of two particles obtained with the FGH method. The exact two-body spectrum is calculated by diagonalization of Equation (1) represented in the single-particle basis, as explained in Appendix A. Figure 2a shows that for $A_{max}=0$, we recover the well-known spectrum of two non-interacting bosons in a harmonic trap, i.e., $E_n^{2P}(A_{max}=0)=n+1$, with $n=n_1+n_2$ and degeneracy $g=n/2+1$ ($g=(n+1)/2$) for even (odd) $n\ge 0$. Here again, raising the central barrier gradually introduces a further degeneracy in the spectrum: The first three lowest-lying states become (nearly) degenerate when increasing $A_{max}$. This effect, also discussed in Ref. [38], is a direct consequence of the two-fold degeneracy of the single-particle ground state of the double-well, since all the eigenstates $|{\Psi }_0{\Psi }_0\rangle$, $|{\Psi }_0{\Psi }_1\rangle$ and $|{\Psi }_1{\Psi }_1\rangle$, with

$$|{\Psi }_n{\Psi }_m\rangle \equiv \frac{|{\Psi }_n\rangle \otimes |{\Psi }_m\rangle +|{\Psi }_m\rangle \otimes |{\Psi }_n\rangle }{\sqrt{2}\sqrt{1+\langle {\Psi }_n|{\Psi }_m\rangle }},$$

(11)

acquire the same energy value at large $A_{max}$ (see Equations (A9) and (A10)). For higher excitations, an analogous effect is observed: e.g., the energies of the states $|{\Psi }_2{\Psi }_0\rangle$, $|{\Psi }_2{\Psi }_1\rangle$, $|{\Psi }_3{\Psi }_0\rangle$ and $|{\Psi }_3{\Psi }_1\rangle$, respectively given by the sums of single-particle energies, $E_2^{1P}+E_0^{1P}$, $E_2^{1P}+E_1^{1P}$, $E_3^{1P}+E_0^{1P}$, and $E_3^{1P}+E_1^{1P}$, converge when increasing $A_{max}$, since $E_1^{1P}\simeq E_0^{1P}$ and $E_3^{1P}\simeq E_2^{1P}$. Therefore, the entire spectrum of two non-interacting bosons, plotted in Figure 2a, can be understood solely in terms of the single-particle spectrum. The emergence of the sequence of quasi-degenerate states is clearly observed below the separatrix $E^{2P}=2A_{max}$, plotted in red in Figure 2a.

Turning on the interaction changes the structure of the energy spectrum, as shown in Figure 2b. The calculation of the energy spectrum in the general case $A_{max}\ne 0$ requires a numerical treatment, whereas an analytical solution exists for the harmonic trap with $A_{max}=0$ and $N=2$ [59,66]. The most striking feature is the opening of an energy gap, clearly observed at large $A_{max}$: At the ground-state level, the three-fold degenerate states for $\lambda =0$ split into a unique ground state which remains unperturbed by the interaction, plus two (nearly) degenerate excited states which are affected by the non-vanishing interaction strength $\lambda \ne 0$. This behavior was already discussed in Ref. [38] for a polynomial double-well. Our present results show that this effect is also observed in the excitation spectrum below the separatrix $2A_{max}$. For instance, the first excited state manifold of the $\lambda =0$ limit (see Figure 2a, in the range $A_{max}\ge 10$), which is four-fold degenerate, splits (for $\lambda =1$, Figure 2b) into two (nearly) degenerate states unperturbed by the interaction, plus two (nearly) degenerate states slightly shifted by the interaction. The presence of these energy gaps in the spectrum will be essential for our understanding of the many-particle dynamics discussed in the next sections.

Consideration of a deep double-well, e.g., $A_{\max }=30$, allows for a better understanding of interaction-induced spectral features, as shown in Figure 3.

Indeed, for energies $E_n^{2P}\ll 2A_{\max }$, one can approximate the two wells by two decoupled harmonic traps with vanishing tunneling coupling. Flat energy levels correspond to the situation where the particles are almost completely localized in opposite wells and, consequently, do not interact. The remaining energy levels represent configurations where both particles occupy the same well. The spectral lines then approach the next higher-lying manifold at strong interaction, e.g., $\lambda \simeq 10$. In the limit of $\lambda \to \infty$, one recovers the Tonks-Girardeau (or fermionization) limit where these states become degenerate [38,43,44] with the second excited state manifold. Please note that by construction, this limit is out of reach for the single-band (or two-mode) approximation widely used in the literature. Figure 3 shows that the trend towards degeneracy between even and odd states with increasing $\lambda$ (fermionization process) is not restricted to the first spectral manifolds, but clearly manifests itself in the entire spectral range $E_n^{2P}\ll 2A_{\max }$.

The situation is (again) very different for three interacting particles [42]: Figure 4 shows the three-particle energy levels, for $A_{\max }=30$, as a function of the interaction strength U. All states are sensitive to the interaction and we observe two manifolds of states—states which exhibit interactions of two particles (full lines), and states which exhibit interactions of three particles (dashed lines). In contrast to the two-particle case, the ground state remains two-fold quasi-degenerate at large $\lambda$. Please note that the present three-particle results were obtained with the BH method (see Appendix B), since the Hamiltonian matrix is sparse in the BH representation, and therefore allows for computationally more efficient handling than the FGH method, for which the eigenenergies converge only slowly as a function of $N_{\mathrm{cut}}$ [56]. Furthermore, in the BH method $U\equiv \lambda {\sum }_i{|w_{0i}|}^4$, cf. Equation (A21), substitutes for $\lambda$ used in the FGH calculations.

3.2. Eigenstate Structure and Few-Body Correlations

Let us now inspect the associated many-particle eigenstates and the spatial correlations encoded into them, again as a function of both the central barrier height $A_{\max }$ and the interaction strength $\lambda$. The probability density, Equation (5), provides useful intuition. For two non-interacting bosons, the probability densities ${\left|{\psi }_n(x_1,x_2)\right|}^2$ are plotted in Figure 5, for energetically low- and high-lying eigenstates, as well as for different choices of the barrier height $A_{\max }$.

At low energies $(n=0,1,2)$, and with increasing barrier height $A_{\max }\to \infty$, ${\left|{\psi }_n(x_1=0,x_2)\right|}^2\to 0$ and ${\left|{\psi }_n(x_1,x_2=0)\right|}^2\to 0$. Consequently, the maxima of the probability density symmetrically split into the two or four corners of configuration space [37,38,67]. For $n=1$, the nodal line $x_1=-x_2$ originates from the superposition of even and odd (nearly) degenerate single-particle states. Please note that the associated eigenenergies are quasi-degenerate at $A_{\max }=30$: $E_{n=0,1,2}^{2P}\simeq 11.34$. At higher excitations, where the spectrum must progressively approach that of a harmonic oscillator (recall Figure 1b), the eigenstates exhibit a metamorphosis, sometimes even displaying a maximum at the saddle point, see, e.g., $n=76$, $A_{\max }=10$, and thus reminiscent of barrier states of the single-particle problem.

Interactions affect the spatial correlations in many ways, as shown in Figure 6 for $\lambda =1$: Comparison to Figure 5 shows that for $n=0-3$, the interaction slightly stretches the maxima of the eigenstates along the anti-diagonal $x_2=-x_1$ [38], and in some cases suppresses the amplitudes for double-occupancy of either site or that of delocalization over both sites. In a deep double-well, e.g., $A_{\max }=30$, the three-fold (nearly) degenerate non-interacting eigenstates $n=0-2$ of Figure 5 split into a unique ground state and two (nearly) degenerate eigenstates $n=1,2$. At higher excitations ($n=76$), we observe product states in the relative $\propto x_1-x_2$ and center-of-mass $\propto x_1+x_2$ coordinates (see Figure 6 for $A_{\max }=0$), and, therefore, also for these states correlated tunneling is expected, as opposed to the independent tunneling imprinted into the eigenstates in Figure 5. The impact of interactions on states in the vicinity of the separatrix is mainly highlighted by a suppression of the density maximum around $x_1=x_2=0$, see the result for $A_{\max }=10,n=76$ in Figure 6.

Next, let us have a closer look at the three-body probability density ${\left|{\psi }_n(x_1,x_2,x_3)\right|}^2$ of the ground state $(n=0)$ in a deep double-well, $A_{\max }=30$. Figure 7a,d show the three-body probability density (5) for non-interacting, $U=0$, and interacting, $U=1$, particles, respectively (see Equation (A21)). Since all particles occupy the same single-particle orbital $|{\psi }_0\rangle$, the non-interacting ground state covers all eight octants of configuration space in Figure 7a.

Like in the two boson case, the three-body wave function develops a nodal line along the main diagonal $x_1=x_2=x_3$ for non-vanishing $U>0$. At strong interaction, the maxima of the wave function are additionally shifted towards the corners of configuration space, along the diagonals $x_1=x_2=-x_3$, $x_1=-x_2=x_3$ and $-x_1=x_2=x_3$. Using a two-mode description, the ground state for sufficiently strong interactions is given by two particles at the same site and one on the opposite site. Therefore, the ground state, illustrated in Figure 7d, has no density in the areas associated with three particles at the same site ($x_1,x_2,x_3>0$ and $x_1,x_2,x_3<0$). Moreover, the two-mode description in the Fock basis $|n_L,n_R\rangle$ helps to understand the structure of the doubly degenerate ground state, since both states

$$\begin{array}{cc}\hfill |{\psi }_1\rangle & =|2,1\rangle ,\hfill \\ \hfill |{\psi }_2\rangle & =|1,2\rangle ,\hfill \end{array}$$

(12)

give rise to the same energy. The degenerate first and second excited states are then given by

$$\begin{array}{cc}\hfill |{\psi }_3\rangle & =|3,0\rangle ,\hfill \\ \hfill |{\psi }_4\rangle & =|0,3\rangle ,\hfill \end{array}$$

(13)

which are energetically even more sensitive to the interaction than the ground state doublet. Therefore, the four-fold degenerate ground state in the non-interacting case evolves into two doublets of states which further separate as a function of the interaction strength, as illustrated in the spectrum in Figure 4.

Finally, we inspect how the correlation information imprinted into the three-particle state is reduced when subsequently integrating out degrees of freedom. Averaging over one degree of freedom leads to the diagonal of the reduced two-body density matrix $|{\psi }_0(x_1,x_2){|}^2=\int \mathrm{d}{x}_3{|{\psi }_0(x_1,x_2,x_3)|}^2$, plotted for $U=0$ and for $U=1$ in Figure 7b,e, respectively. The impact of interaction becomes clearly visible by the reduction of the density along the diagonal $x_1=x_2$, tantamount of reduced correlations—as already observed in Figure 5 and Figure 6. Please note that in some contrast to the density of the two-particle state $n=0$, for $\lambda$ = 1 and $A_{\max }=30$ in Figure 6, the probability to detect two particles in the same well is not fully suppressed at interaction strength $U=1$.

Averaging over the second degree of freedom leads to the diagonal of the reduced one-body density matrices, $|{\psi }_0(x_1){|}^2=\int \mathrm{d}{x}_2\mathrm{d}{x}_3{|{\psi }_0(x_1,x_2,x_3)|}^2$, displayed in Figure 7c,f. The profile of $|{\psi }_0(x_1){|}^2$ for $U=0$, cf. Figure 7c, is exactly the same as the one obtained for the non-interacting two-particle case (red line), as expected. Only a small difference between the one-body densities $|{\psi }_0(x_1){|}^2$ associated with interacting and non-interacting (red line) bosons, respectively, is detectable, cf. Figure 7f (please note that the two-mode approximation (i.e., the double-well Bose–Hubbard model) is not sensitive to changes of the intra-well correlations—which here manifest themselves in the changed one-body density profile). This analysis therefore indicates that even if the interaction strongly affects the correlations, this information is not reflected by the one-body density profile.

4. Dynamics in the Double-Well

4.1. Static Potential: Two-Body Excited State Dynamics

Given the above phenomenology of spectra and eigenstates, we now explore how the tunneling dynamics of two interacting particles in a static double-well depends on the choice of the initial state. To this end, we consider a system initially prepared in a superposition of excited states, such that both particles are localized on the right-hand side of the double-well, at fixed barrier height $A_{\max }=10$. This localized state can be constructed by coherent superposition of (non-interacting) adjacent, even and odd one-body eigenstates:

$$|{\Psi }_n^{\mathrm{loc}}(t=0)\rangle =\frac{1}{2}\left(|{\Psi }_{2n+1}^{1P}\rangle +|{\Psi }_{2n}^{1P}\rangle \right)\otimes \left(|{\Psi }_{2n+1}^{1P}\rangle +|{\Psi }_{2n}^{1P}\rangle \right).$$

(14)

The dynamics is deduced from a spectral decomposition of the many-body Hamiltonian (1) with the FGH method, and we compare the dynamics seeded by a low-lying initial state $|{\Psi }_{n=0}^{\mathrm{loc}}(t=0)\rangle$ to that of an initial state $|{\Psi }_{n=3}^{\mathrm{loc}}(t=0)\rangle$ with energy close to the potential’s saddle point, i.e., $E^{2P}\simeq 20$, see Figure 1a.

In the non-interacting case, the wave function always remains separable and, therefore, one can straightforwardly express the probabilities (8) in terms of the single-particle density, which yields

$$\begin{array}{cc}\hfill P_{(LL)}(t)& =P_L^2(t)={\left[{\int }_{x_{\min }}^0\mathrm{d}x{|\psi (x;t)|}^2\right]}^2,\hfill \\ \hfill P_{(LR)}(t)& =2·P_L(t)P_R(t),\hfill \\ \hfill P_{(RR)}(t)& =P_R^2(t)={\left[{\int }_0^{x_{\max }}\mathrm{d}x{|\psi (x;t)|}^2\right]}^2.\hfill \end{array}$$

(15)

Applying a simplified three-level model for $n=0$ [65], Equation (15) can be rewritten as

$$\begin{array}{cc}\hfill P_{(LL)}(t)& ={sin}^4\left(\frac{\Delta }{2}t\right),\hfill \\ \hfill P_{(LR)}(t)& =\frac{1}{2}{sin}^2\left(\Delta t\right),\hfill \\ \hfill P_{(RR)}(t)& ={cos}^4\left(\frac{\Delta }{2}t\right).\hfill \end{array}$$

(16)

Due to the equidistance between the low-lying energies $E_2^{2P},E_1^{2P},$ and $E_0^{2P}$ for $\lambda =0$, the uncorrelated tunneling dynamics is governed by a single Rabi frequency $\Delta =E_2^{2P}-E_1^{2P}=E_1^{2P}-E_0^{2P}$ [43,44]. In particular, for $A_{\max }=10$, $P_{(LL)}(t)$ and $P_{(RR)}(t)$ oscillate with the period $T(\lambda =0)=2\pi /\Delta \simeq 12·{10}^3$. In the non-interacting case, tunneling is thus a first-order single-particle process.

A finite interaction strength perturbs the equidistance between the low-lying energies and, therefore, two distinct periods emerge from the dynamics: $T_{21}=2\pi /(E_2^{2P}-E_1^{2P})$ and $T_{10}=2\pi /(E_1^{2P}-E_0^{2P})$, in qualitative agreement with experimental observations [31].

The evolution of these periods with $\lambda$, plotted in Figure 8, shows a rapid increase (decrease) of $T_{21}$ ($T_{10}$) for weak interactions $\lambda <0.5$, and a monotonous decrease of $T_{21}$ for $\lambda >0.5$, while $T_{10}$ saturates at $T_{10}\simeq 3$ for $\lambda \to \infty$. Please note that for $\lambda =0.5$, the Josephson oscillation period $T_{21}\sim 1750·T(\lambda =0)$ is much larger than the one for non-interacting particles—but finite. This corresponds to the self-trapping regime [4]. Interestingly, the Josephson oscillation period $T_{21}$ converges to the non-interacting period, $T_{21}\sim T(\lambda =0)$, in the Tonks-Girardeau limit $\lambda \to \infty$. This effect agrees with the fermionized pair-state dynamics discussed in Refs. [43,44].

In the two-mode approximation (i.e., the double-well Bose–Hubbard model) for the present scenario, the dynamics is fully described by the amplitudes of the Fock basis states $|n_L,n_R\rangle \in \{|2,0\rangle ,|1,1\rangle ,|0,2\rangle \}$, with degenerate $|2,0\rangle$ and $|0,2\rangle$. Two correlated two-particle tunneling processes are then possible in this simplified picture: a first-order two-particle tunneling process which corresponds to the direct tunneling of both bosons along the diagonal $x_1=x_2$ (i.e., the transition $|2,0\rangle \to |0,2\rangle$), or a second-order two-particle process (i.e., the transition $|2,0\rangle \to |1,1\rangle \to |0,2\rangle$). We now elucidate the actual nature of the tunneling process for weak interactions.

Starting in the initial state $|{\Psi }_{n=0}^{\mathrm{loc}}(t=0)\rangle$ as defined by (14), with $\lambda =0.005$, the dynamics clearly exhibits the Josephson oscillation period $T_{21}\simeq 65·{10}^3$, garnished by a small amplitude beat frequency associated with $T_{10}\simeq 2·{10}^3$. These oscillations are observed in the time evolution of the detection probabilities (8) in Figure 9a, with the Josephson oscillation period $T_{21}\simeq 5.5·T(\lambda =0)$ strongly enhanced with respect to the non-interacting value $T(\lambda =0)$. This is in good qualitative agreement with experimental observation [31]. One also encounters a strongly reduced probability to observe the bosons in opposite wells, signaled by $\max (P_{(LR)})<0.1$ in Figure 9a. The reduction of $\max (P_{(LR)})$, arising from the interaction between the particles, suggests a direct tunneling along the diagonal $x_1=x_2$, i.e., a first-order tunneling process. Such a reduction, which is a corollary of $P_2(t)\equiv {\int }_{x_1·x_2\ge 0}\mathrm{d}{x}_1\mathrm{d}{x}_2{|\psi (x_1,x_2;t)|}^2=P_{(LL)}(t)+P_{(RR)}(t)=1-P_{(LR)}(t)⪅1$, was previously discussed in Refs. [43,44]. However, its interpretation as evidence of first-order tunneling is in contradiction with the time dependence of the integrated probability current $J_{(RR\to LR)}(t)$ also shown in Figure 9a, which clearly indicates a transport across the domain $(LR)$. Indeed, $J_{(RR\to LR)}(t)$ records all probability which passes $(LR)$ and excludes the tunneling along the diagonal $x_1=x_2$. This quantity thus allows us to discriminate sharply the two types of two-particle tunneling. By virtue of Figure 9a, $J_{(RR\to LR)}(t)\sim P_{(LL)}(t)$ implies that almost all probability that oscillates between regions $(LL)$ and $(RR)$ passes region $(LR)$. This confirms a second-order rather than a direct first-order two-particle tunneling process from region $(LL)$ to $(RR)$.

An explanation of the underlying mechanism follows from the expansion coefficients of $|{\Psi }_{n=0}^{\mathrm{loc}}(t=0)\rangle$ in the interacting two-particle basis. The inset in Figure 9b shows that for non-interacting particles, only three coefficients—associated with equidistant energies—are non-zero, giving rise to the single frequency $\Delta =E_2^{2P}-E_1^{2P}=E_1^{2P}-E_0^{2P}$ oscillations described above. Turning on a weak interaction (e.g., $\lambda \le 0.005$, in Figure 9a), the initial state’s overlap with the ground state decreases, while at the same time, the coefficients of the first two excited states pick up comparable weights (squares and diamonds in the inset). The mechanism behind the observed tunneling process is straightforward: in the previous Section, we showed that the first two excited states stick together to form a doublet with an energy which increases with $\lambda$, while the energy of the ground state—one particle localized on each well—does not depend on the interaction, cf. Figure 3. Therefore, the ground state corresponding to a balanced population in region (LR), see Figure 6, becomes off-resonant. Thus, if a boson tunnels from the right- to the left-hand side, it can populate the ground state only for very short times. The associated timescale is determined by the energy gap between the ground state and the degenerate excited states’ energy. Subsequently, the boson tunnels either back to the right well, or the other boson tunnels from the right to the left well, to re-establish energy conservation. It follows from this latter argument that the involved frequencies can be inferred from a three-level model [65]. Increasing further the interaction, the excited states turn resonant with the next higher-lying band (recall Figure 2b and Figure 3), such that additional transitions kick in, and the tunneling dynamics exhibits more frequencies, with no simple representation in the above three-level model. In terms of the expansion coefficients, this boils down to an increasing number of contributing eigenstates as illustrated, for $\lambda =1$, by the triangles in Figure 9b.

Considering now the non-interacting, excited initial state $|{\psi }_{n=3}^{\mathrm{loc}}(t=0)\rangle$ (see Equation (14)) with energy close to the saddle point, i.e., $E^{2P}\simeq 20$, the uncorrelated tunneling dynamics (not shown) is that of a separable wave function with a single Rabi frequency $\Delta =E_{59}^{2P}-E_{52}^{2P}=E_{52}^{2P}-E_{51}^{2P}$, and period $T(\lambda =0)=2\pi /\Delta \simeq 19.5$. This monochromaticity again is a consequence of the equidistant level spacing of the high-lying energies $E_{59}^{2P},E_{52}^{2P}$, and $E_{51}^{2P}$, for $\lambda =0$ (see circles inset Figure 10b). Please note that the Rabi period $T\simeq 19.5$ is much smaller than the one observed for the initial condition $|{\psi }_{n=0}^{\mathrm{loc}}(t=0)\rangle$, for which $T\simeq 12·{10}^3$, since $E_{52}^{2P}-E_{51}^{2P}>E_1^{2P}-E_0^{2P}$, and the detection probabilities, Equation (8), oscillate with reduced amplitude (smaller than 1), due to a less pronounced localization of $|{\psi }_{n=3}^{\mathrm{loc}}(t=0)\rangle$ in either one of the individual wells.

How do interactions affect the evolution of the initial state $|{\Psi }_{n=3}^{\mathrm{loc}}(t=0)\rangle$? As expected from our above spectral analysis, much stronger interactions than $\lambda =0.005$ must be considered to induce visible effects in the dynamics, since the impact of interactions is comparable for all eigenstates (cf. Figure 6, for $A_{\max }=30$ and $n=76$) which exhibit a large overlap with the initial state. Figure 10a shows the time evolution of the detection probabilities (8) for $\lambda =0.1$. The oscillation period seeded by $|{\Psi }_{n=3}^{\mathrm{loc}}(t=0)\rangle$ appears to be much less sensitive to interactions than for $|{\Psi }_{n=0}^{\mathrm{loc}}(t=0)\rangle$ (recall Figure 9): the oscillation periods of $P_{(LL)}(t)$ and $P_{(RR)}(t)$ almost coincide with the non-interacting period $T(\lambda =0)\simeq 19.5$ indicated by vertical dashed lines. Nevertheless, a small shift is visible after seven periods around $t\simeq 136.5$. This small shift can be understood by inspection of the expansion coefficients of the initial state in the interacting two-body eigenbasis, Figure 10b. In contrast to $\lambda =0$, where only three energy levels contribute to the dynamics (circles, inset Figure 10b), an interaction $\lambda =0.1$ redistributes the amplitudes over four dominant states with a weight larger than 5% (squares, inset Figure 10b). The interactions slightly modify the energy gaps, leading to a small modification of the Josephson period, and give finite weight to one additional eigenstate, leading to a modulation of the plotted observables with period $T\simeq 394$. This additional modulation of the signal must not be confused with the damping of density oscillations as observed for large particle numbers in bosonic Josephson junctions [40,68]. As indicated by the time-integrated probability current which roughly follows $P_{(LL)}(t)$ in Figure 10a, we again witness a second-order tunneling across region (LR), instead of direct first-order tunneling along the diagonal $x_1=x_2$. When further increasing the interaction, see, e.g., the diamonds for $\lambda =1$ in Figure 10b, significantly more states contribute to the time evolution (not shown). The inter-particle interaction enforces mixing of the dynamics in the reduced single-particle subspace, and, accordingly, increases the single-particle entropy.

Before the investigation of the time-dependent double-well, we stress here that improved two-mode models for modeling the dynamics of interacting ultracold bosons confined in double-well potentials [69,70] are not sufficient to capture the dynamics as seeded by highly excited initial states. Furthermore, the effective Hamiltonians in [69,70] contain strongly initial-state-dependent parameters and thus a comprehensive comparison of different initial states is considerably complicated.

4.2. Time-Dependent Double-Well Potential: From Few- to Many-Body Dynamics

We have seen in the previous sections how the barrier height affects the impact of interactions on the many-particle dynamics. We now generalize this analysis by considering a time-dependent switching of the barrier according to Equations (2) and (3), with $A_{\max }=30$. Before this quench, the bosons are prepared in the interacting many-particle ground state of a harmonic trap. Our purpose is here to examine how the reduced one-body density matrix evolves for (quasi)-adiabatic vs. diabatic switching. Extrapolation to larger particle numbers using the MCTDH-X method relates our observations to previous studies of the splitting of a BEC by a laser sheet [34,39,71]. Please note that while quenches can be efficiently simulated with the help of the FGH method, we employ the MCTDH-X method (see Appendix C) for finite switching times, to deal with the time-dependent Hamiltonian (1).

We start with the time evolution of the many-body wave function when the tunneling barrier is suddenly quenched from $A_{\max }=0$ to 30 (i.e., $T_{\mathrm{ramp}}\to 0$) [35,72,73]. Figure 11a–d shows the behavior of the two-particle density for $\lambda =1$, during the initial stage of the quench-induced dynamics. The initial wave packet is split along the diagonal $x_1=x_2$, and spreads towards the outer edges of the double-well, until its reflection after half a period $t\simeq 1.9$. Since all the injected energy, i.e., $A_{\max }=30$, is suddenly transferred to the two bosons, the turning point $x_{\mathrm{turn}}\sim \pm 7.75$ in Figure 11c, where the reflection takes place, corresponds to $V(x_{\mathrm{turn}})\simeq A_{\max }=30$ (see Figure 1b). We observe (not shown) that the higher the tunneling barrier $A_{\max }$, the longer the oscillation period. On its way back, the wave packet broadens more and more due to reflections at the central barrier. Finally, after one period $t\simeq 3.6$, Figure 11d, a large fraction is again located in the vicinity of the saddle point, which, subsequently, splits once more.

In contrast, for a long ramping time $T_{\mathrm{ramp}}=30$, see Figure 11e–h, the wave function has enough time to adapt to the new boundary conditions, such that it rather smoothly follows the minima of the dynamically created double-well potential. The dynamics are still garnished, for this long but finite ramping time, by excitations of the first band, as identifiable by additional nodal structures in Figure 11g,h.

Comparison of the nodal structures of the two-particle densities observed in Figure 11 for $T_{\mathrm{ramp}}=0$ and for $T_{\mathrm{ramp}}=30$, respectively, suggests that less energy is absorbed by the center-of-mass degree of freedom in the latter case (as expressed by considerably fewer nodal lines, indicative of smaller momenta). To corroborate this conjecture (which is based on evidence exclusively gathered from short time dynamics), we plot the two-particle energy expectation value

$$\begin{array}{c}\hfill E^{2P}(T_{\mathrm{ramp}})=\langle \Psi (t_0)|{\mathcal{H}}^{2\mathrm{P}}(T_{\mathrm{ramp}})|\Psi (t_0)\rangle ,\end{array}$$

(17)

at $t_0=200\gg T_{\mathrm{ramp}}$, for variable $T_{\mathrm{ramp}}\in [0,30]$, in Figure 12. We observe a quick initial drop of the energy followed by a long tail approaching smoothly the energy of the ground state, for $A_{\max }=30$ and $\lambda =1$, i.e., $E_0^{2P}\simeq 11.34$. The inset zooms into the range $T_{\mathrm{ramp}}\in [7.5,30.5]$, where the horizontal dashed lines indicate the eigenenergies of the two-particle system, with $A_{\max }=30$ and $\lambda =1$. The evolution of $E^{2P}(t_0)$ implies that for $T_{\mathrm{ramp}}\ge 19$, only transitions between the ground state and the first degenerate (recall Figure 6) excited states occur. Thus, indeed, (quasi-)adiabatic switching does perform essentially no work on the many-particle system.

The static double-well’s entropy of the reduced single-particle density matrix increases from zero at $\lambda =0$ and saturates at $ln2$ [37,38,56] with our definition (6) for $\lambda \to \infty$, $\forall A_{\max }$. In our present, dynamical scenario—where the harmonic trap is split into a double-well during a time $T_{\mathrm{ramp}}$—we also expect the entropy to increase with the interaction. Figure 13 shows the time evolution of the entropy for two ramping times (red/blue) and for two values of the interaction strength, (a) $\lambda =1$ and (b) $\lambda =0.1$.

For short ramping time, $T_{\mathrm{ramp}}=0.001$ (red lines), the entropy increases and saturates at $\approx 2.51$ which is well below the maximal value $S^{\max }=log(M)\approx 2.77$ and which we verified with respect to the time evolution for $T_{\mathrm{ramp}}=0$ using the spectral decomposition based on our FGH computations from Section 3. In agreement with the asymptotic behavior of the energy expectation value observed in Figure 12, the entropy oscillates with a single frequency for large ramping time, e.g., $T_{\mathrm{ramp}}=30$ (blue lines in Figure 13). The stronger the interaction, the larger the frequency as well as the offset of the minima of the entropy oscillations.

Monitoring the time evolution of the entropy over a broad interval of $T_{\mathrm{ramp}}$ allows us to map out the different dynamical regimes for two bosons with $\lambda =1$ and $\lambda =0.1$, respectively, see Figure 14a,d.

In full agreement with what we observed above for the dependence of the energy expectation value on $T_{\mathrm{ramp}}$, the transition from diabatic to (quasi-)adiabatic dynamics is also here the primary feature: For short ramping times, the entropy rapidly saturates at its equilibrium value, whereas, for a sufficiently slow ramp, an oscillation emerges, with a single, well-defined frequency (and decreasing amplitude, for increasing $T_{\mathrm{ramp}}$). A discrete Fourier transform of the signal for $T_{\mathrm{ramp}}\approx 19-30$ shows that this frequency is determined by the energy gap (see Section 3.1) between the ground and first excited state,

$$\begin{array}{c}\hfill \nu (\lambda )=\frac{E_1^{2\mathrm{P}}(\lambda )-E_0^{2\mathrm{P}}}{\pi }.\end{array}$$

(18)

Indeed, only two eigenstates of the reduced single-particle density matrix—with opposite parity and densities which closely resemble the typical structure of the double-well ground state doublet [56]—contribute to the dynamics in this oscillating region (not shown). For intermediate ramping times ($T_{\mathrm{ramp}}\approx 10-19$), the structures observed in Figure 14a,d still express the switching-induced, coherent coupling of more than just two interacting eigenstates, because in this regime the dynamics are not yet (quasi-) adiabatic (in agreement with our discussion of Figure 12).

Remarkably, although the detailed spectral structures are rather different for two and three particles (see Figure 3 and Figure 4), the ramping-induced time dependence of the von Neumann entropy is qualitatively similar for $N=2,3$, and even $N=10$ (where we cannot access the spectral structure, with our presently available numerical resources) see Figure 14a–c, for $\lambda =1$, and Figure 14d–f, for $\lambda =0.1$. We attribute this feature to the coarse-graining effect of a diabatic switch, where only the effective density of states must be gauged against the spectral width of the time-dependent perturbation. Closer inspection suggests that efficient excitation is achieved for slightly longer switching times with increasing particle number, which is consistent with the increase of the density of states with N. The frequency $\nu \sim E_1^{\mathrm{NP}}-E_0^{\mathrm{NP}}$ of the entropy oscillations slowly decreases with the number of particles, since the energy gap $\Delta E=E_1^{\mathrm{NP}}-E_0^{\mathrm{NP}}$ between the first excited state and the ground state decreases with N, i.e., $\Delta E(N=2)>\Delta E(N=10)$. Also note that the oscillating regime of Figure 14c,f corresponds to the two-fold fragmented BEC discussed in Refs. [34,39,40]. Similar results are observed for different barrier heights (not shown) [56].

Let us conclude this section with a remark on the convergence of the MCTDH-X results reported in Figure 14c,f: For moderate and large $T_{\mathrm{ramp}}\gtrsim 7$, only two orbitals of the employed $M=8$ orbitals have a significant population and the entropy S remains significantly smaller than the maximal value $S^{\max }=log(M)$. From this fact it can be inferred that the wave function is accurately described in these MCTDH-X computations at $T_{\mathrm{ramp}}\gtrsim 7$. However, for small ramping times ($T_{\mathrm{ramp}}\lesssim 7$) all $M=8$ employed orbitals in the computation were populated. Consequently, the entropy reaches its maximum $S^{\max }$. This maximal entropy for small $T_{\mathrm{ramp}}$ implies that the Hilbert space provided by MCTDH-X is not large enough to host the complete dynamics of the many-body wave functions and more orbitals ($M>8$) would therefore be necessary to achieve convergence. Based on the FGH results for sudden switches of the potential barrier, to cover the subspace of the Hilbert space more than $M=16$ (corresponding to a maximal entropy of $S^{\max }\approx 2.77$) orbitals are necessary, which exceeds the typically employed number of orbitals ($M\in \{12\dots 16\}$) for $N\le 10$ bosons reported in the literature [74,75]. While the quantitative behavior of the entropy $S(T_{\mathrm{ramp}},t)$ at small $T_{\mathrm{ramp}}\lesssim 7$ in Figure 14c,f therefore cannot be considered fully converged, the observed behavior is qualitatively equivalent to that resulting for smaller particle numbers, where convergence of MCTDH-X could be achieved with a smaller number $M=2,4,6$ of orbitals, and is also consistent with our FGH-based analysis for $N=2$ particles (see Figure 12). This suggests that the results reported in Figure 14c,f correctly indicate the qualitative trend of the evolution also for short ramping times.

5. Conclusions

We analyzed the spectral structures and the dynamics of a few interacting bosons in a one-dimensional double-well potential, for both a static and a time-dependent potential barrier, beyond the two-mode approximation. To this end, we used three complementary numerical methods. The Fourier Grid Hamiltonian method was employed to extract the full spectral information for two interacting bosons, whereas a Bose–Hubbard representation of the continuous double-well potential was found to be more efficient to describe the spectral structure of the three-particle case. Furthermore, we used the MCTDH-X method to simulate the dynamical evolution of $N=2,3$, and 10 interacting bosons in a potential with time-dependent barrier strength.

Our spectral analysis highlights the dependence of the energy spectrum on the interaction strength, on the one hand, and on the potential barrier height, on the other. Ramping up a barrier in the center of an initially harmonic potential introduces a metamorphosis of state space from a simple, highly degenerate harmonic oscillator progression, into a sequence of states which exhibit the characteristic degeneracies associated with tunneling between symmetric wells, below the barrier energy, and a harmonic-like spectrum sufficiently high above the barrier, separated by a range around the barrier energy which mediates between both classes. Interactions lift many of the energetic degeneracies and eventually induce mixing of energetic manifolds which otherwise remain well-separated.

While for two (on-site interacting) particles distributed over two (deep) wells eigenstates exist which remain unaltered by finite interactions, this is no longer true for three particles in the same potential, since at least two particles then must interact: two manifolds of states emerge corresponding to states where two or three particles are interacting. We supplemented our spectral analysis by inspecting many-particle probability densities in configuration space, which directly exhibit the spatial correlations inscribed into the many-body tunneling dynamics, for both energetically low- and high-lying states. For three particles, we visualized the loss of information about correlations when tracing from the three-body density to the two-body, and, eventually, to the one-body density.

We used that spectral information to decipher the tunneling dynamics of two interacting particles in a static double-well. In particular, we compared and characterized Josephson oscillations of two interacting bosons prepared in a superposition of excited states with energies either well below or close to the potential’s saddle point. Inspection of the expansion coefficients of the evolved two-particle state in the interacting two-particle basis provided evidence that a simple three-level description of the dynamics fails at sufficiently strong interactions. The Josephson period at energies close to the saddle point is much smaller and robust with respect to the interaction. In agreement with observations in Ref. [31], we confirm a second-order pairwise tunneling process.

Finally, we investigated the spreading behavior of the many-particle state, when initially prepared in the many-particle harmonic oscillator ground state, under diabatic vs. (quasi-) adiabatic switching of a central barrier—transforming the potential into a double-well. Diabatic switching leads to efficient energy transfer through the population of many many-particle excited states, as quantified by the time evolution of the von Neumann entropy of the reduced single-particle density matrix, while a (quasi-) adiabatic ramp only populates weakly excited states. This phenomenology emerges already for two interacting particles and—due to the increasing spectral density—gets more pronounced for ten particles, the largest particle number here considered.