Longitudinal and transversal resonant tunneling of interacting bosons in a two-dimensional Josephson junction

We unravel the out-of-equilibrium quantum dynamics of a few interacting bosonic clouds in a two-dimensional asymmetric double-well potential at the resonant tunneling scenario. At the single-particle level of resonant tunneling, particles tunnel under the barrier from, typically, the ground-state in the left well to an excited state in the right well, i.e., states of different shapes and properties are coupled when their one-particle energies coincide. In two spatial dimensions, two types of resonant tunneling processes are possible, to which we refer to as longitudinal and transversal resonant tunneling. Longitudinal resonant tunneling implies that the state in the right well is longitudinally-excited with respect to the state in the left well, whereas transversal resonant tunneling implies that the former is transversely-excited with respect to the latter. We show that interaction between bosons makes resonant tunneling phenomena in two spatial dimensions profoundly rich, and analyze these phenomena in terms of the loss of coherence of the junction and development of fragmentation, and coupling between transverse and longitudinal degrees-of-freedom and excitations. To this end, a detailed analysis of the tunneling dynamics is performed by exploring the time evolution of a few physical quantities, namely, the survival probability, occupation numbers of the reduced one-particle density matrix, and the many-particle position, momentum, and angular-momentum variances. To accurately calculate these physical quantities from the time-dependent many-boson wavefunction, we apply a well-established many-body method, the multiconfigurational time-dependent Hartree for bosons (MCTDHB), which incorporates quantum correlations exhaustively. By comparing the survival probabilities and variances at the mean-field and many-body levels of theory and investigating the development of fragmentation, we identify the detailed mechanisms of many-body longitudinal and transversal resonant tunneling in two dimensional asymmetric double-wells. In particular, we find that the position and momentum variances along the transversal direction are almost negligible at the longitudinal resonant tunneling, whereas they are substantial at the transversal resonant tunneling which is caused by the combination of the density and breathing mode oscillations. We show that the width of the interparticle interaction potential does not affect the qualitative physics of resonant tunneling dynamics, both at the mean-field and many-body levels. In general, we characterize the impact of the transversal and longitudinal degrees-of-freedom in the many-boson tunneling dynamics at the resonant tunneling scenarios.

In general, we characterize the impact of the transversal and longitudinal degrees-of-freedom in the many-boson tunneling dynamics at the resonant tunneling scenarios. * abhowmik@campus.haifa.ac.il

I. INTRODUCTION
Tunneling is a purely quantum mechanical phenomenon which takes place in classicallyforbidden region, originally intended to account for α-decay, fusion, and fission in nuclear physics [1,2]. Apart from nuclear physics, tunneling occurs naturally in photoassociation and photodissociation processes, and in solid-state structures [3][4][5][6][7]. Ultracold quantum gasses have been the subject of research to simulate solid-state systems [8][9][10][11][12][13]. In the context of the study of tunneling dynamics of ultracold atoms, a double-well potential with static barrier is a standard example.
Most of the tunneling phenomena were investigated with ultracold atoms in one dimension, either in a double-well potential or in periodic optical lattices, and focused on the atomic motion in the lowest band, i.e., particularly with the ground state [14][15][16][17][18][19]. Moreover, it was shown in a double-well potential that the inclusion of higher single-particle levels is fundamentally important for correlations [20,21]. Josephson-like dynamics of a Bose-Einstein condensate of rubidium atoms was investigated in the second Bloch band of an optical square lattice [22]. Unconventional orbital superfluidity in the P -and F -bands of a bipartite optical square lattice was achieved in [23,24].
Recently, we explored the tunneling dynamics of ultracold bosons in a symmetric double-well potential where the atomic motion lies in the ground as well as excited bands [25].
The interband quantum tunneling of ultracold atoms between the ground and excited bands can occur by breaking the symmetry of the two wells of a symmetric double-well potential, say, by generating an asymmetric double-well potential. This asymmetric double-well potential can be understood as a unit cell of a tilted optical lattice [26][27][28][29]. Also, an external electromagnetic field can induce interband transitions as demonstrated in the spectroscopy of Wannier-Stark levels [30]. Here, we are interested in the interband quantum tunneling in an asymmetric double-well potential, typically describes in the tilted optical lattice. Per definition, the interband tunneling in an asymmetric double-well potential occurs when the energy of the ground state on one side of the double-well coincides with the energy of an excited state on the other side. This leads to tunneling between those states which is resonantly enhanced by the energy matching, usually referred to as resonantly enhanced tunneling [31,32].
Resonantly enhanced tunneling phenomenon has been observed in various fields of research, such as in a nonlinear effect in a Mott insulator by creating a particle-hole excitation [33], in the presence of many-body coherences [34], in waveguide arrays [35], in an optical lattice with an external magnetic field producing a Zeeman splitting of the energy levels [36], in the twoterminal current-voltage characteristics of a finite superlattice [37], in solid states systems such as superlattice [5,38], and in accelerated optical lattice potentials [27,28].
Hitherto, resonantly enhanced tunneling is studied when the resonant condition of energy matching of the two wells is along the direction of the barrier [31,32,39]. We name it longitudinal resonant tunneling. In Ref. [32], the condensates were prepared in a one-dimensional optical lattice with an additional Stark force and determined the tunneling rate by the Landau-Zener formula [40][41][42][43][44]. Zenesini et al. experimentally investigated the impact of atom-atom interactions on the resonantly enhanced tunneling process and presented a complementary theoretical description of Landau-Zener tunneling for ultracold atoms in periodic potentials [31]. The out-of-equilibrium quantum mean-field and many-body dynamics of interacting bosons were recently explored in a onedimensional asymmetric double-well potential [39] using the multi-configurational time-dependent Hartree for bosons (MCTDHB) method [45][46][47]. Precisely, the interacting bosons were prepared in the ground state of the left well (harmonic potential) and the mechanism of the tunneling process in the resulting double-well was studied by analyzing the time-evolution of the survival probability, depletion and fragmentation, and the many-particle position and momentum variances [39].
Resonant tunneling in a two-dimensional (2D) double-well setup brings interesting new questions, especially on the role of transverse excitations. Also, in 2D there can be a resonant tunneling phenomenon which does not exist in one-dimension, namely, transversal resonant tunneling when the energy-matching condition of the two wells of an asymmetric double-well is satisfied along the transverse direction of the barrier. In this work, we focus on the longitudinal and transversal resonant tunneling for the initial bosonic structures of the ground and excited states in a 2D asymmetric double-well potential. Here the bosons are loaded in the left well of the double-well potential and the barrier is formed along the x-direction. In order to investigate the tunneling dynamics in the longitudinal resonant scenario, we consider two different structures of bosonic clouds, i.e., the ground and transversely excited (y-excited) states. Although, the ground state has a one-dimensional analog, to create a transversely excited state, one requires at least a 2D geometry [25]. Therefore, to investigate the effect of the longitudinal resonant tunneling on the very basic state which includes transverse excitations is of fundamental interest. Moreover, we investigate the impact of the transverse direction on the tunneling dynamics of the ground state at the longitudinal resonant tunneling condition. In the transversal resonant scenario, the bosonic structures are assumed to be the ground and longitudinally excited (x-excited) states. These two states have one-dimensional analogs but the transversal resonant tunneling does exist only in a 2D geometry. Therefore, the role of longitudinal excitations in transversal resonant tunneling can be examined. In order to accurately explore the out-of-equilibrium tunneling dynamics for both resonant tunneling scenarios, we solve the underlying time-dependent many-boson Schrödinger equation using the MCTDHB method [46,47]. We solve the quantum dynamics of all the bosonic clouds at the mean-field and many-body levels.
Tunneling dynamics for all the bosonic structures are analyzed by the time evolution of various physical quantities, namely, the survival probability, loss of coherence, depletion and fragmentation, and the many-particle position, momentum, and angular-momentum variances. Even when the bosons are fully condensed, the variance is a sensitive probe of correlations [48]. Therefore, to display the effects of the quantum correlations on the variances of different operators, we compare the mean-field and many-body variances in addition to the corresponding comparison of the survival probabilities. We notice that the rate of growth of the quantum correlations depends on the shape of the bosonic clouds, presence of transverse excitations in the system, and the resonant tunneling condition. The interconnection between the density oscillations and the variances are discussed both at the longitudinal and transversal resonant tunneling scenarios. We show that correlations have different consequences on the various quantities discussed in this work. In general, we ask how the transverse degrees-of-freedom, perpendicular to the junction, can influence the time evolution of various physical quantities in a 2D asymmetric double-well potential at the resonant tunneling scenario. We find that the time evolutions of the variances behave completely differently in the longitudinal and transversal resonant tunneling conditions.

II. THEORETICAL FRAMEWORK
The dynamics of N interacting bosons in a two-dimensional trap can be described by the timedependent many-body Shrödinger equation, [T (r j ) +V (r j )] + j<kŴ (r j − r k ), (2.1) whereT (r) andV (r) represent the kinetic and potential energy terms, respectively. The interaction between the bosons is repulsive and considered as a Gaussian function [49] W (r 1 − r 1 ) = λ 0 e −(r 1 −r 2 ) 2 /2σ 2 2πσ 2 with σ = 0.25 √ π. In Appendix B, we investigate the role of the width of the interparticle interaction potential and demonstrate explicitly that the mean-field and many-body physics found in this work do qualitatively not depend on this width. λ 0 is the interaction strength and Λ the interaction parameter. Λ = λ 0 (N − 1) where N is the number of bosons. Note that the model inter-bosons interaction does not have any qualitative impact on the physics described here.
The numerical implementation of the real-time dynamics employed in this work can be found in [70,71].
Here, we investigate the dynamics of a few bosonic clouds prepared as either the ground, transversely-excited or longitudinally-excited state for longitudinal and transversal resonant tunneling phenomena in two different double-well potentials. In case of longitudinal resonant tunneling, we consider that the bosons are initially prepared in the left well, of a double-well potential, where c is the asymmetry parameter. The bosons are taken in the state of either as the non-interacting ground, Ψ G = 1 √ π F (x, y), or as the transversely-excited, Fig. 1 shows the initial density distributions of Ψ G and Ψ Y . To explore the dynamics of the bosonic clouds, we suddenly quench the inter-particle interaction at t = 0 from Λ = 0 to Λ = 0.01π and simultaneously change the trapping potential from V L (x, y) to the longitudinally-asymmetric 2D double-well potential, Here the asymmetry implies the right well is lower than the left well. This is achieved by adding a linear term in the longitudinal direction. The mathematical form of V T (x, y) is (see Fig. 1(a)): In order to investigate the transversal resonant tunneling phenomenon, we assume the bosons are prepared either as the non-interacting Ψ G or longitudinally-excited state, in the left well, V L (x, y), of a transversely-asymmetric 2D double-well potential, V T (x, y). Here the asymmetry means the right well is wider than the left well. This is generated by making the transverse frequency of the trap spatially-dependent. V L (x, y) is functioned as Here, S(x) = 1 + (ω n − 1)sin 2 (x + 1)π 4 is a switching function for the transversal frequency, where ω n is the minimal value of the transversal frequency. The initial density distribution of Ψ X is presented in Fig. 1. For investigating the time evolution, as described in the longitudinal resonant tunneling, we quench the inter-particle interaction at t = 0 from Λ = 0 to Λ = 0.01π.
At the same time, we convert V L (x, y) to the transversely-asymmetric 2D double-well potential V T (x, y) (see Fig. 1 To accurately tackle the many-body physics involved in the two-dimensional bosonic Josephson and (e) show the initial density distributions for the ground, transversely-excited, and longitudinally-excited states, respectively. Panel (e) shows the x-and y-axes.

III. RESULTS AND DISCUSSIONS
Here, we divide the section into two parts, namely, longitudinal and transversal resonant tunneling. In each part, we investigate the time evolution of various physical quantities, such as the survival probability, loss of coherence, and the variances of the position, momentum, and angularmomentum many-particle operators when the symmetry of the double-well potential is broken and it reaches to the resonant tunneling condition. The quantities discussed here incorporate a detailed information of the time-dependent many-boson wave function, explicitly, the density, reduced one-particle density matrix, and reduced two-particle density matrix.

A. Longitudinal resonant tunneling
At first we investigate how the longitudinal resonant tunneling phenomena of Ψ G and Ψ Y are achieved by gradually increasing the asymmetry parameter, c, in the longitudinally-asymmetric 2D double-well potential described in Eq. 2.4. Further, we discuss the time-evolution of the different physical quantities, mentioned above, for Ψ G and Ψ Y at the resonant tunneling condition and compare them with the corresponding results of symmetric 2D double-well potential, which can be obtained at c = 0 from Eq. 2.4. The investigation is performed and compared at the mean-field and many-body levels of theory.
In this work, Ψ G and Ψ Y are prepared in the left well of a 2D double-well potential. When we gradually increase the value of c starting from zero, the symmetry of the double-well potential breaks and the left well becomes the upper well and, consequently, the right well becomes the lower well. For some special values of c, the one-body energy of the left well coincides with one of the higher one-body energy levels of the right well, resulting in an enhanced tunneling of bosons from the left well to the right well. This enhancement of tunneling is usually referred to as the resonant tunneling in a double-well potential, see in this context [31,32]. From Eq. 2.4, we can find out that at any asymmetry parameter c, the energy difference of one-body spectrum between the two wells becomes 4c. Analytically, for harmonic left and right wells, one can realize that the resonant tunneling occurs when 4c will be equal to an integer. Now, we show how the gradual increase of asymmetry between the two wells of a double-well potential influences the tunneling of particles. In Fig. 2, we present, at the mean-field level, the variation of the maximal number of particles tunneling from the left to the right well for the states Ψ G and Ψ Y with the asymmetry parameter c. This maximal number of particles in the right well is determined from +∞ x=0 +∞ y=−∞ dxdy ρ(x, y; t) N , where ρ(x, y; t) is the density of the bosonic cloud. Here, we start from the symmetric double-well potential. At c = 0, when the bosons are allowed to evolve in time, we observe that 100% of the bosons can tunnel from the left to right well, signifying the delocalization of the one-particle eigen-functions in a symmetric double-well. Now, a small increase of asymmetry between the wells makes the one-particle state to be only partially delocalized over both wells, leading to a suppression of tunneling of particles. Further increase of the asymmetry, beyond c = 0.1, the one-particle state slowly begins to be delocalized in one of the wells, and the tunneling of particles increases. At c = 0.25, we find well-delocalization of the one-particle state in one of the well which leads to a complete tunneling of particles or specifically, refer to as resonant tunneling condition, when one-body energy of the left well matches with the energy of one of the one-body excited states of the right well. Similarly, the second resonant tunneling appears at c = 0.5. As the system is weakly interacting, we find that the maximum tunneling of particles is 100% both for Ψ G and Ψ Y at c = 0, c = 0.25, and c = 0.5. Also, we observe a small difference of maximum tunneling of particles between the results of Ψ G and Ψ Y at off-resonant tunneling. Here, it is noted that if the bosons are prepared in the right well (lower well), the gradual increase of asymmetry between the wells reduces the maximal number of bosons tunneling to the left well. Further increase of asymmetry between the wells, at c 0.1, the bosons become trapped in the right well, and produces only a mild breathing motion, mainly due to the quenching the interaction.
Next, we focus on both the resonant tunneling values of c and study how the many-body correlations affects the dynamical behaviors of the survival probability, loss of coherence, and the many-particle position, momentum, and angular-momentum variances. Moreover, we will display the corresponding analysis of the results obtained for the symmetric double-well potential to serve as a reference. In addition, we will present a comparison study of the mean-field and many-body results both for Ψ G and Ψ Y .
In order to adequately capture the time evolutions of Ψ G and Ψ Y between the two wells of a double-well potential, we examine the survival probability in the left well, is the density of the bosonic cloud. We compare between the results of the mean-field and many-body dynamics of P L (t) in Fig. 3 for the asymmetry parameters, c = 0, 0.25, and 0.5. Here, we want to mention that, to have a proper comparison among all the results presented in this paper, the time-scale for the dynamics is set to be equal to the Rabi oscillations (t Rabi ) of the symmetric double-well trap (when c = 0). We find that t Rabi = 132.498 for the symmetric 2D double-well potential [25].
In Fig. 3, we observe, for both Ψ G and Ψ Y , back and forth of the density between the left and right wells with the essentially same frequency of oscillations at a particular value of c. In the noninteracting system, at t = 0, Ψ G and Ψ Y lie approximately on the lowest band along the direction of the barrier and when there is an interaction quench, more bands are coupled. Both the frequency and amplitude of the tunneling oscillations of Ψ G are essentially similar to those of Ψ Y at the meanfield level. But, at the many-body level, only the frequency of the tunneling oscillations remains very similar for Ψ G and Ψ Y at a particular value of c, while the amplitudes show a substantial difference as time progresses. The difference in the many-body tunneling amplitude between Ψ G and Ψ Y occurs as more states start to couple with time to Ψ Y in comparison to Ψ G . At c = 0.25 and 0.5, the ground state energy of the left well matches with energies of the second (which is two-fold degenerate, xΨ G and yΨ G ) and third (which is three-fold degenerate, xyΨ G , (x 2 − 1)Ψ G , and (y 2 − 1)Ψ G ) excited states of the right well, respectively. But, in the process of tunneling, Ψ G couples to xΨ G and (x 2 − 1)Ψ G of the right well at c = 0.25 and 0.5, respectively. Similarly, Ψ Y couples respectively to xΨ Y and (x 2 − 1)Ψ Y of the right well. As Ψ G and Ψ Y mainly couple to those excited states which have excitations only along the x-direction, we call this tunneling the longitudinal resonant tunneling.
For both Ψ G and Ψ Y , the short-time mean-field and many-body dynamics of P L (t) overlap for the interaction strength taken here. This condition imitates the so-called infinite-particle limit of the time-dependent many-boson Schrödinger equation [72,73]. Contrary to the mean-field dynamics, one can find a signature of the quantum correlations in the many-body dynamics of P L (t) in terms of incomplete tunneling of the densities and, consequently, the amplitude of the oscillations of P L (t) gradually decays. The decay in the amplitude of the many-body P L (t) signifies a collapse in the density oscillations, which is more pronounced for Ψ G than for Ψ Y . However, as c increases, the decay rates of Ψ G and Ψ Y become smaller. For a certain value of c, the decay rate of Ψ G is larger compared to the corresponding rate of Ψ Y , suggesting that in the tunneling process, the many-body effects develops in a different rate for different initial structures of the bosonic clouds. Explicitly, we observe that having transverse excitation delays the many-body process of the density oscillations collapse. Moreover, at c = 0.5, the time evolution of P L (t) exhibits a partial revival for both the states.
Signature of a growing degree of quantum correlations is already found in terms of decay in the amplitude of the time evolution of many-body P L (t). Now, we will discuss how this gradual increase of the many-body correlations can affect the coherence of the condensates, i.e., bosonic clouds of Ψ G and Ψ Y , when the resonant tunneling occurs. To compare the results, we also show the loss of coherence for the tunneling of the above considered bosonic clouds in the symmetric double-well potential as a reference. In Fig. 4, we present the time evolution of the condensate fraction, n 1 (t) N , obtained by diagonalizing the reduced one-particle density matrix of the timedependent many-boson wave-function (Eq. 2.2) [74,75]. The general feature of the occupation of the first natural orbital is found to be decreasing with time having a weak oscillatory background for all values of asymmetry parameters. It is observed that Ψ Y loses coherence faster than Ψ G for a particular value of c, suggesting that inclusion of the transverse excitation enhances fragmentation in Ψ Y [25]. We find that the rate of loss of coherence becomes slower when one move from c = 0 to the first resonant tunneling value at c = 0.25 and further slower at the second resonant tunneling.
In the context of occupations of the first natural orbital, it is worthwhile to mention the occupancy of higher natural orbitals. We find that, as the time passes by and fragmentation of the condensate grows, the occupancy of all higher natural orbitals gradually increases. The occupancy of higher natural orbitals have direct impact on the different many-body quantities discussed later, i.e, on the mechanism of many-body resonant tunneling in 2D double-wells. The detailed process of fragmentation with their convergences are discussed in the supplemental materials.
Having explicated the development of many-body correlations and their effects on the time evolution of the survival probability and coherence of the bosonic clouds, we find further information of the time-dependent many-particle wavefunction and especially, the signature of the depletion and fragmentation on the condensate. As the variance is a sensitive probe of correlation [48], we  mathematical form of the many-particle variance and its basic properties). As a general feature, we find that the many-body and mean-field values of 1 N ∆ 2X (t) are oscillatory in nature. Due to the correlations, the average value of the many-body 1 N ∆ 2X (t) deviates -it rather significantly increases -from the corresponding mean-field results for all asymmetry parameters and bosonic clouds. Although, we found in Fig. 4 that Ψ Y is always more fragmented than Ψ G , the deviations in the values of 1 N ∆ 2X (t) are always smaller for Ψ Y compared to Ψ G at a particular value of c, apart for short-times at c = 0.5. Therefore, one can say that the many-body dynamics is complicated and there is no one-to-one correlation between the different many-body properties in the junction, such as fragmentation and  Next, we consider the momentum variance per particle along the x-direction, 1 N ∆ 2P X (t), and present its time evolution at the mean-field and many-body levels for Ψ G and Ψ Y in Fig. 6. We explicitly show the possible quantum correlations effects on the many-body dynamics of at the resonant tunneling condition and compare with the corresponding results obtained for the symmetric double-well potential. Recall that 1 N ∆ 2P X (t) is relatively a more complex quantity than 1 N ∆ 2X (t) in the junction as the former one is controlled by and more sensitive to the shape of the orbitals. In Fig. 6, we observe, also for 1 N ∆ 2P X (t), two types of oscillations, smaller amplitude with higher frequency and larger amplitude with lower frequency. But there is a difference with respect In this context, it is worthwhile to briefly discuss the position and momentum variances along , respectively. We found in our work that the mean-field and many-body values of The results are compared at the mean-field and manybody levels for Ψ G and Ψ Y . At t = 0, in case of the symmetric double-well trap, can be calculated analytically and is found to be 2 and 7 for Ψ G and Ψ Y , respectively. At the first and second resonant tunneling, the initial values of 1 N For the symmetric double-well potential, we find that 1 N shows a decay in amplitude as fragmentation grows, which is more evident at c = 0.25. The rate Until now, we determined the effect of the transverse excitations on the dynamics of Ψ G and Ψ Y at the resonant tunneling condition within the mean-field and many-body levels of theory. The effects of the transverse degrees-of-freedom are analyzed in terms of various physical quantities.
The results prove that the transverse direction have substantial influences on the dynamics of the ground as well as excited states at the longitudinal resonant tunneling scenario in two spatial dimensions.

B. Transversal resonant tunneling
Proceeding to the transversal resonant tunneling, we want to emphasize that this scenario does only exist in 2D, thereby bringing to the front new degrees-of-freedom with the transverse direction, which we would like to explore. We have found that transverse excitations play a substantial role in longitudinal resonant tunneling. Here, we ask how and in what capacity longitudinal excitations play a role in transversal resonant tunneling.
At first, we examine how the transversal resonant tunneling can be achieved for two initial bosonic structures, Ψ G and Ψ X , in a transversely-asymmetric 2D double-well potential presented in Eq. 2.6. Here, we gradually decrease the transverse frequency in the right well, ω n , from an initial value (ω n = 1) which represents a symmetric double-well potential. Similarly to the longitudinal resonant tunneling, here also, the initial states are prepared in the left well of the double-well potential. Fig. 8 presents the maximal number of particles tunneling to the right well for Ψ G and Ψ X . The maximal number is extracted from the time-dependent solution of the Gross-Pitaeveskii is the density of the bosonic cloud. One can see from Fig. 8, at ω n = 1, that the maximal number of particles tunneling to the right well is 100%. Decreasing the value of ω n reduces the number of maximum particles which can tunnel to the right well for both initial bosonic clouds. The rate of decay of the maximal number of particles tunneling to the right well with ω n is smaller for Ψ X as it essentially lies in the first excited band along the direction of the barrier and thus feels a smaller barrier when it tunnels. Further decrease of ω n , the maximal number of particles tunneling to the right well grows. When the one-body energy of the left well coincides with one of the one-body higher energy levels of the right well, again the maximal number of particles tunneling to the right well reaches to almost 100%. We find only one resonant tunneling for the span of ω n considered here. For Ψ G and Ψ X , the transversal resonant tunneling occurs at ω n = 0.19 and ω n = 0.18, respectively. Also, it can be observed from the figure that for slightly off-resonant values of ω n , the maximal number of particles tunneling to the right well significantly decreases in comparison to the respective value of ω n at resonant tunneling for Ψ G , but it varies much weaker with ω n for ψ X .
The transversal resonant tunneling for Ψ G found in Fig. 8 is when Ψ G of the left well couples to (y 2 − 1)Ψ G of the right well. Similarly, Ψ X of the left well couples to (y 2 − 1)Ψ X of the right well to facilitate the transversal resonant tunneling of Ψ X , see Fig. 9. The resonant tunneling channels Ψ G → yΨ G and Ψ X → yΨ X are symmetry forbidden, at least at the mean-field level, as the trapping ω n of the right well for Ψ G and Ψ X (see Fig. 1(b)). The interaction parameter is Λ = 0.01π.
The same plot is highlighted in the inset by enlarging the resonance region. The data is generated from the time-dependent solution of the Gross-Pitaeveskii equation. We show here dimensionless quantities. 0.20, the tunneling frequencies for Ψ G are 1.41t Rabi , 2.30t Rabi , and 1.55t Rabi , respectively. While for Ψ X , the tunneling frequency is around 0.38t Rabi for all three ω n values.
At a fixed value of ω n , both the frequency and amplitude of the survival probability at the mean-field and many-body levels are practically the same for Ψ X . Whereas for Ψ G , although the frequency of the survival probability are very similar at the mean-field and many-body levels, the amplitudes of the survival probability show a significant deviation which is a maximal at ω n = 0.19. The gradual decay at the many-body level of dynamics of the survival probability signifies the growing degree of quantum correlations which will be quantitatively discussed later in terms of the loss of coherence in the system. All in all, Fig. 10 reflects that correlations develops faster for Ψ G compared to Ψ X . Now, we demonstrate the effect of growing degree of correlations on the loss of coherence of the bosonic clouds, Ψ G and Ψ X , at the transversal resonant tunneling and also at its vicinity.
In Fig. 11, we display the time-dependent occupation of the first natural orbital, namely the condensate fraction, n 1 (t) N . As observed in the longitudinal resonant tunneling, here also, n 1 (t) N For both the states, Ψ G and Ψ X , it is found that the frequency of oscillation of 1 N ∆ 2 Y (t) overlaps with the corresponding frequency of the survival probability. The amplitude of the many-body time evolution of 1 N ∆ 2 Y (t) decays with time, originating from the growing degree of fragmentation which is more evident for Ψ G . It is clear from Fig. 13 t) by the breathing oscillations. We notice that, for Ψ G , the transversal resonant tunneling gives rise to a beating pattern in the many-body dynamics . This beating pattern may be the consequence of the combination of different breathing frequencies. A dedicate study of many-body excitations in the transverselyasymmetric 2D double-well potential could resolve the many-boson states.
Proceeding, we examine the time evolution of the angular-momentum variance per particle, (t), for Ψ G and Ψ X at the transversal resonant tunneling scenario (see Fig. 14). The mean- (t) which is significantly more prominent for Ψ G and hardly visible for Ψ X .

IV. SUMMARY AND CONCLUSIONS
The paradigm of a bosonic Josephson junction, in which bosons can tunnel back and forth between two potential wells, is shown to be very rich at resonant tunneling condition in two spatial Moreover, at c = 0.5 we find a partial revival pattern in the mean-field as well as the many-body survival probabilities. Although the rate of collapse of density oscillation for Ψ Y is slower compared to Ψ G at a fixed value of c, the presence of transverse excitation in Ψ Y makes it more fragmented than Ψ G at a given time. Compared to the symmetric double-well potential, we observe that both initial states become less fragmented at c = 0.25 and further less fragmented at c = 0.5.
Referring to the transversal resonant tunneling scenario, we scan the frequency along the transverse direction in the right well, ω n , starting from the value ω n = 1 (symmetric 2D double-well) and reaching down to ω n = 0.17, and find that the first resonant tunneling in the transverse direction occurs at ω n = 0.19 and 0.18 for Ψ G and Ψ X , respectively. This resonant tunneling is achieved when the bosons tunnel from Ψ G of the left well to (y 2 − 1)Ψ G of the right well. Similarly, for Ψ X , the bosons tunnel from Ψ X to (y 2 − 1)Ψ X . The time evolution of the many-body survival probability of Ψ G and Ψ X signifies that the rate of collapse of density oscillations is faster for Ψ G compared to Ψ X . This collapse of density oscillations occurs due to the growing degree of many-body correlations which is graphically shown in terms of the loss of coherence in the system.
To obtain further information of the time-dependent many-particle wavefunction, we present the many-particle variances of the observables, The present work can inspire several promising and interesting future research directions. An immediate extension would be the dynamics of even more intricate bosonic structures, e.g., vortex state or a mixture of vortex states induced by Raman process using light with orbital angular momentum [77][78][79][80][81] or by synthetic magnetic fields [82] in the longitudinal and transversal resonant tunneling conditions. Also, to apply a linear response theory [76] to accurately calculate the different breathing frequencies involved in the time evolution of different quantities at the resonant tunneling scenario would be challenging. Finally, in the long run we could envision that additional geometries would open up in three spatial dimensions.

APPENDIX A: MANY-PARTICLE VARIANCE
The variance of an observable,Â, is determined by the combination of the expectation values ofÂ andÂ 2 . Here the expectation value ofÂ = N j=1â (r j ) depends on the one-body operators while the expectation ofÂ 2 is a mixture of one-and two-body operators,Â 2 = N j=1â 2 (r j ) + j<k 2â(r j )â(r k ). The variance can be written as [66] 1 where {φ j (r; t)} are the natural orbitals, {n j (t)} the natural occupations, and ρ jpkq (t) are the elements of the reduced two-particle density matrix, ρ(r 1 , r 2 , r 1 , For one-body operators which are local in position space, the variance described in Eq A.1 boils down to [47] 1 N ∆ 2Â (t) = dr ρ(r; t) Nâ 2 (r)−N ρ(r; t) Nâ (r) 2 + dr 1 dr 2 ρ (2) (r 1 , r 2 , r 1 , r 2 ; t) N a(r 1 )a(r 2 ), (A.2) where ρ(r; t) is the time-dependent density. Eq A.1 describes the variances when the center-of-mass of the bosonic clouds are at the origin and the wavefunction is Ψ(0, 0) (or denoted in short by Ψ).
As the initial states considered in this work are prepared at the position (a, b) = (−2, 0), here we present a general form of variances by incorporating the translated wavefunction, i.e., Ψ(a, b) (or denoted in short by Ψ ab ). For the position and momentum operators, the variances do not change due to the translated wavefunction [66], therefore, . But, for the angular-momentum variance, the situation becomes intricate and expressed as [66] 1 N ∆ 2  In Fig. B2, we present the many-body dynamics of at c = 0.25. From Fig. B2 (a), one can observe that the rate of loss of coherence is slightly faster for σ = 0.25/ √ π compared to other values of σ. It is found that decreasing the width of σ of the interparticle interaction potential gradually diminishes the quantitative difference between the dynamical behaviors of a particular quantity. Interestingly, the variances which are known to be sensitive many-body quantities due to depletion are hardly influenced by the width of σ. Although there is a small quantitative difference in the long-time fragmentation dynamics for the choices of σ, the many-body physics described for the longitudinal resonant scenario is robust which can be seen in the dynamical behavior of the different quantum mechanical properties, see Fig. B2 The dynamics are computed with N = 10 bosons and the interaction parameter Λ = 0.01π. The many-body dynamics are computed with M = 6 time-dependent orbitals.
The quantitative many-body results are found to be weakly dependent on σ. We show here dimensionless quantities. Color codes are explained in each panel.  use the numerical implementation in [4,5]. We demonstrate the details of the mechanism of fragmentation in terms of the occupancy of the higher natural orbitals for both resonant tunneling conditions which were not discussed in the main text. Moreover, we graphically display the position and momentum variances [6,7] along the y-direction for the longitudinal resonant tunneling at the mean-field and many-body levels (not shown in the main text). We have found in the main text that for Ψ G and Ψ X , the transversal resonant conditions are satisfied at ω n = 0.19 and ω n = 0.18, respectively, and demonstrated the variances only at the resonant values of ω n . Here, we present all the many-body variances [8,9], discussed in the main text, along with their convergences  Having analyzed the convergence of n 1 (t) N , we investigate the further details of how the coherence is lost in the system by examining how fragmentation is developed. Fig. S2 displays the time evolution of the occupancies of the higher natural orbitals per particle, n j>1 (t) N , for Ψ G and Ψ Y with their convergences with the number of time-dependent orbitals. Hence, convergence of the As presented in the discussion of loss of coherence, here also, we plot n j>1 (t) N for the asymmetry parameters, c = 0, 0.25, and 0.5. The higher natural orbitals which have significant amount of i.e., c = 0.5. In the main text, we discussed that, at c = 0.5, 1 N ∆ 2X (t) and 1 N ∆ 2P X (t) have two kinds of oscillations, namely, small frequency with large amplitude (due to the density oscillation) and high frequency with small amplitude (due to the breathing mode oscillation). The panels in Fig. S3 show that both variances are well converged with the orbital numbers. Even the oscillation which has high frequency with small amplitude, found with lower orbital numbers, overlaps with the corresponding results calculated from the higher orbital numbers.
As discussed in the main text, the many-particle position and momentum variances per particle along the y-direction,     In order to learn how the fragmentation develops in the transversal resonant tunneling scenario and at its vicinity, we explore the time-dependent occupancies of higher natural orbitals per particle, n j>1 (t) N , along with their convergences with the number of orbitals for Ψ G and Ψ X . We graphically present only n 2 (t) N and n 3 (t) N in Fig. S8 as they have occupations in appreciable amount. Among the three different values of ω n , n 2 (t) N is maximally occupied at ω n = 0.19 for Ψ G and at ω n = 0.18 for Ψ X , as these values of ω n produce the transversal resonant tunneling for the respective states.
We observe the oscillatory nature in the the dynamics of n 2 (t) N and n 3 (t) N for Ψ G , whereas it is hardly visible for Ψ X .
Figs. S9 and S10 show the many-body position and momentum variances per particle along the x-direction, Unlike Ψ X , we find significantly different many-body dynamics of In Figs. S11 and S12, we display the many-body position and momentum variances per particle along the y-direction, the amplitude of the many-body time evolution of 1 N ∆ 2 Y (t) decays with the growing degree of correlations, see Fig. S11, and the decay rate is maximal at ω n = 0.19 and minimal at ω n = 0.20.
At the resonant value for Ψ G , ω n = 0.19, the amplitude of 1 N ∆ 2 Y (t) reaches upto around 15, while at the other frequencies they are comparatively small, i.e., around 4.5 at ω n = 0.18 and around 6.5 at ω n = 0.20. For the slower rate of growth of the fragmentation for Ψ X in comparison to Ψ G , the decay rate of the many-body 1 N ∆ 2 Y (t) is hardly visible even at the resonant value of ω n , see Fig. S11. By comparing the many-body 1 N ∆ 2 Y (t) at different frequencies for Ψ X , we observe that at the resonant condition for Ψ X , i.e., ω n = 0.18, the amplitude of oscillation reaches upto 12 whereas at ω n = 0.19 and ω n = 0.20, the respective amplitudes reach upto around 10.5 and 9, respectively.
Focusing on the momentum variance along the y-direction for Ψ G , we notice that at ω n = 0. The main text and the supplemental material so far describe all quantities, i.e., the survival probability, loss of coherence, fragmentation, and the variances for the transversal resonant scenario with 128 × 128 grid points. In order to verify the convergence with the grid points, we repeat our computation with 256×256 grid points for both objects, Ψ G and Ψ X and of course the computations become demanding due to the increased grid density. To demonstrate the convergence with the grid points, we display 1 N ∆ 2 L Z (t) which is the most sensitive quantity discussed in this work. Fig. S14 displays the convergence of the many-body 1 N