Many-body quantum dynamics of an asymmetric bosonic Josephson junction

The out-of-equilibrium quantum dynamics of an interacting Bose gas trapped in a one-dimensional asymmetric double-well potential is studied by solving the many-body Schrödinger equation numerically accurately. We examine how the gradual loss of symmetry of the confining trap affects the macroscopic quantum tunneling dynamics of the system between the two wells. In an asymmetric double well, the two wells are not equivalent anymore, say, the left well is deeper than the right one. Accordingly, we analyze the dynamics by initially preparing the condensate in both the left and the right wells. The dynamics of the system is characterized by the time evolution of a few physical quantities of increasing many-body complexity, namely, the survival probability, depletion and fragmentation, and the many-particle position and momentum variances. In particular, we have examined the frequencies and amplitudes of the oscillations of the survival probabilities, the time scale for the development of fragmentation and its degree, and the growth and oscillatory behavior of the many-particle position and momentum variances. There is an overall suppression of the oscillations of the survival probabilities in an asymmetric double well, except for resonant values of asymmetry for which the one-body ground state energy in the right well matches with one of the one-body excited states in the left well, thereby resulting in resonantly enhanced tunneling from the right well ground state. Overall, depending on whether the condensate is initially prepared in the left or right well, the repulsive inter-atomic interactions affect the survival probabilities differently. For a sufficiently strong repulsive interaction, the system is found to become fragmented. The degree of fragmentation depends both on the asymmetry of the trap and the initial well in which the condensate is prepared in a non-trivial manner. Furthermore, we show that the phenomenon of resonantly enhanced tunneling can be accompanied by a large degree of fragmentation (depletion) for the strong (weak) interaction. The many-particle position and momentum variances follow the density oscillations of the system in the asymmetric double well and bears prominent signatures of the degree of depletion or fragmentation, depending on the strength of the interactions. These quantities further exhibit a fine structure signifying a breathing-mode oscillation. Finally, a universality of fragmentation for systems made of different numbers of particles but the same interaction parameter is also found and its dependence on the asymmetry is investigated. The phenomenon is robust despite the asymmetry of the junction and admits a macroscopically-large fragmented condensate characterized by a diverging many-particle position variance. This is as far as one can get from the dynamics of the density in the junction.


Introduction
The dynamics of ultra-cold quantum gases has attracted a lot of interest since the experimental observations of Bose-Einstein condensation (BEC) [1][2][3]. The advent of advanced trapping techniques and controlling of interparticle interactions has made it possible to experimentally study several problems which were elusive until the density oscillations for the right well. Further, the repulsive inter-atomic interaction facilitates the tunneling between the two wells when the initial condensate is prepared in the left well. On the other hand, if the initial BEC is prepared in the right well, the repulsive interaction suppresses the oscillations further. For a stronger interaction, the BEC becomes fragmented with time and the degree of fragmentation is found to depend on the initial well. Moreover, the phenomenon of resonantly enhanced tunneling is found to be accompanied by a large degree of fragmentation (depletion) for the strong (weak) interaction when the initial condensate is prepared in the right well. A universality of the fragmentation dynamics is also observed, though again the degree of the universal fragmentation differs for the left and the right well. We also found prominent signatures of density oscillations as well as breathing-mode oscillations in the time evolution of the variances of the many-particle position and momentum operators. Note that for the description of the breathing mode oscillations, one needs to take into account the coupling with higher energy bands and, therefore, it is beyond the scope of Bose-Hubbard dimer.
This paper is organized as follows. In section 2, we provide a brief outline of the theoretical framework used in this work. In section 3, we present and discuss our findings. Finally, we summarize and put our concluding remarks in section 4. Physics of the one-body dynamics necessary for understanding the outcomes of this work is given in appendix A. Further details about the theoretical framework including the physical quantities used to characterize the dynamics, the in-principle numerically-exact many-body method used to solve the timedependent many-body Schrödinger equation, and discussion about numerical convergences are given in appendix B.

Theoretical framework
Here we are interested in the dynamics of a system of N interacting structureless bosons in a one-dimensional asymmetric double well which is governed by the time-dependent many-body Schrödinger equation: Here x j is the coordinate of the jth boson, is the one-body Hamiltonian containing kinetic energy T(x) and an asymmetric double-well trapping potential V T (x) terms, and W(x j −x k ) is the pairwise interaction between the jth and kth bosons. Dimensionless units are employed throughout this work.
Here the symmetry in V T (x) is broken by adding a linear slope of gradient C to the spatially symmetric doublewell. Further discussions about the system can be found in appendix A. We used dimensionless unit in this work.
The dynamics of the system is studied by exploring the time evolution of different physical quantities such as the survival probability, depletion and fragmentation, and the many-particle position and momentum variances. While some of these quantities can be studied both at the mean-field and the many-body levels, others can only be studied at the many-body level. We refer to appendix B.2 for the definitions and discussions about these quantities.
The time period of the Rabi oscillations in the double well, = p t E E Rabi 2 2 1 (E 1 and E 2 being the lowest two energy levels in the double well trap), provides a natural choice for the time scale of the dynamics and therefore, we will use the time period of Rabi oscillations t Rabi as a unit of time for the description of the dynamics in a particular asymmetric double well trap. However, as shown in appendix A, t Rabi varies with C and therefore is not suitable for comparing the dynamics in different asymmetric traps. Thus, for comparing the dynamics in different traps, we will use t 0 =100 as a unit of time.

Results
In this section, we discuss the outcome of our investigation of the dynamics of a BEC in an asymmetric double well. Specifically, we are interested to understand how the asymmetry between the two wells influences the overall dynamics of the BEC for different interaction strengths. In this work, we consider the dynamics of systems made of N=100-10 000 bosons interacting via a contact δ interaction of strength λ 0 which corresponds to the interaction parameter Λ=λ 0 (N−1). We changed the asymmetry between the two wells by varying C from a very small value to a large value for which t Rabi becomes very small and corresponds to the time period of breathing mode oscillations, see appendix A. In experiments, it is possible to vary the asymmetry of the trap by continuously tuning the energy gap between the two local minima of the double well such that once the left well is lower than the right one and in the other limit the right well is the lower one [53]. Therefore, in experiments, one can have a double well with marginally small asymmetry. For our choices of parameter C, the left well is always lower than the right one. For an asymmetric double well trap, one can prepare the initial state either in the left well V L (x) or in the right well V R (x), and then allow the system to evolve in time in the double well V T (x). Accordingly, we will study the dynamics of the system once starting from V L (x) and then from V R (x).
3.1. Quantum dynamics in an asymmetric double-well As already mentioned above, we will characterize the dynamics in an asymmetric double well trap by the time evolution of a few physical quantities such as the survival probability, depletion and fragmentation, and the many-particle position and momentum variances. The corresponding dynamics in the symmetric double well will serve as a reference for our analysis. We studied the time evolution of these quantities at the many-body levels for a weak as well as a strong interaction strength Λ. We also studied the corresponding dynamics at the mean-field level, wherever applicable, to explicitly highlight the many-body effects in the dynamics.

Survival probability
We start with the survival probability p(t) in the initial well (left well) of a symmetric double well which will serve as the reference for our subsequent analysis of the survival probabilities in the left (p L (t)) and the right (p R (t)) wells of an asymmetric double well. As discussed in appendix B.2, the survival probabilities can be studied both at the mean-field and many-body levels. Accordingly, in figure 1(a) we plot the mean-field results of p(t) for different Λ. We see that p(t) performs smooth oscillations back and forth between the two wells. For a symmetric double well, the one-body Hamiltonianˆ( ) h x is invariant under parity and therefore its eigenstates are also parity eigenstates: the ground state has even parity while the first excited state is odd. Accordingly, the superposition of these two states are localized in one or the other well. Therefore, when a one-particle state initially localized in one well is allowed to evolve in time, it keeps on tunneling back and forth between the two wells. However, in case of systems with a finite number of interacting particles like a BEC, there will be an effect of inter-particle interactions on this tunneling dynamics. Such effects are manifested through the frequency of oscillations of p(t) in figure 1(a). We observe that, as the inter-atomic interaction Λ increases, the frequency of oscillations of p(t) decreases. In the same figure, we also plot the many-body results of p(t) for N=1000, Λ=0.01, and M=2 orbitals. The complete overlap between the mean-field and the many-body results of p(t) confirms that for these parameters, the density per particle of the system and hence the survival probability can be accurately described by the mean-field theory.
Next, we consider an asymmetric double well potential with a very small asymmetry, C=0.001. Due to the presence of asymmetry, the parity symmetry ofˆ( ) h x is now lifted and therefore, the eigenstates ofˆ( ) h x are no more parity eigenstates. Accordingly, the superpositions of the first two eigenstates ofˆ( ) h x (see figure A1) are no longer well localized in one or the other well. Therefore, if a one-particle state initially localized in one well is allowed to evolve in time, it will become partially delocalized over both wells and hence there will never be full oscillations of the density of the system between the two wells.
However, for such a small asymmetry C=0.001 and a weak interaction Λ=0.01, we did not find any visible suppression of oscillations of p L (t) and p R (t) at the mean-field level (not shown here). Moreover, both p L (t) and p R (t) lie on top of each other. Therefore, we conclude that C=0.001 is too small of an asymmetry to have any visible impact on the tunneling dynamics of the system. We also repeat our calculations of p L (t) and p R (t) for a system of N=1000 interacting bosons by the MCTDHB method with M=2 orbitals and confirm that these mean-field descriptions of p L (t) and p R (t) are accurate.
Next, we enhance the asymmetry to C=0.01 keeping Λ=0.01 fixed. The mean-field p L (t) and p R (t) are shown in figure 1(b). Now, we observe the expected suppression of tunneling between the two wells. The amplitudes of oscillations of both p L and p R have decreased by nearly 40% indicating that almost 40% of the system does not tunnel out of the initial well. Moreover, though p L (t) and p R (t) practically overlap with each other, a small phase difference is found to develop with time after a few oscillations. This small phase difference is the combined effect of asymmetry and such weak interaction on the dynamics. In the same figure, we also plot the many-body p L (t) and p R (t) obtained with M=2 orbitals. That the respective mean-field and many-body curves for p L (t) and p R (t) again lie atop each other confirms that the mean-field description of the system is accurate for such a weak Λ. Thus, here we observe that even an asymmetry as small as C=0.01 has a prominent effect on the macroscopic tunneling dynamics of the system.
To systematically understand this effect of asymmetry on the macroscopic tunneling dynamics, we studied the time evolution of the survival probabilities for various C keeping Λ=0.01 fixed. We present our mean-field results in figure 1(c) for the left well and in figure 1(d) for the right well. We have checked that our mean-field results perfectly matches with the corresponding MCTDHB results obtained with M=2 orbitals. As an example, we have also shown the many-body results for C=0.25 only in figures 1(c) and (d). Again, both the many-body and mean-field results lie atop each other and hence, it reassures that the mean-filed analysis at such weak interaction is accurate even for the strong asymmetry. We observe that with increasing asymmetry C, the density oscillations are increasingly suppressed for both wells. This leads to practically self-trapping of the system in the left well for large asymmetry C0.1. However, for the right well, we observe sudden enhancement in the density oscillations for C=0. 25. With further increase of C to 0.3, the density oscillations of p R (t) is damped again though the amplitudes of oscillations are still larger than those for smaller C values. These enhanced oscillations of p R (t) for C=0.25 are the signature of resonantly enhance tunneling. As explained in appendix A, for C equal to integer multiples of 0.25 the one-body ground state energy of the right well matches with one of the one-body excited states of the left well and we refer to them as resonant values of C. For C=0.25, the one-body ground state of the right well matches with the one-body first excited state of the left well leading to the resonant tunneling. For C=0.3, there is a mismatch between the one-body ground state of the right well and the one-body first excited state of the left well leading to reduced oscillations in p R (t). However, for such small increase in C from the resonance value C=0.25, the difference between the one-body energy levels is quite small and therefore the oscillations are larger than those for other non-resonant C values. For C=0.5 (not shown), we again have large oscillations in the time evolution of p R (t) confirming that these are indeed signatures of resonant tunneling.
To further probe the effect of interaction on the tunneling dynamics between the two wells of an asymmetric double well, we next increase Λ, keeping the asymmetry C=0.01 fixed. The mean-field results of p L (t) and p R (t) for different Λ are shown in figures 1(e) and (f), respectively. We find a complementary effect of Λ on p L (t) and p R (t) at the mean-field level. While both p L (t) and p R (t) still oscillates back and forth, their amplitudes and frequencies vary with Λ in an opposite fashion. While stronger Λ facilitates oscillations of p L (t), it suppresses the oscillations of p R (t). Also, the frequencies of oscillations are found to decrease with increasing Λ for starting the dynamics from the left well, whereas when started from the right well, the frequencies increase with increasing Λ.
Qualitatively, the repulsive interaction pushes up the energy of the BEC with respect to the barrier height. For a sufficiently high barrier, the energy levels of the ground state and the first excited state of an asymmetric double well approximately coincide with the ground states of the left (lower) and the right (upper) wells, respectively. Therefore, when the initial state is prepared in the right (upper) well, the energy of the initial condensate becomes closer to the energy level of the first excited state of the asymmetric double well with increasing repulsive interaction Λ. Thus, the system tends more to remain in that state and the tunneling is more and more suppressed with increasing Λ. On the other hand, when the initial condensate is prepared in the left (lower) well, its energy is pushed away from the ground state of the asymmetric double well by the repulsive interaction and therefore, it becomes more prone to tunneling with increasing Λ.
At stronger Λ keeping N fixed, a mean-field theory may not be sufficient to describe the system. Accordingly, we again refer to the symmetric double well case. Figure 2(a) exhibits the time evolution of p(t) for a symmetric double well calculated by MCTDHB with M=2 orbitals. In the inset, the corresponding mean-field result of p (t) is provided for comparison. We clearly see that, contrary to the mean-field result, the many-body result for p (t) exhibits a collapse of the oscillations. Such collapse of oscillations in the many-body dynamics of p(t) has been reported earlier also [49,50]. For a symmetric double well, such collapse of oscillations has recently been observed [71] in the oscillations of the population imbalance which is directly related to the survival probability, see appendix B.2. This makes a many-body calculation necessary for Λ0.1 for a system of N=1000 bosons. Therefore, next, we calculate the p L (t) and p R (t) by MCTDHB method with M=2 orbitals for the same parameters as in figures 1(c) and (d) for N=1000 bosons. As an example, here we present the many-body results only for Λ=0.1 in figure 2(b). The collapse of oscillations for both p L (t) and p R (t) can be seen on top of the overall mean-field effects described above. However, the collapse time differs: While the collapse for p L (t) is quicker compared to the symmetric double well, it is delayed for p R (t). We further found that, with an increase in Λ, the collapse is quicker for both the p L (t) and p R (t) for a fixed asymmetry C. On the other hand, for a fixed Λ, the collapse of both p L (t) and p R (t) is deferred with increasing C in terms of the number of Rabi cycles (recall that t Rabi depends on C).
Next, to explore how the resonantly enhanced tunneling is affected in the strong interaction regime, we have studied the time evolution of p L (t) and p R (t) by MCTDHB for various C and Λ=0.1 and presented our results in figures 2(c) and (d), respectively. We again see that with increasing C the density oscillations are suppressed right from the beginning for both wells, except for the values of C where there is resonant tunneling for p R (t). However, now we see that even for resonant values of C, though there are large density oscillations at the beginning, they are gradually damped out with time for such large interaction. This marks the transition of the resonant tunneling phenomenon into the many-body regime. We have checked that such damping of oscillations also happens for the stronger Λ.

Depletion and fragmentation
Having seen that a many-body calculation exhibits new features already for the time evolution of the survival probabilities at stronger Λ, next we would like to examine the time development of the depletion and fragmentation f (depending on Λ) of the condensate which is a purely many-body quantity. We found the BEC to become depleted with time for weak interactions such as Λ=0.01. In figure 3 the time development of the depletion f of a BEC made of N=1000 bosons is shown for different asymmetries C. We observe that f is very small and therefore the system is practically condensed. Explicitly, n 1 >999.999 for the times shown in figure 3. Thus for such interaction strengths, the density per particle of the system is accurately described by the meanfield theory. Even then, time development of f exhibits some interesting features. First, the depletion is found to be maximum for the symmetric double well (C = 0) and gradually decreases with increasing C, except for resonant values of C. Further, we observe an interesting difference between the time development of f for the left and right wells. Whereas for very small C the respective time evolutions of f are essentially same, the difference starts to develop with C and, while for C0.005, the depletion for the left well is larger than that for the right well, the situation reverses for C=0.01. For stronger values of C, f is negligibly small except for the resonant  values of C when the depletion is sharply enhanced for the right well. With further increase in C, as one moves away from the resonance the depletion for the right well starts to decrease again.
Next we consider a stronger interaction Λ=0.1. In figures 4(a) and (b), we plot the natural occupations n N j for a system of N=1000 interacting bosons as a function of time for different asymmetries C. We also plot the corresponding results for the symmetric double well (C = 0) in both panels as a reference. The results presented here are obtained with M=2 orbitals. For all cases, we observe that starting from » 1 n N 1 , the occupation in the first orbital n N 1 gradually decreases with time. Simultaneously, the occupation in the second orbital n N 2 slowly increases with time starting from a negligibly small value. Thus, with time the condensate becomes fragmented with a fragmentation fraction = f n N 2 . Finally, as the density oscillations collapse (see figure 2(b)), f reaches a plateau at f=f col . Moreover, we see small oscillations in f prior to attaining the plateau. Such oscillations are the signatures in f of the time-dependent density oscillations. As the density oscillations collapse by the time f reaches the plateau, the oscillations in f are also heavily damped and remain so at the plateau.
Comparing the results for different asymmetries C, we find that both the growth rate of f and f col depend on C. However, there is a crucial difference between the cases when the initial BEC state is prepared in the left (figure 4(a)) and the right wells (figure 4(b)). For the left well, f col first increases with increasing C and the condensate becomes more fragmented with reference to the symmetric double well until C=0.002. With further increase in C, f col decreases and the condensate becomes less fragmented. For stronger asymmetry C0.1, the the occupation of the first orbital n N 1 is nearly 100% and the system is essentially condensed. On the other hand, for the right well, f col is found to decrease monotonically with increasing C, as far as C0.005. For stronger asymmetry C=0.1, the system remains essentially condensed for a long time also for the right well. However, for resonant values of C such as C=0.25, the system becomes fragmented again with the degree of fragmentation f being larger than that in the symmetric double well. Even for a small increase of asymmetry C from the resonant values, f decreases sharply as shown for C=0.3 in figure 4 We may understand the findings for small asymmetry C qualitatively by treating the small asymmetry as a perturbation. In [49] for each eigenstate ñ |E n of the Bose-Hubbard dimer, its fragmentation f n as a function of the eigenstate energy per particle E n /N has been discussed. It has been shown both analytically and numerically that f n first increases with E n /N, reaches a maximum of 50%, and then decreases with further increase of E n /N. For a small perturbation, such qualitative functional dependence is expected to remain valid. Comparing the results in figures 4 with 3 of [49], we can infer that, for the parameters used in this work, the initial state for the symmetric double well lies on the upper part of the right-hand-branch of the f n versus E n /N curve (figure 3 of [49]). Now, the introduction of a small asymmetry C pulls down the energy for the left well and pushes up for the energy for the right well. Accordingly, f col for the left well initially increases with C and then, with further increase of C and consequent decrease of the initial state energy, it crosses to the left branch of the f n versus E n /N curve (see figure 3 of [49]) and starts to decrease. On the other hand, with the increase in C the initial eigenstate energy for the right well monotonically increases resulting in a monotonic decrease of f col (for not too large C). One would need to go beyond such a perturbation-based analysis to understand the behavior of fragmentation for C>0.005 in the right well, see figure 4(b). 3.1.3. Many-particle position and momentum variance Next, we consider the time evolution of the many-particle position and momentum variances. These quantities characterize the fluctuations in the particles' positions and momenta in the junction. Although not easy to measure, they are fundamental quantum mechanical observables. Since these quantities depend on the actual number of depleted or fragmented atoms, it is expected that prominent signatures of the depletion and fragmentation of the condensate would show up in these variables. For Λ=0.01, it has been shown above that the system remains practically condensed for a long time (see above) and therefore, its out-of-equilibrium dynamics should be adequately described by the mean-field theory. So, first we study the time evolution of the many-particle position variance DN X 1 2 at the mean-field level. We present our results for different asymmetries C and a fixed Λ=0.01 in figures 5(a) and (b) for preparing the initial BEC state in V L (x) and V R (x), respectively. For comparison, we also plot DN X 1 2 for the symmetric double well in both panels.
We observe that for both V L (x) and V R (x) of the asymmetric double well trap, DN X 1 2 oscillates with a frequency which equals to the Rabi frequency. This is in contrast to the case of the symmetric double well in which DN X 1 2 oscillates with a frequency equal to twice the Rabi frequency. This is due to the incomplete tunneling between the two wells of the asymmetric double well trap, and that there is always a remnant in the each well which is further manifested in the irregular peaks of oscillations of DN X 1 2 . We observe that with increase in C starting from the symmetric double well (C = 0), the peaks of the oscillations of DN X 1 2 first split into two sub-peaks which gradually turn into broad peaks for C=0.01. For resonant values of C such as C=0.25, while the sub-peaks of the oscillations of DN X 1 2 become sharper for the right well, DN X 1 2 for the left well practically remains constant approximately at 0.5. This is consistent with earlier observation that for large asymmetries C0.1, the system is self-trapped for the left well while for the right well it exhibits large density oscillations. Further, we observe high-frequency small-amplitude oscillations on top of the peaks of the large-amplitude oscillations. Such high-frequency oscillations are because of the breathing-mode oscillations of the system in the asymmetric double well and can be seen more vividly in the many-particle momentum variance (see below). Also, the minima of the oscillations are slightly higher than 0.5 for all times t>0. Moreover, comparing the panels (a) and (b), we see that the peak values of the oscillations are slightly higher for the right well V R (x). All of these quantify the fluctuations in the particles' positions in the asymmetric double well at the mean-field level.
As discussed in section 3.1, the many-particle position variance can deviate from their corresponding meanfield results even when the mean-field theory is expected to accurately describe the density per particle of the system. So, we now study the time evolution of the many-particle position variance DN X 1 2 at the many-body level. For all cases, DN X 1 2 is found to grow in an oscillatory manner. For the symmetric double well (C = 0), the maxima of DN X 1 2 grow approximately quadratically, also see [52]. This growth is slower for an asymmetric double well, where the growth rate decreases with increasing C, except for the resonant values. This is consistent with our earlier observation that the depletion of the condensate is maximal in a symmetric double well. As observed at the mean-field level, here also, the oscillations of DN X 1 2 are irregular in nature. However, now the two sub-peaks are of unequal heights and the difference between them grows with time t for both V L (x) and V R (x). Comparison between the left (panel (c)) and the right wells (panel (d)) shows that, while the higher sub-peaks is on the left side for V L (x), it appears on the right side for V R (x). Further, while for C=0.005 the maximal values for the DN X 1 2 in left well are larger than those in the right well, the situation reverses for C=0.01. This is again consistent with our earlier observation (in figure 3) that the system in the left well is more depleted until C=0.005, whereas the system in the right well is more depleted for C=0.01. For the resonant values of C such as C=0.25, the growth rate of the maxima of DN X 1 2 is much higher for the right well than that for non-resonant (and smaller) values of C. For the left well, DN X 1 2 for such large C is again practically constant approximately at 0.5. This is consistent with the fact that the system becomes more depleted for the right well for resonant values of C.
Next, in figure 6, we plot the many-particle momentum variance DN P 1 2 of the system for starting the dynamics from V L (x) and V R (x), respectively. We studied DN P 1 2 both at the mean-field and the many-body levels.
In figures 6(a) and (b) we present the mean-field results of DN P 1 2 for the left and the right wells, respectively. In each panel, we also plot DN P 1 2 for the symmetric double well for comparison. Generally, we observe two oscillations associated with the time evolution of DN P 1 2 : the first, with a larger amplitude and frequency equal to twice the Rabi frequency and, the second, with a smaller amplitude but a higher frequency. The first one is a manifestation of the density oscillations, whereas the second one is due to the breathing oscillations of the system. However, while in the symmetric double well the amplitude of the breathing mode oscillations are larger than those of the density oscillations, the situation is reversed in the asymmetric double well with small C. For large C0.1, only the breathing oscillations are prominent in the time evolution of DN P 1 2 for the left well as the system is practically self-trapped in the left well. Further, for resonant values of C such as C=0.25, only the density oscillations are prominent in the time evolution of DN P 1 2 for the right well. This is because of the much larger density oscillations due to the resonantly enhanced tunneling at such C values. A closer examination of the high frequency breathing mode oscillations suggests that these may arise due to the transition of two bosons from the lowest energy band to the second band or one boson from the lowest band to the third band. An analysis by a linear-response theory in the line of [89] is required to attribute such highfrequency oscillations to a particular transitions unambiguously and accurately. In any case, it can be safely said that one needs to consider higher bands to take into account such high-frequency breathing mode oscillations and this, of course, is beyond the scope of the standard Bose-Hubbard model. Moreover, though DN P 1 2 starts from the same value for both the left and right well, there is a Π-phase difference between the oscillations for the left well and the right well. Explicitly, the momentum variance first decreases when starting from the lower (left) well, whereas it first increases when starting from the higher (right) well. This can be understood from an energetic point of view. The BEC tunneling from the lower to upper well initially loses kinetic energy (momentum) and gains kinetic energy when tunneling from the higher to lower well. The momentum variance behaves accordingly. Finally, we show the corresponding MCTDHB results of DN P 1 2 with M=2 orbitals in figures 6(c) and (d) for the left and the right wells, respectively. We find that the MCTDHB dynamics of DN P 1 2 is similar to the corresponding mean-field dynamics. Actually, the many-particle momentum variance depends on the derivatives of the orbitals. For a weak asymmetry and a weakly interacting system, the shape of the orbitals deviate only slightly from their corresponding (non-interacting and) mean-field shape, and this leads to even smaller derivatives thereby producing practically same DN P 1 2 both at the mean-field and the many-body levels.

Universality of the fragmentation dynamics in an asymmetric double well
A unique many-body feature predicted in the dynamics of BECs in a symmetric double well is the universality of the degree of fragmentation with respect to N for a fixed Λ [49]. It was first established by solving the many-body Schrödinger equation and then using the Bose-Hubbard dimer, it was also shown that the universality of fragmentation in a symmetric double well is a general many-body phenomenon [49]. Also, in the previous subsection, we have already found a significant effect of the asymmetry of the trap on the time evolution of the survival probabilities, fragmentation, and the many-particle position and momentum variances of BEC in an asymmetric double well trap. Naturally, questions arise if the universality of the fragmentation exists in an asymmetric trap and if so, how it is affected by the asymmetry of the trap. Once again we start with the corresponding symmetric double well as a reference. In figure 7(a) we have plotted the natural occupations for different N keeping Λ fixed. As discussed above, we see that initially only one natural orbital is occupied with » 1 n N 1 and negligibly small f for all cases shown in figure 7. However, with time the second natural orbital starts to be occupied, the system becomes fragmented and, during the collapse of the density oscillations, the occupations of the natural orbitals reach the same plateau for different numbers of bosons N keeping Λ fixed. The values at the plateau are about = 60% n N 1 and = 40% n N 2 , respectively. Hence for all cases, after the collapse of the density oscillations the system becomes f col ≈40% fragmented irrespective of N, showing a universal fragmentation dynamics [49].
Next, we consider an asymmetric double well with a very small asymmetry C=0.001. Figure 7(b) shows the results for V L (x) and figure 7(c) for V R (x). For both wells, qualitatively, we see the same dynamics as in the symmetric double well. We observe that following an oscillatory growth, f reaches the same plateau f col irrespective of the number of particles N for a fixed Λ. Therefore, the universality of the fragmentation dynamics also persists in an asymmetric double well. However, quantitatively f col for the left well differs from that for the right, f col ≈45% versus f col ≈35%, respectively. Here, an interesting point is that f col for the symmetric well is actually the mean of f col for the two wells of the asymmetric double well. As shown in [49], the fragmentation depends on the energy per particle. Since we have introduced the asymmetry by adding a linear slope of a fixed gradient, it pushes up the right well by about the same amount as it pulls down the left well. Therefore, the changes in f col for both wells are expected to be similar but in opposite directions, leading to the above relation between the fragmentation values.
As discussed earlier, the many-particle position variance DN X 1 2 bears prominent signatures of the fragmentation. Hence, next, we study the time evolution of the many-particle position variance DN X 1 2 , both at the mean-field and the many-body levels, to explore the possible manifestation of the universality of the fragmentation dynamics. In figure 8 we plot the MCTDHB results with M=2 orbitals of DN X 1 2 , for different N but the same Λ, as a function of time for starting the dynamics from both the left (panel (a)) and the right (panel (b)) wells. We also plot the corresponding mean-field results in both panels for comparison. In the mean-field theory, there is only one parameter Λ and therefore, we have only one curve for the time development of DN X 1 2 for a particular Λ irrespective of N. On the other hand, at the many-body level, we find different time for N=1000 and 10000, respectively. These observations can be understood as follows. In figure 7, we have seen that f col for the left well is only about 10% higher than that of the right well, for all N. Naturally, the actual occupation numbers n 2 are of the same order of magnitude for both wells for all N. Since DN X 1 2 depends on the actual value of n 2 [90], its saturation values for a particular N are of the same order of magnitude (as a power of 10) for both the wells. Similarly, due to the universality of fragmentation dynamics, f col corresponding to different N and same Λ have the same value for a particular well. Therefore, the actual number of fragmented atoms n 2 increases by a factor of N also exhibit an oscillatory growth followed by an equilibration after the collapse of the density oscillations. However, the important point is that the curves for different N, keeping Λ same, saturate to the same mean value about which DN X 1 2 2 keeps on oscillating. This is the signature of the universal fragmentation dynamics. Therefore, the universality of fragmentation is a quite robust many-body phenomena and its signature appears in all many-body quantities that depend on the occupation numbers of the natural orbitals.

Summary and concluding remarks
Summarizing, we have examined how the BJJ dynamics is affected by the gradual loss of symmetry of the confining double well potential for different interaction Λ. In an asymmetric double well, the two wells are no longer equivalent, say the left well is deeper. Therefore, we have studied the dynamics by initially preparing the condensate in both the left and right wells. We have analyzed the dynamics by examining the time evolution of the survival probability, depletion or fragmentation, and the many-particle position and momentum variances.
We find that the impact of the asymmetry of the trap depends on the interaction Λ and the initial well. Overall, the tunneling between the two wells becomes increasingly suppressed with growing asymmetry. Moreover, the repulsive inter-atomic interaction facilitates the tunneling between the two wells when BEC is initially in the left well whereas the tunneling is further suppressed for starting the dynamics from the right well. However, for the resonant values of C, unusually large density oscillations due to the resonantly enhanced tunneling are observed for preparing the initial BEC state in the right well. Further, while we observed persistent large oscillations in p R (t) due to the resonantly enhanced tunneling for weak interaction, for large Λ such oscillations are gradually damped out with time. We point out that such damping of the oscillations is purely many-body in nature.
For a sufficiently strong interaction Λ, the condensate becomes fragmented with time and the degree of fragmentation f depends on the asymmetry C and the initial well. For preparing the initial BEC in the right well, f decreases with increasing asymmetry C between the two wells of an asymmetric double well (except for resonant values of C) and for a sufficiently strong non-resonant values of C, the system remains essentially condensed for a long time. For resonant values of C, we observed a sharp increase in f and the system becomes more fragmented than in the symmetric double well. In case of the left well, the system first becomes more fragmented for very small degrees of asymmetry, then with further increase in asymmetry, it again become less fragmented than in a symmetric double well and finally for sufficiently large asymmetries, it does not fragment at all. Further, that f reaches the saturation value around the same time when the density oscillations are damped out indicates that fragmentation and inhibition of tunneling goes hand in hand. This holds for both the symmetric as well as the asymmetric double well.
The time evolution of the DN X 1 2 bears prominent signature of the depletion of the system and deviates from its corresponding mean-field dynamics even for weak interaction Λ. In an asymmetric double well, both the frequencies and the amplitudes of the oscillations of DN X 1 2 are found to be affected by the asymmetry. In particular, the growth rate of the maxima of DN X 1 2 decreases with increasing C for both wells. However, this growth rate for the right well increases sharply for the resonant values of C. The dynamics of DN X 1 2 in an asymmetric double well trap bears signatures of the breathing-mode oscillations in addition to the density oscillations. The signatures of the breathing-mode oscillations are more prominent in the time evolution of the many-particle variance DN P 1 2 . In the time evolution of DN P 1 2 , while breathing-mode oscillations are more prominent than the density oscillations for the symmetric double well, both are distinctly visible in case of the asymmetric double well. Further for the resonant values of C, while only breathing oscillations are prominent for the left well, in case of right well, signatures of the density oscillations are more prominent due to the resonance enhanced tunneling. Since breathing-mode oscillations arise from coupling to higher energy bands, such features are beyond the scope of the Bose-Hubbard dimer. An important observation of our study is the universal fragmentation dynamics of asymmetric BJJ. However, the degree of universal fragmentation for BECs consisting of different N corresponding to the same Λ, depends on the initial well. Universality of fragmentation is found to manifest in the same mean saturation value of the DN X 1 2 2 for different N corresponding to the same Λ at the many-body level. This means that the fluctuations of the positions of the particles in the junction show a universal behavior.
In conclusion, our work has revealed that the many-body dynamics of BEC in an asymmetric double well depend on the asymmetry C of the trap and the interaction strength Λ in a complicated manner and different quantities characterizing the dynamics responds differently to the many-body effects. Transition of the resonantly enhanced tunneling into the many-body regime, enhanced degrees of fragmentation and depletion (depending on Λ) for resonant values of C, and a universality of fragmentation are some of the important manybody effects found in this work.
Performance Computing Center Stuttgart (HLRS) is gratefully acknowledged. We acknowledge valuable comments by the referees. SKH also gratefully acknowledges the continuous hospitality at the Lewiner Institute for Theoretical Physics (LITP), Department of Physics, Technion-Israel Institute of Technology.

Appendix A. One-body dynamics
The one-body dynamics is completely determined by the trapping potential V T (x). As mentioned in the main text, the asymmetric double well V T (x) is constructed by adding a linear slope of gradient C to the symmetric double well which itself is obtained by fusing two slightly shifted harmonic potential The symmetric double-well part is taken from [50]. This will allow us to relate and compare results in the asymmetric junctions to that in the symmetric one. In this work, as mentioned in the main text, C is varied from a very small value to a large value such that t Rabi becomes very small and corresponds to the time period of breathing mode oscillations (see below). The shape of an asymmetric double well for C=0.01 used in this work along with its first few energy levels E n and eigenstates j n (x) are shown in figure A1(a). One can see that it is hardly distinguishable from the symmetric double well. Also, even for such a small asymmetry C, the superposition of the first two eigenstates are not completely localized in one or the other well, thereby affecting the density oscillations between the two wells. Moreover, the spacing between the two successive energy levels increases as one goes up the spectrum: while the lowest two energy levels lie very close to each other and form the lowest energy band, the higher energy levels from E 5 onward are practically unaffected by the barrier between the two wells and form an almost uniform spectrum. Therefore, the dynamics of the system in such an asymmetric double well is primarily controlled by the lowest energy band. However, with increasing C, the spectrum starts to be affected more prominently by the asymmetry and the higher energy levels begin to play more important role in the dynamics.
To highlight the point further, we compute the ratio of the inter-band spacing to the intra-band spacing of the lowest band, viz. = D E E 52 21 are quite large and of the order of their values for the C=0. Therefore, in the regime of our interest, the lowest energy band is expected to play the lead role in the dynamics of the system. However, the next nearest band may influence the dynamics by giving rise to the breathing mode oscillations on top of the Rabi oscillations controlled primarily by the lowest band. Also, in our present study, the inter-atomic interaction W(x j −x k ) may lead to a coupling with the higher energy levels. Therefore, even for such a small asymmetry, it is necessary to effectively take all bands into account. Only then, one can be sure that the lowest band is the dominant one.
The time period of the Rabi oscillations in the double well, = , provides a natural choice for the time scale of the dynamics. t Rabi as a function of the asymmetry C is shown in figure A1(c) with the region of our interest being highlighted in the inset. We note that t Rabi also decreases exponentially with C. Actually, for a small asymmetry C, the ground state and the first excited state in the asymmetric double well are delocalized, and t Rabi gives the time period of Rabi oscillations in the double well. However, for large C, the barrier becomes very high and there is no tunneling back and forth between the two wells leading to 'self-trapping' of the one-particle system in the initial well. Then E 1 and E 2 become the lowest two energy levels in the lower well V L (x) and is associated with the time period of breathing mode oscillations.
However, for certain values of C, it can happen that the one-body ground state energy of the right (upper) well matches with one of the one-body energy levels of the left (lower) well. In that case, there will be a pronounced tunneling from the ground state of the right well although the corresponding tunneling from the ground state of the left well will still be practically blocked. Such conditions are known as the resonantly enhanced tunneling [60,61]. In our case, we can see that the local minima of the two wells are vertically shifted by 4C. Therefore, the one-body spectrum of the individual wells will also be separated vertically by 4C. So, for the specific C values such that 4C is integer, the energy of the ground state of the right well matches with one of the excited states of the left well leading to resonantly enhanced tunneling. We refer to these values of C as the resonant values. For C=0.25, the one-body ground state energy of the right well matches with the one-body first excited state of the left well while for C=0.5 the one-body ground state of the right well matches with the one-body second excited state of the left well leading to the resonantly enhanced tunneling for the right well.

B.1. The MCTDHB
The time-dependent many-boson Schrödinger equation (1) cannot be solved exactly (analytically), except for a few specific cases only, see, e.g. [91]. Hence, to solve equation (1) in-principle numerically exactly, Below we briefly describe the basic idea behind the method.
In MCTDHB, the ansatz for solving equation (1)   preserve the total number of bosons N. For an exact theory, M should be infinitely large. However, for numerical computations one has to truncate the series at a finite M. In actual calculations, we keep on increasing M until we reach the convergence with respect to M and thereby we obtain a numerically-exact result. In the context of bosons in a double-well, the latter implies that the MCTDHB theory effectively takes all required bands into account. Here we would like to point out that for M=1, the ansatz equation (B1) gives back the ansatz for the Gross-Pitaevskii theory [73].
Therefore, solving for the time-dependent wavefunction Ψ(t) boils down to the determination of the timedependent coefficients  { ( )} C t n and the time-dependent orbitals {f k (x, t)}. Employing the usual Lagrangian formulation of the time-dependent variational principle [92,93] subject to the orthonormality between the orbitals, the working equations of the MCTDHB are obtained as follows nn . It is convenient to define the different quantities of interest in terms of the one-body and the two-body reduced density matrices [94][95][96][97] instead of the full many-body wavefunction. Given the normalized manybody wavefunction Ψ(t), the reduced one-body density matrix can be calculated as (M is the number of single particle orbitals used to construct the many-boson wavefunction, see section B.1). If only one macroscopic eigenvalue  » ( ) ( ) n t N 1 exists, the system is condensed [98] whereas if there are more than one macroscopic eigenvalues, the BEC is said to be fragmented [99][100][101][102]. The diagonal of the r ¢ ( | ) ( ) x x t ; 1 1 1 gives the density of the system r r º ¢ = ( ) ( | ) ( ) x t x x x t ; ; 1 . Similarly, the two-body density can be calculated as Therefore, the matrix elements of the two-body reduced density matrix are given by r = áY Yñ where b k and † b k are the bosonic annihilation and creation operators, respectively.

B.2. Physical quantities
In this work, we will study the dynamics of the system by exploring the time evolution of different physical quantities defined as follows. While some of these quantities can be studied both at the mean-field and the manybody levels, others can only be studied at the many-body level.
(a) Survival probability. Survival probability measures the fraction of the system left in the initial well at any instant during the tunneling dynamics in a double well. In a double well, assuming that the BEC is initially in the left well, it can be defined as x t N 0 , ρ(x, t) being the density of the system at time t. Another common quantity to characterize the density oscillations in experiments is the population imbalance [71] which is defined as = -( ) n t ; N N N L R N L,R being the atom numbers in the left (L) and the right (R) well at time t, respectively. Note that both p(t) and n(t) are closely related and have similar qualitative features as they both characterize the same density oscillations. However, there are quantitative differences. While p(t) can oscillate about a mean value 0.5 between a maximum 1 and a minimum 0, n(t) can oscillate about a mean-value 0. Also in case of oscillations of n(t), the peak values depend on the experimental set up and can be controlled by varying the vertical shift between the two wells, see [71]. Therefore, the amplitudes of oscillations are different for p(t) and n(t), though the frequencies are in-principle the same as the frequencies depend on the tunneling rate between the two wells. In the dynamics of BEC in an asymmetric double well following a trapping quench from a harmonic well to an asymmetric double well at t=0, we can prepare the initial BEC state either in the left well (L) or in the right well (R). Accordingly, we can calculate two types of survival probabilities. For starting with the initial BEC state in the left well, left well [p L (t)] as we can define the survival probability in the where ρ L (x; t) is the density when the initial BEC state is prepared in the left well. Similarly, when the initial state is prepared in the right well, the survival probability in the right well [p R (t)] can be defined as is the density when the initial condensate is in the right well. For a symmetric double well, both ρ R (t) and ρ L (t) are equivalent and therefore we have only a single type of survival probability p(t), i.e.
. Since the density can be studied both at the mean-field and the many-body levels, survival probabilities can also be calculated both at the mean-field and the many-body levels.
(b) Depletion and fragmentation. As discussed above, when the system is condensed and the sum over all the microscopic fractions of occupations in the higher orbitals = å = f j M n N 2 j (M is the number of orbital, see above) is known as the depletion per particle. On the other hand for a fragmented system, the macroscopic occupation of a higher natural orbital, viz , is called fragmentation. From the definition, it is clear that one needs more than one orbital to study the depletion and fragmentation and hence, these quantities can only be studied by at least a two-orbital (many-body) theory [103,104] and preferably a multi-orbital (many-body) theory [51,97,105]. We remark that the depletion of a BEC is usually small and may not have a prominent effect on the density per particle and energy per particle which, in effect, can be accurately described by a mean-field theory. However, fragmentation can have a dominant effect on the energy per particle and the density per particle of the system. Moreover, though the depletion and the fragmentation are physically different quantities and appear under different conditions, for a two-orbital theory they have the same mathematical expression. Accordingly, for the computation with M=2 orbitals only we will refer to both of them by f, see section B.1 below.
(c) Many-particle position and momentum variance. The quantum variance of an observable is a fundamental quantity in quantum mechanics due to its connection with the uncertainty principle. It gives a measure of the quantum resolution with which an observable can be measured. For any many-body operator = å =ˆ( ) A a x j N j 1 whereˆ( ) a x j is a Hermitian operator and local in position space, the variance per particle D ( ) t N A 1 2 [52,89,90,106] is given by , takes into account all other contributions to the many-particle variance. Similar expressions hold for operators which are local in momentum space. We point out that one can, in principle, study the variance of any operator at the mean-field level by substituting the many-body wavefunction Ψ(t) by the corresponding mean-field wavefunction. However, in the mean-field theory, only D ( ) t a,density 2 has a non-zero contribution while D ( ) t a,MB 2 is identically equal to zero. Therefore, even for the interaction strengths for which the mean-field theory is expected to accurately describe the density per particle of the system, the many-body variance can deviate from its mean-field result. Accordingly, in this work we will consider the variances of the many-particle position and momentum operators at the many-body level only.

B.3. Numerical computations and their convergence
Here we discuss the details of our numerical computations by MCTDHB.A parallel version of MCTDHB has been implemented using a novel mapping technique [107,108]. We mention that by propagating in imaginary time the MCTDHB equations also allow one to determine the ground and excited state of interacting manyboson systems, see [51,77]. In our present work we have performed all computations with M=2 time-adaptive orbitals. By repeating our computations with M=4, 6, and 8 orbitals the results have been verified and found to be highly accurate for the quantities and propagation times considered here.
In our numerical calculations, the many-body Hamiltonian is represented by 128 exponential discretevariable-representation grid points (using a fast Fourier transformation routine) in a box size [−10,10). We obtain the initial state for the time propagation, the many-body ground state of the BEC either in the left well or in the right well, by propagating the MCTDHB equations of motion (equation (B2)) in imaginary time [51,77]. For our numerical computations, we use the numerical implementation in [107,108]. We keep on repeating the computation with increasing M until convergence is reached and thereby obtain the numerically accurate results.
Below we demonstrate the numerical convergence of the many-particle position and momentum variances. We already discussed in the text that the variance of any quantum operator is much more sensitive to the manybody effects compared to the oscillations in the survival probabilities and the fragmentation f. Actually, it is seen that the convergence of the momentum variance requires more numerical resources than the convergence of the position variance. Therefore, demonstration of convergence of the position and momentum variances will automatically imply the convergence of the survival probabilities and the fragmentation f with respect to M.
In an asymmetric double well, the two wells are not equivalent and therefore, in the text we discussed the dynamics of the system separately for preparing the initial BEC in the left well and in the right well. Accordingly here also, we discuss the convergence for both the cases separately, first when the initial BEC is prepared in the left well and then when the condensate is initially in the right well. We consider a system of N=10 interacting bosons in an asymmetric double well with the asymmetry C=0.01. Since in the limit  ¥ N , keeping Λ fixed, the many-body effects diminish and the density per particle of the system converge to its corresponding mean- Figure B1. Convergence of variances of the many-particle position and momentum operators with respect to the orbital number M for a system of N=10 interacting bosons and Λ=0.01 in the asymmetric double well trap of asymmetry C=0.01. (a) Timeevolution of the many-particle position variance DN X field values [109][110][111][112], the convergence of the quantities for higher N values considered in the text (but same Λ values considered here) are actually better than what is shown below. We first consider Λ=0.01. In figures B1(a) and (b) we plot DN X 1 2 computed with M=2, 4, 6, and 8 for starting the dynamics from the left and right well, respectively. We see that, as discussed in the main text, DN X 1 2 exhibits a slow oscillatory growth for both wells. Furthermore, we see that there is an overall oscillation of DN X 1 2 Figure B2. Convergence of variances of many-particle position and momentum operators with respect to the orbital number M for a system of N=10 interacting bosons and Λ=0.1 in the asymmetric double well trap of asymmetry C=0.01. (a) Time-evolution of the many-particle position variance DN X  in time with a frequency equal to twice the Rabi frequency. Also, on top of the peaks of these oscillations, there is another oscillation with a higher frequency but smaller amplitude. The origin of these oscillations are discussed in the main text, see section 3.1.3. Here, we observe that, for both cases, the results for M=2, 4, 6, and 8 are in very good agreement with each other, such that not only the overall oscillations of DN X 1 2 due to the density oscillations but also the small amplitude high-frequency oscillations on top of the peaks of the first ones are accurately described with M=2.
To demonstrate the convergence of DN P 1 2 , in figures B1(c) and (d), we plot the DN P 1 2 computed with M=2, 4, 6, and 8 orbitals for preparing the condensate initially in the left and the right well, respectively. We clearly see that there are two oscillations associated with the time evolution of DN P 1 2 on top of one another, one with a larger amplitude and frequency equal to twice the Rabi frequency and the other one with a smaller amplitude but higher frequency. Once again, we see that the results of DN P 1 2 for different M practically overlap with each other, and the M=2 is sufficient to describe both oscillations accurately. Next, we consider the convergence for the stronger interaction Λ=0.1. We plot DN X Finally, to explicitly show that the convergences of other quantities are also achieved with the same M as the variances, for the same system parameters, below as an example, we consider the convergence of the natural occupation numbers n N j . In figure B3 we plot n N j for a system of N=100 and Λ=0.1 computed with M=2 and 4 orbitals for preparing the initial condensate in the left (panel (a)) and right (panel (b)) wells of an asymmetric double well with C=0.001. We found that the computation with M=4 reproduces the same f as with M=2 for both initial wells: the curves for the two largest occupation numbers show that convergences improve with N keeping Λ fixed for a repulsive interaction. Finally and importantly, it also demonstrates that the universality of fragmentation is a robust many-body phenomenon and does not fissile out by using larger numbers M of orbitals.