Many-Body Quantum Dynamics of a Bosonic Josephson Junction with a Finite-Range Interaction

TL;DR

This paper investigates how finite-range interactions influence the out-of-equilibrium quantum dynamics of a Bose gas in a double well, revealing that interaction range affects system behavior differently in each well and requires many-body theory for accurate description.

Contribution

It provides a detailed numerical analysis of finite-range interaction effects on Bose gas dynamics in a double well, highlighting the necessity of many-body approaches like MCTDHB.

Findings

Interaction range impacts dynamics differently in each well.

Finite-range effects are significant even in the infinite-particle limit.

Many-body theory is essential to accurately describe the system.

Abstract

The out-of-equilibrium quantum dynamics of a Bose gas trapped in an asymmetric double well and interacting with a finite-range interaction has been studied in real space by solving the time-dependent many-body Schr\"odinger equation numerically accurately using the multiconfigurational time-dependent Hartree method for bosons (MCTDHB). We have focused on the weakly interacting limit where the system is essentially condensed. We have examined the impact of the range of the interaction on the dynamics of the system, both at the mean-field and many-body levels. Explicitly, we have studied the maximal and the minimal values of the many-body position variance in each cycle of oscillation, and the overall pace of its growth. We find that the range of the interaction affects the dynamics of the system differently for the right well and the left well. We have also examined the infinite-particle limit and find that even there, the impact of the range of the interaction can only be described by a many-body theory such as MCTDHB.

Figures

Equations

$$V_T(x)=\{\begin{array}{c}\frac{1}{2}(x+2{)}^2+Cx,x<-\frac{1}{2}\\ \\ \frac{3}{2}(1-x^2)+Cx,|x|\le \frac{1}{2}\\ \\ \frac{1}{2}(x-2{)}^2+Cx,x>\frac{1}{2}\end{array}.$$

$$\frac{1}{N}{\Delta }_{\hat{X}}^2$$

$${\rho }^{(1)}(x_1|x_1^{\prime };t)$$

$${\rho }^{(2)}(x_1,x_2|x_1^{\prime },x_2^{\prime };t)=N(N-1)\int d{x}_3\dots d{x}_N{\Psi }^{\ast }(x_1^{\prime },x_2^{\prime },x_3,\dots ,x_N;t)\Psi (x_1,x_2,x_3,\dots ,x_N;t).$$

$$|\Psi (t)⟩=\sum _{\vec{n}}{C}_{\vec{n}}(t)|\vec{n};t⟩.$$

Full text

1 Introduction

Following the experimental observation by the MIT group [1], the Bose-Einstein condensate (BEC) in an ultra-cold atomic gas has emerged as a popular test-bed for the interacting many-particle system [2]. One of the main reasons for such huge interest in this system is its high degree of tunability and thereby its ability to simulate various condensed matter systems. One such well-studied system in condensed matter physics is the Josephson junction [3]. An interacting Bose gas in a symmetric double well provides a paradigmatic model for the Josephson junction and is popularly known as the bosonic Josephson junction (BJJ) [4]. The BJJ has been extensively studied both theoretically and experimentally [4-9]. Theoretically, two commonly used methods for the study of the BJJ dynamics are the mean-field method [2] and the Bose-Hubbard model [10]. While there are several features of the BJJ dynamics, such as the collapse and revival of the density oscillations [11] and the fragmentation dynamics [12], which can not be described by the mean-field method but can be described by the Bose-Hubbard model, the latter also fail to grasp the full many-body features [13].

Recently, a numerically-exact many-body method has been used to study the BJJ dynamics for both the repulsive as well as attractive interactions [14]. Specifically, fragmentation and the uncertainty product of the many-body position and momentum operators have been studied [15, 16]. It has also been shown that though there is a symmetry between the time evolution of the Bose-Hubbard dimer for the repulsive and attractive interaction of the same interaction strength, no such symmetry exists for the many-body Hamiltonian [17]. Moreover, a universality of the fragmentation dynamics in the sense that the condensates with different number of bosons but same interaction parameter fragment to the same degree has been found [12].

Though all these works considered a contact $\delta-$interaction, a typical finite-range interaction has also been considered in a recent study and the impact of the range of the interaction on the BJJ dynamics has been examined [18]. Such studies are important in view of the recent experiments with the ultra-cold dipolar atoms where it has been clearly shown that both the short-range as well as the long-range terms are necessary to accurately describe the observed physics. BJJ dynamics was observed to intricately depend on the competitive effects of the range of the interaction and the trapping potential. The tail of the interaction potential was found to effectively enhance the effect of the repulsive interaction until the range became comparable with the length scale set by the trap.

Further, in a recent work, an asymmetric double well has also been considered and the effect of the loss of symmetry on the BJJ dynamics has been thoroughly studied [19]. It has been found that there is an overall suppression of the density oscillations in the asymmetric double well. Also, for a sufficiently strong interaction, the condensate becomes fragmented with time and the degree of fragmentation depends on the asymmetry as well as the initial well in a non-trivial manner. The universality of the fragmentation is also observed in the asymmetric double well though the degree of fragmentation depends on the initial well.

In this work, we look back at the problem of the BJJ dynamics and combine the issues of the finite-range interaction and the asymmetric trap in order to understand their combined impact on the BJJ dynamics. Hence, in this work, we consider the BJJ dynamics in an asymmetric double well for a typical finite-range interaction. We examine how the BJJ dynamics is affected by the tail of the interaction in presence of an asymmetry. We focus on the weakly interacting regime where the system is essentially fully condensed. Further, we examine the large particle limit and look for the combined effect of the range of interaction and the asymmetry of the trap at the many-body level. We observe that even at the effectively infinite-particle limit, where the system is $100\%$ condensed, there are interesting effects at the many-body level. The structure of the paper is as follows. In section 2, we describe the system considered here. We present our results in Section 3 and summarize in Section 4. The methodology and its accuracy are demonstrated in the appendices.

2 System

Here we are interested in a gas of $N$ interacting structureless bosons confined in a one-dimensional (1D) asymmetric double well $V_{T}(x)$. $V_{T}(x)$ is constructed by fusing two shifted harmonic wells $V_{L,R}(x)=\frac{1}{2}(x\pm 2)^{2}+Cx$, i.e.,

[TABLE]

The time period of Rabi oscillation $t_{Rabi}=\frac{2\pi}{\Delta E}$ [$\Delta E$ being the energy gap between two lowest levels of $V_{T}(x)$] in $V_{T}(x)$ provides the natural time scale of the dynamics. In this work we consider only a small asymmetry $C=0.01$. For such a trap $t_{Rabi}=102.0618$. The spatial separation $l=4$ between the two local minima of the asymmetric double well trap sets the length scale for the system. The two-body inter-atomic interaction of strength $\lambda_{0}$ is modeled as $W(x_{j}-x_{k})=\frac{\lambda_{0}}{\sqrt{(|x_{j}-x_{k}|/D)^{6}+1}}$ of half-width $D$. $\lambda_{0}$ corresponds to the interaction parameter $\Lambda=\lambda_{0}(N-1)$. For a fixed $\Lambda$, the range of the interaction can be tuned by varying $D$ as a factor of the length scale $l$. Note that for $\frac{|x_{j}-x_{k}|}{D}\gg 1$, $W(x_{j}-x_{k})$ behaves as the long-range dipolar interaction while in the limit $\frac{|x_{j}-x_{k}|}{D}\rightarrow 0$, $W(x_{j}-x_{k})$ reduces to a soft-core interaction having many similar effects of a short-range $\delta$-interaction.

We characterize the dynamics of the system mainly by the time evolution of the many-body position variance of the system. The main reason for focusing on the time evolution of the many-body position variance is the recent observation that it can differ from the predictions of the mean-field Gross-Pitaevskii (GP) theory even when the system is fully condensed [20]. Therefore, the many-body position variance can grasp those many-body effects which are otherwise washed out in the infinite-particle limit when the system is $100\%$ condensed and its density per particle and the energy per particle can be accurately described by the mean-field GP theory. Further, we also study the growth in the number of depleted particles with time to check the degree of condensation of the system.

Given the many-body wavefunction $\Psi(t)$, one can define the variance of the many-body position operator $\hat{X}=\sum_{j=1}^{N}\hat{x}_{j}$ as

[TABLE]

Here, $\rho(x)$ is the density of the system and is given by the diagonal element of the reduced one-body density matrix, i.e, $\rho(x;t)\equiv\rho^{(1)}(x|x^{\prime}=x;t)$. Similarly, $\rho_{ksql}$ is the matrix elements of the two-body reduced density matrix, viz., $\rho_{ksql}=\left<\Psi\left|b_{k}^{\dagger}b_{s}^{\dagger}b_{q}b_{l}\right|\Psi\right>$ where $b_{k}$ and $b_{k}^{\dagger}$ are the bosonic annihilation and creation operators, respectively. The reduced one-body and two-body density matrices $\rho^{(1)}(x_{1}|x_{1}^{\prime};t)$ and $\rho^{(2)}(x_{1},x_{2}|x_{1}^{\prime},x_{2}^{\prime};t)$, respectively, are defined as

[TABLE]

and

[TABLE]

We point out that the eigenvalues and eigenvectors of $\rho^{(1)}(x_{1}|x_{1}^{\prime};t)$ are the time-dependent natural occupation numbers $n_{j}(t)$ and the time-dependent natural orbitals $\phi^{NO}_{j}(x_{1},t)$, respectively. Further, $n_{j}(t)$ characterize the (time varying) degree of condensation in a system of interacting bosons [21] and satisfy $\sum_{j=1}^{M}n_{j}=N$ ($M$ is the number of single particle orbitals used to construct the many-boson wavefunction, see A). If only one macroscopic eigenvalue $n_{1}(t)\approx{\mathcal{O}}(N)$ exists, the system is condensed [21] whereas if there are more than one macroscopic eigenvalues, the BEC is said to be fragmented [15,22-26]. For a condensed system, the microscopic occupations in all the higher orbitals $f=\sum_{j=2}^{M}\frac{n_{j}}{N}$ gives the depletion. On the other hand, the macroscopic occupation of a higher natural orbital $f=\frac{n_{j>1}}{N}$ gives the fragmentation. Finally, we remark that substituting the many-body wavefunction $\Psi(t)$ in Eq. (2) by the corresponding mean-field wavefunction, one can in study the variance of the many-body position operator at the mean-field level.

3 Results

In this section we describe the findings of our study of the 1D BJJ dynamics in an asymmetric trap. In an asymmetric double well trap, two wells are no longer equivalent - for the particular double well trap considered here, the left well $V_{L}(x)$ is lower than the right well $V_{R}(x)$. Accordingly, one can prepare the initial condensate state either in $V_{L}(x)$ or in $V_{R}(x)$. In this work, we start the out-of-equilibrium dynamics of the system once from $V_{L}(x)$ and then from $V_{R}(x)$. As mentioned above, we consider the time evolution of the many-body position variance up to about $30$ Rabi periods, i.e., more than 3000 in absolute time units. We consider systems of $N=100$ and $10000$ bosons interacting via a finite-range interaction of strength $\lambda_{0}$.

As already mentioned, we are interested only in the very weakly interacting limit. Accordingly, we consider $\Lambda=0.01$. At such weakly interacting limit, the system is expected to be essentially condensed and, therefore, its out-of-equilibrium dynamics should be adequately described by the mean-field GP theory. So, first we consider the time evolution of the position variance at the mean-field level. The corresponding symmetric double well case, as discussed in Ref. [18], will serve as a reference for our study.

We plot the time evolution of $\frac{1}{N}\Delta_{X}^{2}(t)$ for different ranges of interaction for the left and the right well in Fig 1(a) and (b), respectively. We observe that for both wells, $\frac{1}{N}\Delta_{X}^{2}(t)$ oscillates smoothly with a frequency equal to the Rabi frequency. This is in contrast to the symmetric double well where the frequency of the oscillations in $\frac{1}{N}\Delta_{X}^{2}(t)$ is twice the Rabi cycle. Actually, unlike in the symmetric double well, the density never completely tunnels out of the initial well of the asymmetric double well trap. This is manifested in the broadening of the peaks of oscillations in $\frac{1}{N}\Delta_{X}^{2}(t)$. Also, the minima of the oscillations, after starting from $0.5$ is always slightly higher than $0.5$. Further, the curves for different ranges $D$ of the interaction overlap with each other for both the wells, though a very small deviation develops for the right well at a longer time. This shows that there is no significant effect of the ranges of interaction at the mean-field level.

Next, we study the time evolution of the $\frac{1}{N}\Delta_{X}^{2}(t)$ at the many-body level. We computed the many-body results by the MCTDHB method with $M=2$ orbitals, see A for details. In Fig. 2(a) and (b) we have plotted the time evolution of the many-body position variance of a system of $N=100$ interacting bosons for various ranges of the interaction, for the left well and the right well, respectively. We observe that $\frac{1}{N}\Delta_{X}^{2}$ oscillates with the Rabi frequency between a growing maximum and a minimum which is approximately equal to $0.5$. The deviation from $0.5$ is because of the remnants in each well as discussed above. This is further manifested in the splitting of the peaks into two sub peaks. Moreover, the difference between the two sub peaks grows with time for both wells. We also see a clear impact of the range $D$ of the interaction on the time evolution of $\frac{1}{N}\Delta_{X}^{2}(t)$ though $D$ affects the the time evolution of the many-body position variance $\frac{1}{N}\Delta_{X}^{2}(t)$ differently for the left $V_{L}(x)$ and the right $V_{R}(x)$ wells. While for $V_{L}(x)$, the peak values of the oscillations in $\frac{1}{N}\Delta_{X}^{2}(t)$ decrease with increasing $D$ until $D=l/2$ for which the peak values are minimum. On the other hand, similarly to the symmetric double well, the peak values of $\frac{1}{N}\Delta_{X}^{2}(t)$ increase with increasing $D$ until $D=l/2$ for which it reaches the maximum. However, here, the minima also increases with time for all $D$. Further, while the higher sub peak appears on the left side for $V_{L}(x)$, it is on the right side for $V_{R}(x)$.

Now, to confirm that at this interaction strength for $N=100$, the system is essentially condensed, we have also plotted the number of depleted atoms outside the condensed mode as a function of time for the right as well as the left wells in Fig. 3(a) and (b), respectively. We observe that the number of depleted atoms grows with time in an oscillatory manner. We also note that the growth rates of the number of depleted atoms are different for $V_{L}(x)$ and $V_{R}(x)$. Moreover, for $V_{L}(x)$ itself, these growth rates are different depending on the range $D$ of the interaction. Accordingly, at short time, the system is the most depleted for $D=l/2$ whereas at longer time, the depletion is the least for $D=l/2$. Also for the left well, the depletions of the system for all $D$ are approximately of the same order, though they vary a bit depending on $D$. On the other hand, for $V_{R}(x)$ the growth rate of the depletion is independent of $D$. However, the number of depleted particles varies widely for different $D$. Also, at any point of time, the number of depleted particles increases with $D$ and is maximum for $D=l/2$. With further increase in $D$, the number of depleted particles decreases. Further, for all $D$, the system is more depleted for $V_{R}(x)$ than $V_{L}(x)$. Finally, for all cases, the number of depleted particles are negligibly small implying that the system is essentially condensed and the mean-field GP theory should be adequate to describe the dynamics of the system. Thus, the effects of the range $D$ of the interaction observed above in the time evolution of $\frac{1}{N}\Delta_{X}^{2}$ appear only at the many-body level, even when the mean-field theory should be applicable for the description of the out-of-equilibrium dynamics of the system’s density per particle.

Rigorous results exist in the literature proving that, in the infinite-particle limit keeping $\Lambda$ fixed, the density per particle and the energy per particle of the ground state of the system coincide with the corresponding mean-field GP results [27-30]. Therefore, next we examine if the effects of the ranges of interaction and the asymmetry of the trap, observed at the many-body level for a finite size system, still persist for a much larger system size en route to the infinite-particle limit. Hence, next we consider a system of $N=10000$ bosons keeping $\Lambda=0.01$ fixed. We again study the time evolution of the many-body position variance of the system by MCTDHB method with $M=2$ orbitals. The $\frac{1}{N}\Delta_{X}^{2}(t)$ for various ranges of interaction $D$ for both $V_{L}(x)$ and $V_{R}(x)$ are shown in Fig. 2(c) and (d), respectively. We observe essentially identical behavior in the oscillatory growth of $\frac{1}{N}\Delta_{X}^{2}(t)$ with time for $V_{L}(x)$ as for the system with $N=100$ particles. However for $V_{R}(x)$, though $\frac{1}{N}\Delta_{X}^{2}(t)$ exhibit an oscillatory growth similarly to that of a system with $N=100$ bosons, here the minima of $\frac{1}{N}\Delta_{X}^{2}(t)$ do not grow with time and is only slightly higher than $0.5$. We also note that for the right well, the time evolution of $\frac{1}{N}\Delta_{X}^{2}(t)$ for all $D$ are qualitatively similar to the corresponding symmetric double well case except that the growth rate, frequencies, and amplitudes are different.

Since, as already mentioned, the many-body position variance depends on the actual number of depleted atoms, next we also plot the corresponding number of depleted particles for $N=10000$ as a function of time for both the left and the right wells in Fig. 3(c) and (d), respectively. We find that the oscillatory growth of the number of depleted particles are practically identical for both $N=100$ and $N=10000$ bosons. Further, at any point of time, the corresponding numbers of depleted atoms for all $D$ and wells are also essentially identical for both $N=100$ and $N=10000$. This implies that at any point of time the system with $N=10000$ is actually two orders of magnitudes less depleted and therefore closer to the mean-field approximation. Thus the practically identical behaviors of $\frac{1}{N}\Delta_{X}^{2}(t)$ for $N=100$ and $N=10000$ for all cases indicate that the effects of the range $D$ of interaction have converged with respect to $N$ and will, therefore, persist in the infinite-particle limit also. This implies that the effects of the range $D$ of the interaction on the asymmetric BJJ observed here are many-body in nature and can not be described by the mean-field theory even when the system is $100\%$ condensed.

4 Summary

In summary, we have studied the impact of a finite range interaction on the dynamics of an interacting Bose gas in a one dimensional asymmetric double well trap. We have characterized the dynamics by the time evolution of the many-body position variance $\frac{1}{N}\Delta_{X}^{2}(t)$ which can deviate from the corresponding mean-field results even when one out of a million atoms is outside the condensed mode. We focused only on the weakly interacting regime where the system is essentially fully condensed and the mean-field theory should be adequate for the description of its the out-of-equilibrium dynamics. Also, we examined the infinite-particle limit where the system is $100\%$ condensed and its ground state energy per particle and the ground state density per particle coincide with the corresponding mean-field GP results. We find that, even for a weakly interacting system, the effect of the range of interaction on the dynamics of a one dimensional asymmetric BJJ can be described only at the many-body level. Further, such many-body effects are also found to persist in the large particle systems when the infinite-particle limit is essentially achieved. Also, the range of the interaction affects the dynamics differently depending on whether the initial condensate state is prepared in the left well or in the right well. While the dynamics for the right well is qualitatively similar to the symmetric BJJ dynamics, differences appear in the growth rate, frequencies, and peak values of $\frac{1}{N}\Delta_{X}^{2}(t)$. On the other hand, completely new features are observed for the left well.

\ack

This research was supported by the Israel Science Foundation (Grant No. 600/15). Computation time on the High Performance Computing system Hive of the Faculty of Natural Sciences at University of Haifa is gratefully acknowledged. We also acknowledge computation time on the BwForCluster and the Cray XC40 system Hazelhen at the High Performance Computing Center Stuttgart (HLRS). SKH gratefully acknowledges the continuous hospitality at the Lewiner Institute for Theoretical Physics (LITP), Department of Physics, Technion - Israel Institute of Technology. SKH also acknowledges fruitful discussions with Anal Bhowmik.

Appendix A Multiconfigurational time-dependent Hartree method for bosons (MCTDHB)

In this work, we have used the multiconfigurational time-dependent Hartree method for bosons (MCTDHB) for studying the dynamics of an interacting Bose gas in a one dimensional asymmetric double well. MCTDHB is an in-principle numerically exact method which has been developed [14, 31] for solving the time-dependent many-body Schrödinger equation and benchmarked with an exactly-solvable model [32, 33]. This method has already been extensively used in the literature [13,34-42].

MCTDHB is a time-dependent variational method in which the ansatz is chosen as the superposition of all possible configurations obtained by distributing $N$ bosons in $M$ time-dependent single-particle orbitals, viz.,

[TABLE]

Here the occupations $\vec{n}=(n_{1},n_{2},\ldots,n_{M})$ preserve the total number of bosons $N$. With this ansatz, the time-dependent action is minimized to obtain the working equations for the orbitals and the coefficients. For details, we refer to Ref. [14]. Here we point out that for $M=1$ orbital, we get the ansatz for the GP equation while in the limit $M\rightarrow\infty$ the ansatz expands over the full Hilbert space and hence gives an exact theory. However, in practice, one has to use a finite number $M$ of orbitals. In our numerical calculations, we keep on repeating the computations with increasing $M$ until the convergence is reached and thereby we obtain the numerically highly accurate results, see B. We point out that because of the time-dependent permanents, one can use a much shorter expansion in Eq. (5) than if only the coefficients are allowed to be time-dependent, and this leads to a significant computational advantage. The ground state properties of the system can be obtained by propagating the MCTDHB equations in imaginary time [15].

Appendix B Details of the numerical computations and their convergence

For our numerical computations, we have used the numerical implementation of the parallel version of MCTDHB [43, 44]. Here we have represented the many-body Hamiltonian by 128 exponential discrete-variable-representation (DVR) grid points using a fast Fourier transformation (FFT) routine within a box of size [-10:10). We obtain the initial many-body ground state of the BEC in both the left and right wells by propagating the MCTDHB equations in imaginary time. As already mentioned above, in our numerical calculations we keep on increasing the number of time-adaptive orbitals $M$ to obtain numerically exact results. Below, we explicitly show the numerical convergence of our results with respect to the number of orbitals $M$ for a system of $N=10$ bosons and $\Lambda=0.01$. Since increasing $N$, keeping $\Lambda$ fixed, amounts to an effectively weaker interaction, the convergences of our results for the systems of $N=100$ and $N=10000$ bosons are ensured by this procedure.

We demonstrate the numerical convergence of $\frac{1}{N}\Delta_{X}^{2}(t)$ for each $D$ separately for both the wells. This also ensures that the small differences observed among $\frac{1}{N}\Delta_{X}^{2}(t)$ and the number of depleted atoms for different $D$ are well converged. In Fig. 4, each panel shows the results for $\frac{1}{N}\Delta_{X}^{2}(t)$ for each $D$ obtained with $M=2,4,$ and $6$ orbitals for both wells. We find that for all cases, $\frac{1}{N}\Delta_{X}^{2}(t)$ exhibit an oscillatory growth with the frequency of oscillation being equal to the Rabi frequency. However, now the minima of $\frac{1}{N}\Delta_{X}^{2}(t)$ grow with time for the left well also in addition to that for the right well. Also, $\frac{1}{N}\Delta_{X}^{2}(t)$ computed with $M=2,4,$ and $6$ practically overlap with each other for all cases with the results for $M=4$ and $6$ lying atop each other. This shows that the use of $M=2$ orbitals already gives the accurate description of the time evolution of $\frac{1}{N}\Delta_{X}^{2}(t)$ for all cases.

Next, we plot the number of depleted atoms for each $D$ and both the left and right wells in each panel of Fig. 5. Here also, we see that for all cases, the number of depleted atoms grow in an oscillatory manner. As discussed in the text, the system is more depleted for the right well than the left well for all $D$. Also, the depletion of the system is found to depend on the range $D$ of the interaction more prominently for the right well. Here also, we find that the number of depleted particles computed with $M=2,4,$ and $6$ practically overlap with each other while the results for $M=4$ and $6$ are again lying on top of each other. This shows that the results obtained with $M=2$ orbitals are already accurate and aptly describe the impact of the asymmetry of the trap and the range of the interaction at the many-body level.