{'fraction_non_alphanumeric': 0.03645886474169303, 'fraction_numerical': 0.014828418868822909, 'mean_word_length': 4.361779722830051, 'pattern_counts': {'":': 0, '<': 3, '<?xml version=': 0, '>': 2, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 5, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Dynamics in asymmetric double-well condensates

The dynamics of Bose-Einstein condensates in asymmetric double-wells is studied. We construct a two-mode model and analyze the properties of the corresponding phase-space diagram, showing in particular that the minimum of the phase-space portrait becomes shifted from the origin as a consequence of the nonvanishing overlap between the ground and excited states from which the localized states are derived. We further incorporate effective interaction corrections in the set of two-mode model parameters. Such a formalism is applied to a recent experimentally explored system, which is confined by a toroidal trap with radial barriers forming an arbitrary angle between them. We confront the model results with Gross-Pitaevskii simulations for various angle values finding a very good agreement. We also analyze the accuracy of a previously employed simple model for moving barriers, exploring a possible improvement that could cover a wider range of trap asymmetries.


I. INTRODUCTION
The two-mode (TM) model has been studied and extensively applied to symmetric double-well atomic Bose-Einstein condensates in the last years [1][2][3][4][5][6][7][8]. The dynamics of such a model relies on assuming that the condensate order parameter can be described as a superposition of wave functions localized on each well with time-dependent coefficients. The predicted Josephson and self-trapping regimes [1,2] have been experimentally observed by Albiez et. al. [5]. Also ring-shaped condensates with two radial barriers forming a double well have been theoretically studied in the case of dipolar [7] and contact [9] interactions.
On the other hand, the situation is quite different in the case of asymmetric double wells, since no similar theoretical development for such configurations has been reported so far. This kind of systems has become increasingly interesting due to recent experiments [10][11][12], where toroidal traps with two radial barriers moving symmetrically were investigated, which involves a much richer type of dynamics. In fact, in addition to representing the first experimental realization of a SQUID analog with a Bose-Einstein condensate [10,11], such experiments show a critical barrier velocity above which atoms become compressed on one side and expanded on the other side, in close analogy to the transition from dc to ac Josephson effects. Such an ac regime presents a higher degree of complexity since excitations like solitons and vortices may be shed into the rarefied portion of the condensate, giving rise at their eventual decay to an additional resistive current [12]. If we restrict ourselves to small displacements of the barriers from the symmetric configuration, it has been shown in Ref. [10] that a straightforward generalization of the symmetric TM model to such a dynamic configuration works well. However, it is easy to understand that any first step to achieve an analogous model valid for more general barrier movements should necessarily involve the study of an asymmetric TM model, which constitutes the main goal of the present work. Thus, focusing on an arbitrary configuration of fixed barriers, we begin in Sec. II by analyzing the properties of the stationary states of an asymmetric pair of weakly coupled condensates, from which both localized states are derived. In fact, the nonvanishing overlap between the Gross-Pitaevskii (GP) ground state and the excited stationary state turns out to determine the position of the minimum in the phasespace diagram. We derive the full set of TM model parameters and further introduce in Sec. III the corrections in the interaction energy parameter using the proposal of Ref. [9] adapted for an asymmetric system. In Sec. IV A, we describe the system we use in our model applications and simulations following the experimental settings of Ref. [10]. The corresponding phase-space diagram is obtained in Sec. IV B, where we compare the TM results with GP simulations, showing that our model, corrected by the modified effective interaction parameters, yields a much more accurate dynamics. The case of moving barriers is considered in Sec. IV C, where we derive equations of motion similar to those employed in Ref. [10], discuss the importance of terms disregarded in this approach and explore a possible improvement. Finally, in Sec. V we summarize our work giving some concluding remarks.

II. ASYMMETRIC TWO-MODE MODEL
The TM model dynamics has been extensively studied in symmetric double-well potentials [2,3]. The commonly used ansatz for the wavefunction reads where ψ 1 (r) and ψ 2 (r) are real, normalized to unity, localized wave functions at each well. The complex timedependent coefficients are written as b k (t) = N k /N e iφ k (k = 1, 2), where φ k and N k represent the phase and particle number in the k well, respectively, and N denotes the total number of particles. We will use the same ansatz for our generic asymmetric configuration, with localized states constructed from the ground state ψ G (r), which we assume real and positive, and from the excited stationary state ψ E (r), which we assume real and negative (positive) in the '1' ('2') well. Both states are supposed to be normalized to one. Then, the localized states read, where β ≡ ψ G |ψ E = d 3 r ψ G (r)ψ E (r). Here it is important to remark that although the stationary states of the asymmetric case present a nonvanishing overlap, the corresponding localized states are indeed orthogonal by construction. Writing the stationary states in terms of the localized ones we have, Thus, we may see from the above equations that both stationary states have identical populations at each site. As usual, in order to obtain the TM dynamics, we introduce the order parameter into the time-dependent GP equation, Projecting it onto ψ 1 (r) and ψ 2 (r), and integrating both equations using the hopping and on-site energy parameters given in Eqs. (A.1) to (A.8) of the Appendix, one obtains ih These equations include all the terms introduced for a symmetric system in the improved two-mode model [3]. In terms of imbalance Z = (N 2 − N 1 )/N and phase difference φ = φ 1 − φ 2 , we obtain the following equations of motion where with the corresponding parameter definitions given in the Appendix. The above equations of motion can be obtained from the following Hamiltonian, using the fact Z and φ are canonical conjugated coordinates i.e.,Ż = −∂H/∂φ andφ = ∂H/∂Z. The phase-space portrait (Z, φ) of such a Hamiltonian exhibits a minimum (Z 0 , 0), and for strongly interacting systems, a saddle (Z 0 , π) and two maxima (±Z M , π). As can be easily deduced from Eqs. (4)-(5), the overlap between the ground and excited states determines the stationary imbalance Z 0 = β. Thus, we may see that minimum and saddle will be shifted in a Z = β value from the corresponding locations on the φ axis in the case of the symmetric double well. Using the TM model parameters, Z 0 may be approximated disregarding almost negligible terms by The separatrix between closed, with a bounded phase (BP), and open, with a running phase (RP) orbits arises from the condition that the energy corresponds to the saddle point, H(Z, φ) = H(Z 0 , π). Particularly, for φ = 0 the separatrix yields a critical imbalance Z c given approximately by where the plus (minus) sign corresponds to the separatrix above (below) the BP orbits.

III. EFFECTIVE INTERACTION EFFECTS
The inclusion of effective interaction effects in the TM model of a symmetric double-well system has shown to provide an accurate correction to the disagreements with GP simulations [9,13]. Such a correction takes into account the density deformation of the GP mean-field term when varying the imbalance, which results in a net reduction of the interaction energy parameter. Within the Thomas-Fermi approximation, the corresponding reducing factor can be analytically obtained [13], whereas a numerical calculation using the ground state density has shown to provide accurate values in more general cases [9]. We note that a similar correction (about a 20 % reduction) has been introduced in the plasma oscillation frequency of a tunable superfluid junction, which was shown to be crucial to accurately describing experimental results [14].
Here we follow the same procedure of Ref. [13] and revise the term related to the interaction energy: where In asymmetric double-well systems, in principle, the deformation of each localized density is not the same and thus the effective interactions U ∆N1 1 and U ∆N2 2 [9,13] should be treated separately. The idea behind the method is that nonequilibrium states can be well aproximated by localized on-site states corresponding to the instantaneous population at each well, N 0 k + ∆N k , where N 0 k denotes the population of the k-well for N particles. Such localized states can be obtained from the stationary states of systems with a total number of particles different from N , whose localized on-site densities (normalized to unity) are denoted by ρ ∆N k k (r) in (19). More details about this calculation will be given in Sec. IV A.
We will assume that analogously to the symmetric case [9,13], the following first-order approximation remains valid in any case where the parameter α k may be numerically evaluated according to the procedure of Ref. [9]. Using 2 ∆N1 Then, with the following definitions of effective on-site energy dependent parameters, and introducing the correction (21) into the equation of motion (10), we obtain which is consistent with the following 'effective' Hamiltoniañ The model represented by the equations of motion (9)-(25) and the Hamiltonian (26), will be called as the effective two-mode (ETM) model in what follows. With respect to the new phase-space portrait derived from this Hamiltonian, it can be easily verified that the position of the minimum remains located at Z 0 , which may be also approximated by whereas the shape of the orbits may differ from that obtained with the bare parameters, as will be shown in the following sections. We also note thatŨ a andÃ become reduced with respect to U a and A, as seen from Eqs. (22) and (23), provided the parameters α k are positive. On the other hand, we remark that the parameterB arises from the combined effects of interaction and asymmetry.

A. The system
We describe in what follows the system utilized in our simulations and model applications. All the trapping parameters and condensate details have been chosen to reproduce the experimental setting of Ref. [10]. The trapping potential can be written as the sum of a part that depends only on x and y and a part that is harmonic in the tightly bound direction z: The above potential consists of a superposition of a toroidal term V T (r) (r 2 = x 2 + y 2 ) and the radial barrier term V B (x, y). The toroidal potential was modeled through the following Laguerre-Gauss optical potential [15] where V 0 corresponds to the depth of the potential and r 0 the radial position of its minimum. We have used scaled units referenced to a chosen unit of length denoted by L 0 (in our case L 0 =1 µm). Energy and time units were defined in terms of L 0 : where m denotes the mass of a condensate atom. For the present case of 87 Rb atoms we have E 0 /k B = 5.5298 nK and T 0 = 1.3813 ms. The barrier was modeled as where Θ denotes the Heaviside function with θ 1 = θ and θ 2 = π − θ. The parameter θ will be assumed as time dependent in the case of moving barriers. We have utilized, according to Ref. [10], the following trap parameters: .07 nK, and λ b = 1 µm. These barrier parameters yield a full-width at half-maximum of the barrier nearly below 2 µm, which is in agreement with the experimental data leading to a tunnel junction. We have assumed a high λ = 8 value yielding a quasi-bidimensional condensate and allowing a simplified numerical treatment [16]. So, stationary states are written as the product of a two-dimensional (2D) wave function ϕ(x, y) and a Gaussian wave function along the z coordinate, λ 1/2 π 1/2 e − λz 2 2 . Thus, assuming barriers remaining at rest, the GP equation for the former reads [16] − 1 2 where a = 98.98 a 0 denotes the s-wave scattering length of 87 Rb, a 0 being the Bohr radius. In Fig. 1 we depict the 2D particle density |ϕ(x, y)| 2 for the ground state at different positions of the barriers. According to the notation of previous sections we will call the top and bottom wells of Fig. 1 as '1' and '2', respectively. In the following we will restrict our calculations to a system with N = 3000. As regards the excited state, we have obtained its wavefunction by evolving in imaginary time an initial wavefunction identical to that of the ground state in site '2', whereas we introduced a change of sign in the site '1', i.e., a wavefunction positive in site '2' and negative in site '1', a feature that turns out to persist until final convergence to the stationary excited state. This procedure works well for all configurations below θ = 0.394π, while for more asymmetric systems, the energy gap between the excited and the ground states becomes so low that the imaginary-time evolution leads to a 'decay' to the ground state. Here it is interesting to notice that, generalizing the TM result of the symmetric case [14], such a gap reads (J a + Z 0 J b ) 1 − Z 2 0 > 0, which vanishes for Z 0 → 1, as expected. To illustrate the method we used to calculate the coefficients α k in Eq. (20), we depict in Fig. 2 where N 0 k is obtained from the projection of the ground state onto the k-localized state for N = 3000, whereas N k is calculated analogously, but with a different total number of particles. Thus, according to Eq. (20), we have extracted the values of α k from the slope of the lines. We note that α 1 turns out to be larger than α 2 because the smaller condensate should present the larger deformation for an identical change in the particle number.

B. Phase-space portrait and dynamics for static barriers
We first note that for 0 ≤ θ ≤ π/2, the overlap between the ground and excited states verifies 0 ≤ β ≤ 1, and the same occurs for the position of minimum and saddle since Z 0 = β. In Fig. 3 we depict the phase-space diagram of Hamiltonian (15) for θ = π/5. In this case we have Z 0 = 0.5163, which corresponds to the Z coordinate of minimum and saddle in such a figure. The separatrix between BP and RP orbits has been numerically obtained and has been denoted by dashed lines in Fig. 3. We remark that the values Z c = Z 0 ± 0.0358 derived from Eq. (17) are in well accordance with the intersections of the dashed lines and the vertical axis.
In Fig. 4 we depict the phase-space portrait arising from GP simulations and from ETM and TM models. We notice that the minimum of Hamiltonian (26) remains located at Z 0 (cf. Eq. (27)), whereas the shape of the orbits differ from that observed in Fig. 3. Particularly, the separatrix between closed and open orbits covers a wider range of Z values, as shown in Fig. 4 through the locations of the critical imbalance Z c arising from the TM model (cf. Eq. (17)) indicated by red dots, and those arising from the ETM model (blue stars). Here it is worth noticing also that the value of Z c obtained from GP simulations coincides with the corresponding ETM result. Therefore, we may conclude that the matching between ETM model and GP simulation results turns out to be much better than that of the plain TM model. We have depicted in Fig. 5 the time evolution of imbalance and phase difference for θ = π/5 and θ = π/3 arising from GP simulations, together with the corresponding TM and ETM results. It is remarkable that the excellent agreement between GP and ETM time evolutions persists even for most asymmetric configurations.

C. Moving barriers
The experimental results of Ref. [10] were well reproduced from a simple model [17], that adapted the TM equations of motion of a symmetrical configuration to the case of moving barriers by simply taking into account the effect of a θ-dependent equilibrium imbalance Z 0 (θ(t)). Such a simple model for moving barriers (SMMB), can be easily derived from our equations of motion (9)-(25) by approximating all the model parameters by their values of the symmetric configuration, neglecting terms proportional to the small parameter I, and replacingÃ → Z 0 (θ(t))Ũ a N according to (27). Thus we obtain,hŻ where J = J a andŨ =Ũ a respectively denote hopping and on-site energy parameters given by the corresponding values of the symmetric case (θ = 0). The barrier movement in these equations is represented by Z 0 (θ(t)), and of course only small departures from the symmetric configuration should be expected to be well reproduced. Particularly, the experiments in [10] were restricted to θ < π/8, and in this paper we will explore the dynamics for a wider range of asymmetric final configurations. By numerically analyzing Z 0 as a function of θ, we have observed a linear behavior Z 0 = α θ (2πα = 5.1585), except for values reaching θ ≃ 0.4π where both barriers begin to overlap (Z 0 → 1). We note that the approximation 2πα = 4 used in Ref. [10], which amounts to assuming a linear behavior of Z 0 (θ) up to θ = π/2 (Z 0 (θ = π/2) = 1), corresponds to the limit of a negligible barrier width and also neglecting any healing length arising from the presence of barriers. For simplicity, we will assume in this paper barriers moving with a constant angular frequency f b , thus θ(t) = 2πf b t and we may approximate In Figs. 6 and 7 we depict the time evolution of imbalance Z and phase difference φ for barriers moving with f b = 0.1 Hz and 0.5 Hz, respectively. The solid lines correspond to GP simulation results, while the dashed ones correspond to the SMMB results. There we may observe typical behaviors of the so-called dc-and ac-Josephson regimes for barrier frequencies 0.1 and 0.5 Hz, respectively [10,17]. In fact, the dc-Josephson regime, which occurs below certain critical barrier frequency, is characterized by small oscillations around Z 0 (t) and a bounded phase difference, and for this reason we shall call it as the dynamical bounded-phase (DBP) regime. On the other hand, for a large enough barrier velocity a compression dynamics occurs [10,12], which defines the ac-Josephson regime. Such a regime exhibits an unbounded phase and so we shall call it in what follows as the dynamical running-phase (DRP) regime.
It is instructive to analyze the short-time behavior of both regimes using the SMMB equations (35)-(36). For small barrier angular frequencies within the DBP regime, an analytically tractable dynamics arises by linearizing the SMMB for Z ≪ 1 and φ ≪ 1. Thus we obtain the approximate solutions, with ω p = Ũ N J h 2 , that qualitatively resembles the GP dynamics of Fig. 6, mainly the fact that Z oscillates around Z 0 (t) and the bounded phase remains confined to negative values. For longer times, increasing differences between the GP simulation results and those of the SMMB are observed for θ > 0.15 π in the phase evolution at the right panel of Fig. 6.
On the other hand, in the DRP regime, we may approximate for large barrier angular frequencies which yields an unbounded monotonically decreasing phase. This qualitatively reproduces the behavior of the phase difference shown in the right panel of Fig. 7, particularly the increasing negative slope that is observed along the evolution. As regards the imbalance shown in the left panel, we may see that the GP behavior is only reproduced by the SMMB at very short times (t < 60 ms). For longer times the SMMB completely fails to describe the GP dynamics. In particular, from the inset in this figure, it may be seen that the model asymptotically oscillates around Z ≃ 0.02, whereas the GP simulation shows an almost quadratic behavior, approaching Z = 1. Note that any realistic dynamics should tend to Z = 1 where the barriers are superposed. The above failure of the SMMB in the DRP regime is easily understood from a simple inspection of the density and phase snapshots shown in Fig. 8. In fact, it is clear that any variant of the TM model, like that yielding the SMMB, is supposed to rely on assuming that the shape of the density should at any time resemble that of a stationary state (Fig. 1), and the phase should remain almost homogeneous at each well. This is indeed the case for the snapshots of f b = 0.1 Hz in Fig. 8, which are consistent with the quite acceptable results for the SMMB in Fig. 6, and also suggest that an 'improved SMMB' covering a wider range of asymmetric configurations could be eventually devised for the DBP regime. On the other hand, the snapshots of f b = 0.5 Hz show a very different situation, with important deformations in the density, as compared to that of the ground state, and with clear inhomogeneities in the phase at the bottom well, which reflects the formation of excitations like vortices in the rarefied portion of the condensate, as also discussed in Ref. [12]. None of these features could be taken into account in any simple model like the SMMB. In Fig. 9 we depict Z(t) − Z 0 (t) as a function of φ(t) obtained from GP simulations and from the SMMB for several frequencies. It may be seen that the DBP oscillations are confined within an oval-like region, whereas the DRP evolutions are localized on the negative Z − Z 0 half-plane and can acquire any phase value. As regards the DBP oscillations, it is important to remark that they are not closed orbits, since each new loop does not exactly reproduce the previous one. Moreover, after a number of DBP loops, a given orbit (e.g., f b = 0.25 Hz) may 'decay' to the DRP regime. So, for the sake of clarity, we have only plotted in Fig. 9 that part of the trajectory corresponding to the first period of each DBP evolution. In addition, we point out that the asymmetry parameter θ in Fig. 9 does not exceed in any case the experimental limit of π/8.
It is interesting to analyze the turning points of the above trajectories, where the time derivatives of Z − Z 0 or φ vanish. The φ turning points are located, as seen in Fig. 9, approximately at the origin and on the negative φ axis for the DBP loops, while they are absent in the DRP regime, as expected. On the other hand, the Z − Z 0 turning points change from being a couple of almost vertically aligned points for the DBP loops (e.g. the black circles for f b = 0.25 Hz in Fig. 9), to becoming pairs of a maximum at left and a minimum at right, located almost symmetrically with respect to the vertical line φ = −π/2, as indicated by the red stars for f b = 0.2775 Hz. However, there are no such turning points above certain frequency, as observed for f b = 0.5 Hz in Fig. 9. This may be deduced from the following estimate of the phase φ of the Z − Z 0 turning points which stems from Eqs. (35) and (37) using the approximation Z 2 ≪ 1. First it is convenient to compare the above prediction with the phase of the first SMMB turning points for several barrier frequencies in the inset of Fig. 9. We may see that there is an excellent agreement between the above formula and such points, while we have found that this approximation remains quite acceptable for subsequent periods within θ < π/8, as observed for the red stars in the main plot of Fig. 9. Now, one may define two critical frequencies f ≃ J/(h2πα) ≃ 0.38 Hz, there are no more turning points, which corresponds to the barrier frequency leading to the critical current of Ref. [10]. On the other hand, the minimum barrier frequency f above which there are no more φ turning points, with the corresponding lack of any DBP regime (0.2775 Hz for GP results), corresponds in the inset of Fig. 9 to the transition from black circles to red stars.
The increasing discrepancies between the GP simulation results and those of the SMMB above certain asymmetry shown in Fig. 6 arise, as already pointed out, from the use of the parameters of the symmetric system in the SMMB. A first attempt to quantify the extent to which such an approximation could affect the SMMB results, can be evaluated by appreciating the differences in the value of the ETM model parameters in Table I with respect to those of the symmetric case. Note that the only nonvanishing parameters of the symmetric configuration are those of the SMMB, U a =Ũ and J a = J, and the neglected parameter I. On the other hand, although all the remaining parameters become finite for asymmetric configurations, for θ ≤ π/8 only a little effect of such asymmetries on the SMMB accuracy should be expected to occur since,Ũ a and J a vary less than 10 %, we have Z 0 J b ≪ J a , andÃ ≃Ũ a N Z 0 (cf. Eq. (27)) is well fulfilled in Table I. This explains the good agreement between GP and SMMB results observed in Fig. 9 and also the agreement with experimental data reported in Ref. [10]. Now, the simplest improvement to the SMMB in order to take into account the evolving trap configuration may be provided by an immediate generalization of the ETM model with its parameters depending on the instantaneous  trap asymmetry, i.e., for a time-dependent θ. To test such a possibility in the DBP regime, we depict in Figs. 10 and 11 the derivatives of imbalance and phase difference as functions of time arising from GP simulations for f b = 0.1 Hz, and compare such results with those given by the SMMB, and with the corresponding values arising from the ETM equations (9)-(25) for a variable (time-dependent) θ with Z and φ taken from the GP simulation results. Then we may observe in Fig. 10 that the SMMB prediction for the imbalance derivative (blue solid line) shows an increasing departure from the GP results for θ above 0.15 π. On the other hand, the results arising from the ETM model (yellow stars) show a better agreement with the GP simulation values for the whole time evolution.
As regards the calculations of the phase derivative, first we want to remark that we have found that the hopping terms (terms in J a , J b , I, and J) in the equations (10) for TM, (25) for ETM, and (36) of the SMMB turn out to be all negligible with respect to the on-site energy dependent terms, so we have disregarded their contribution in the calculations corresponding to Fig. 11. Thus, again one finds that the SMMB results increasingly differ from the GP results for θ above 0.15 π, and also that the ETM calculation clearly improves the agreement with the simulation results. In addition, it is possible in this case to test the accuracy of the ETM model results versus those of the plain TM model, as observed in Fig. 11.

V. CONCLUSION
We have developed a two-mode model of a Bose-Einstein condensate in an asymmetric double-well. Taking into account effective interaction effects as proposed previously for symmetric configurations, we have introduced corrections to such a model in order to improve the agreement with simulation results. Thus, we applied this formalism to a recently explored experimental setting of a toroidal trap split into an asymmetric double-well condensate by means of a pair of radial barriers. We have found that the qualitative agreement with the GP simulation results arising from the plain asymmetric TM model, became enhanced to the extent of an excellent concordance when the effective interaction effects were considered, for practically the whole range of asymmetric configurations.
We have explored the range of validity of the simplest theoretical model for moving barriers, previously utilized to describe experimental results, finding that it should fail for larger departures from the symmetric configuration than those considered in the experiment, and also for the running-phase regime of higher barrier frequencies. Finally, we have performed an analysis of a possible improvement of this model, which consists in an immediate generalization of the asymmetric ETM model to a moving-barrier configuration, simply by employing parameters depending on the instantaneous trap asymmetry. Our results obtained for the DBP regime pave the way for further studying more complex dynamics driven by different kinds of barrier movements.

Appendix: TM model parameters
The TM model parameters read, We note that according to the definitions (11) and (14), A and U a turn out to be on-site energy dependent parameters, as they are built from the on-site energy parameters ε k and U k . On the other hand, K is the standard hopping coefficient, while the remaining hopping parameters F 12 , F 21 and I have been first introduced by Ananikian and Bergeman in Ref. [3]. Particularly, N F 12 and N F 21 can be interpreted as additional contributions to the tunnelling obtained from a modified Hamiltonian that includes the interaction term as an effective potential [4]. We remark that in both references the authors assume symmetric traps with F 12 = F 21 , a parameter which in our case would only contribute to the hopping parameter J a (12). On the other hand, in an asymmetric trap, the difference between F 21 and F 12 gives rise to the additional hopping parameter J b (13). As regards the parameter I, it corresponds to atom-pair tunnelling processes and, although it turns out to be in most cases negligible, it was recently shown that for enough strong interaction regimes a quantum phase transition driven by such atom-pair tunnelling events should be expected to take place [18].