1 Introduction↩︎

Recently, substantial attention has been devoted to beyond-mean-field effects in ultracold gases. These effects arise from quantum fluctuations (QFs), which can dynamically stabilise an otherwise collapsing Bose-Einstein condensate (BEC). For example, in binary Bose mixtures with weak residual attractive mean-field interactions, QFs can arrest collapse [1], [2]. In three dimensions, the QF contribution is described by the Lee-Huang-Yang (LHY) correction and takes the form of a repulsive term \(\propto n^{5/2}\), where \(n\) is the condensate density. The balance between residual mean-field attraction and QF-induced repulsion produces new localised states, quantum droplets which have been observed in binary BECs [3], [4] and in single-component dipolar BECs [5].

Beyond the study of QFs-induced stationary states, it is important to analyse non-stationary coherent processes, notably the Josephson effect in quantum gases. In bosonic Josephson junctions, phenomena such as macroscopic quantum tunnelling, self-trapping, bifurcations, and homoclinic structures are predicted and have been explored both theoretically and experimentally [6]–[10].

Because macroscopic quantum tunnelling and self-trapping are sensitive to the nonlinear terms in the extended Gross-Pitaevskii equation, QFs are expected to significantly affect dynamical processes in atomic Josephson junctions. To date, only a few studies have examined beyond-mean-field effects in such junctions. In the work [11], a two-mode model for one-component stable condensate with repulsive mean-field interactions and the LHY correction considered as a small perturbation is investigated.

The extended Gross-Pitaevskii equation with an LHY term in a double-well potential has also been investigated: for quasi-one-dimensional Bose-Bose mixtures [12], [13] and for elongated geometries [14]. Stationary states in a four-well potential for a two-dimensional BEC with LHY corrections were obtained in [15]. The Josephson effect between droplets has been explored for two interacting three-dimensional droplets in a symmetric Bose-Bose mixture [16], where the interplay between Josephson oscillations and droplet motion was treated under a weak-overlap approximation. Josephson dynamics of quantum droplets in a quasi-one-dimensional two-core trap were analysed in [17].

Understanding Josephson oscillations between droplets is also important for measuring the superfluid fraction in dipolar-BEC supersolids [18], [19]. Here, we study macroscopic quantum tunnelling and self-trapping of matter waves in a two-dimensional two-core trap while accounting for beyond mean-field effects. We consider two types of states: homogeneous condensates and quantum droplets, and we analyse Josephson oscillations and self-trapping in both settings.

The paper is organised as follows. Section 2 introduces the model. Sections 3 and 4 analyse the dimer model and its bifurcations. Droplet interactions are investigated in Section 5. The Josephson effect for moving droplets is discussed in Section 6.

2 The model↩︎

We consider the dynamics of a two-component BEC loaded in the dual-core cigar-type trap, when the beyond mean-field effects, describing quantum fluctuations, are taken into account. The case when the wave functions of each component inside each core are equal is considered. In this case, the system is described by two coupled modified Gross-Pitaevskii equations of the form [1], [2], [20]: \[\begin{align} & i\hbar\Psi_{1,T}=-\frac{\hbar^2}{2m_0}\nabla^2\Psi_1 +V(\rho)\Psi_1 + \nonumber \\ & \bar g|\Psi_1|^2\Psi_1 \ln(\frac{|\Psi_1|^2}{\sqrt{e}n_0})-\bar K\Psi_2,\nonumber \\ & i\hbar\Psi_{2,T}=-\frac{\hbar^2}{2m_0}\nabla^2\Psi_2 +V(\rho)\Psi_2 + \nonumber \\ & \bar g|\Psi_2|^2\Psi_2 \ln(\frac{|\Psi_2|^2}{\sqrt{e}n_0}) -\bar K\Psi_1 . \end{align}\] where \(m_0\) is atomic mass, \(\bar g=8\pi\hbar\omega_{\perp}/\ln^{2}(a_{12}/a)\) with \(\omega_{\perp}\) denoting the transverse confinement frequency, and where \(a_{11}=a_{22}\equiv a\) and \(a_{12}\) represent the intra- and inter-species \(s\)-wave scattering lengths, respectively. The scattering lengths here are defined in 2D, and their relation to the corresponding 3D scattering lengths is given by \[a_{ij}^{(2D)} = \left(\frac{4\pi}{B}\right)^{1/2} l_0 \exp\!\left(-\gamma - \sqrt{\tfrac{\pi}{2}}\,\frac{l_0}{a_{ij}^{(3D)}}\right),\] where \(B=0.915\), \(l_0 = (\hbar/m_0 \omega_{\perp})^{1/2}\) is the harmonic oscillator length, and \(\gamma\) is the Euler-Mascheroni constant (see Ref. [2]). The scattering length is expressed as \(a_{ij}\) with \(a_{12}=a_{21}\), \(i, j=1\) and \(2\).

The system in dimensionless form is: \[\begin{align} &iu_t + \frac{1}{2}\nabla^2 u - g\ln(|u|^2) |u|^2 u + K v=0, \nonumber \\ &iv_t +\frac{1}{2}\nabla^2 v -g\ln(|v|^2)|v|^2 v + K u=0, \label{eq:coupGPE} \end{align}\tag{1}\] where \[g=\frac{8\pi\sqrt{e}n_0}{\ln^2(a_{12}/a)}, \;n_0 = \frac{e^{-2\gamma -3/2}}{2\pi}\frac{\ln(a_{12}/a)}{a a_{12}}, \;K=\frac{\bar K}{\hbar \omega_{\perp}}.\]

Related work that incorporates quantum fluctuations in a setting closely aligned with our model has been carried out in one dimension, see Ref. [20]. In Ref. [11], the role of dimensionality in atomic Josephson junctions with quantum fluctuations is analysed. The authors employ a two-mode approximation for a two-dimensional Bose-Einstein condensate with repulsive mean-field interactions, supplemented by a weak (perturbative) LHY correction. Owing to the logarithmic renormalisation of the effective coupling in 2D, quantum fluctuations enhance the interaction energy, leading to an upshift of the Josephson oscillation frequency and a lowering of the self-trapping threshold (i.e., facilitating the onset of macroscopic quantum self-trapping).

3 Dimer limit↩︎

Let us consider the case of homogeneous in the space fields i.e.\(u(x,t), v(x,t)\).

Then we have the LHY dimer equations \[\begin{align} iu_t =g|u|^2 u \ln(|u|^2) - K v,\nonumber\\ iv_t = g|v|^2v\ln(|v|^2) - K u. \label{eq:LHYdimer} \end{align}\tag{2}\] Representing fields in the form \[u=A_1 e^{i\phi_1}, \; v = A_2 e^{i\phi2}.\] and introducing the relative variables- \(Z\) is the population imbalance and the relative phase \(\theta\) according to: \[Z=\frac{N_2-N_1}{N}, \;\theta =\phi_2 -\phi_1, N_i=A_i^2, N = N_1 + N_2,\] we obtain the system of equations: \[\begin{align} \label{eq:dimerZtThetat} &Z_t=2 K \sqrt{1-Z^2}\sin(\theta), \\ \nonumber &\theta_t =\frac{gN}{2}[(1-Z) \ln(\frac{N(1-Z)}{2})-\\ \nonumber &-(1+Z) \ln(\frac{N(1+Z)}{2})] -\frac{2 K Z}{\sqrt{1-Z^2}}\cos(\theta). \end{align}\tag{3}\]

The frequency of the Josephson oscillations:

i) zero-phase mode (\(\theta=0\)): \[\label{eq:Jfzero} \omega_J^{(0)} = 2K \sqrt{1+\frac{gN}{2K} \ln\left(\frac{e N}{2} \right)},\tag{4}\]

ii) \(\pi\)-phase mode(\(\theta=\pi\)): \[\label{eq:Jfpi} \omega_J^{(\pi)} = 2K \sqrt{1-\frac{gN}{2K} \ln\left(\frac{e N}{2} \right)}.\tag{5}\] Considering \(Z\) and \(\theta\) as canonical variables, we can write the Hamilton equations: \[Z_t=-\frac{\partial H}{\partial \theta}, \theta_t=\frac{\partial H}{\partial Z}.\] The Hamiltonian for the system (3 ) is: \[\begin{align} \label{HamiltonianPW} H&=&2 K \sqrt{1-Z^2}\cos\theta -\frac{g N}{4} [ 2 Z^2 \ln(\frac{N}{2 \sqrt{e}}) \nonumber\\ &+& (1-Z)^2 \ln(1-Z) + (1+Z)^2 \ln(1+Z) ]. \end{align}\tag{6}\]

The condition for the transition from the macroscopic quantum tunelling regime to the self-trapping mode is: \[H(Z_0,\theta_0)>\mp2K,\] where \(Z_0\) and \(\theta_0\) are the initial imbalance and phase, and we obtain the following expression for the critical value \[\begin{align} \label{CV} & K_{\mathrm{cr}}=(gN(2Z_0^2\ln{\frac{N}{2\sqrt{e}}}+(1+Z_0)^2\ln{(1+Z_0)}+\nonumber\\ &(1-Z_0)^2\ln{(1-Z_0})) \times \nonumber \\ & (8(\sqrt{1-Z_0^2}\cos{\theta_0}\pm1))^{-1}. \end{align}\tag{7}\]

Figure 1 illustrates the phase portraits derived from the Hamiltonian function given by Eq.(6 ). In Figure 1 (a), the zero-phase mode is analysed for various initial values of the population imbalance \(Z_0\). When the initial imbalance is small, e.g., \(Z_0 = 0.04\), Josephson oscillations are observed between the two condensates (dotted line). As \(Z_0\) increases, the character of the dynamics changes, leading to a regime in which atoms remain predominantly localised, a signature of self-trapping behaviour (dashed line). The critical value separating the Josephson and self-trapping regimes is \(Z_{0c} = 0.553\) (solid line), consistent with Eq. (7 ).

Figure 1 (b) presents the corresponding behaviour in the \(\pi\)-phase regime. A similar transition from Josephson oscillations to self-trapping is observed; however, this transition occurs only when the number of atoms is sufficiently small. In this case, the critical imbalance is \(Z_{0c} = 0.6829\), which also agrees with the prediction of Eq. (7 ).

Figure 1 (c) illustrates the case of \(N=0.85\). For \(Z_0=0.004\), a localisation-like mode appears, but it differs from the self-trapping mode. This regime is known as localisation-revival, where the relative imbalance \(Z(t)\) oscillates almost between \(0\) and \(1\). At \(Z_0=0.5\), a self-trapping mode occurs: the relative phase does not grow indefinitely, but instead oscillates harmonically around the \(\pi\)-phase. For \(Z_0=0.99\), Josephson oscillations are observed; however, these oscillations deviate from being perfectly harmonic.

Figure 2 (a), (b), and (c) show the time evolution of the atomic imbalance \(Z(t)\) for the three initial imbalances \(Z_0\) used in Figure 1 (a). Both numerical simulations (solid line) and theoretical predictions are presented (dotted line), demonstrating excellent agreement. The theoretical curves are obtained from the analytical solutions of Eqs. (3 ). A similar dynamical behaviour is observed in the \(\pi\)-phase mode. In Fig. 2 (d), the phase dynamics corresponding to these values of the initial atomic imbalance are presented. The solid line represents the Josephson oscillations at \(Z_0=0.04\), the dashed line indicates the critical trajectory at \(Z_0=0.0553\), and the dotted line denotes the self-trapping mode at \(Z_0=0.065\).

Figure 3 (a) illustrates the dynamics of the atomic population imbalance and the relative phase. The population imbalance oscillates between \(0\) and \(1\), while the relative phase exhibits oscillations around its initial value \(\theta_0\). This regime corresponds to the localisation revival regime. All particles transfer from the first component to the second, then return to the first until the populations equalise, and this cycle repeats periodically, see Fig. 3 (b). Unlike in the self-trapping mode, in the localisation revival regime, the value of the relative phase does not diverge to infinity. This nonlinear oscillation mode was observed in a one-dimensional ultradilute quantum liquid in a double-well potential [13], and we also observed it for a two-dimensional two-core BEC.

Figure 4 shows a graph of the Josephson frequency versus the number of atoms based on Eqs. (4 ) and (5 ). It can be seen from the figures that in the zero-phase mode, Josephson oscillations exist, and as the number of atoms increases, the Josephson frequency also increases. However, for the \(\pi\)-phase regime, Josephson oscillations may exist in a narrow interval (\(N<0.76\)). When \(N\geq0.76\), the Josephson frequency becomes complex, which corresponds to the self-trapping regime. In both panels, the VA and numerical results are in good agreement.

4 Bifurcation↩︎

The fixed points \((Z^*, \theta^*)\) of Eq. (3 ) satisfy \(Z^* = 0\) and \(\theta^* = 0\). Setting \(Z^* = 0\) yields \(\theta^* = n\pi\), or alternatively \(Z^* = \pm 1\). However, \(Z^* = \pm 1\) leads to a singularity in \(\theta_t\), we discard this solution.

From the condition \(\theta_t = 0\), we obtain \[\begin{align} &F(Z) = (1 - Z) \ln\left(\frac{N(1 - Z)}{2}\right) - \\ \nonumber & (1 + Z) \ln\left(\frac{N(1 + Z)}{2}\right)- \cfrac{2K Z \sigma}{\sqrt{1 - Z^{2}}}, \label{eq:FZ} \end{align}\tag{8}\] where \(\sigma = \cos(\theta^*) = \cos(n\pi) = (-1)^n = \pm 1\).

Evaluating the Jacobian matrix at the fixed point \((Z^*, \theta^*)\), we find \[J^* = \begin{pmatrix} 0 & 2K \sigma \sqrt{1 - Z^{*2}} \\ F_Z(Z^*)& 0 \end{pmatrix}. \label{eq:Jacob}\tag{9}\]

The trace and determinant of the Jacobian matrix are given by \(\mathrm{Tr}(J^*) = 0,\) \[\mathrm{det}(J^*) = -2K \sigma \sqrt{1 - Z^{*2}} F_Z(Z^*). \label{eq:detJ}\tag{10}\]

The sign of \(\mathrm{det}(J^*)\) determines the type and stability of the fixed point. The eigenvalues are purely imaginary when \(\mathrm{det}(J^*) > 0\); in this case, the fixed point is a center, which is neutrally stable, see Fig. 1 (a). When \(\mathrm{det}(J^*) < 0\), the eigenvalues are real and have opposite signs, indicating a saddle point, which is unstable, see Fig. 1 (b). The case \(\mathrm{det}(J^*) = 0\) corresponds to a bifurcation.

The bifurcation condition \(\mathrm{det}(J^*) = 0\) leads to the critical particle number \[N_{\mathrm{b}} = -\cfrac{2K \sigma}{g\, W_k \left(-\cfrac{eK \sigma}{g}\right)}, \label{eq:Ncr}\tag{11}\] here, \(W_k(\cdot)\) denotes the Lambert \(W\) function, where the subscripts \(k=0\) and \(k=-1\) refer to the principal and lower real branches, respectively. Since the real branches of the Lambert \(W\) function exist only for arguments \(\ge -\tfrac{1}{e}\), this requires \[\frac{\sigma K}{g} \le e^{-2}. \label{eq:Kg}\tag{12}\]

It follows from Eq. (11 ) that at \(\sigma=1\) the zero-phase admits two real branches of the Lambert \(W\) function, whereas the \(\pi\)-phase admits only one. Consequently, there are two bifurcation points in the zero-phase and a single bifurcation point in the \(\pi\)-phase. The corresponding bifurcation diagrams, obtained by numerically solving Eq. (3 ), are shown in Fig. 5 (a) and (b) for the zero-phase, and in Fig. 5 (c) for the \(\pi\)-phase.

The figure 5 clearly indicates that all bifurcations in the system are of the pitchfork type. Noticing that \(F(-Z)=-F(Z)\) is odd, the equilibrium \(Z=0\) exists for all \(N>0\) and is the unique solution preserving the \(Z\to -Z\) symmetry. In the following, we analyze the bifurcation diagrams for the zero-phase and \(\pi\)-phase separately.

At the zero-phase, i.e., for \(\sigma = 1\), there are two bifurcation points as discussed above. For the chosen parameters \(g = 1\) and \(K = 0.01\), the first bifurcation occurs at \(N_{b,1} \simeq 0.0038\), which is marked in Fig. 5 (a). This corresponds to a super-critical pitchfork bifurcation, in which the stable equilibrium at \(Z^* = 0\) loses stability and gives rise to two stable symmetry-broken equilibria at \(\pm Z^*\). In the bifurcation diagrams, stable and unstable branches are denoted by “S” and “U”, respectively. Beyond this critical value, a second bifurcation of sub-critical type occurs. In this case, the bifurcation point corresponds to \(N_{b,2} \simeq 0.7155\), as shown in Fig. 4(b). Here, the previously existing unbalanced equilibria \(\pm Z^*\) (which are now unstable) collide with \(Z^* = 0\) and disappear, thereby restoring the stability of the origin. In the interval \(N_{b,2} < N < N_s \simeq 0.8665\), different distinct stable equilibria coexist: the symmetric point \(Z^* = 0\) and a pair of large-imbalance fixed points. Consequently, the long-term state reached as \(t \to \infty\) depends on the initial condition \(Z_0\); whichever basin of attraction contains \(Z_0\) determines the outcome. This multistability permits discontinuous jumps and a hysteresis loop when the control parameter \(N\) is varied slowly. Let us start the system in the symmetric state \(Z^{*} = 0\) at some \(N > N_s\) and decrease \(N\) quasi‑statically. The trajectory remains at the origin until \(N\) reaches \(N_{b,2}\), at which point the trivial equilibrium loses stability. An infinitesimal perturbation then sends the system to one of the symmetry‑broken, large‑imbalance branches. Further reduction of \(N\) simply carries the state along this branch. If \(N\) is now increased again, the solution persists on the large‑imbalance branch even after \(N\) rises past \(N_{b,2}\). Only when \(N\) reaches \(N_s\) does the large‑imbalance branch terminate in a saddle‑node bifurcation, forcing the system to jump back to the origin. This path‑dependence, where the response to a cyclic variation of \(N\) is irreversible, is precisely what is known as hysteresis [21], see Fig. 5 (b).

Figure 5 (c) depicts the bifurcation diagram in the \(\pi\)‑phase, which corresponds to the subcritical pitchfork bifurcation. Here, the bifurcation occurs at \(N_{b} \simeq 0.7555\), corresponding to the principal branch of the Lambert \(W\) function for \(\sigma = -1\) in Eq. (11 ).

Figure 6 plots the bifurcation points (in terms of \(N\)) against the coupling constant \(K\), see Eq. (11 ). The solid red curve denotes the \(\pi\)-phase; the dashed blue curve represents the principal branch of the Lambert \(W\) function in the zero-phase; and the dotted blue markers indicate its lower branch.

In the zero-phase, bifurcations occur only within finite regions of the \((N, K)\)-plane, with the turning point at \(K = 0.135335\), see Eq. (12 ). By contrast, in the \(\pi\)-phase, bifurcations persist for all \(K\). One can analyse the bifurcations in the system by its normal form. To obtain the normal form, we expand Eq. (3 ) in a Taylor series about \(Z^*=0\). Carrying the expansion out to the fifth order, one finds \[F(Z)\simeq \alpha(N) Z +\beta(N) Z^3 +\gamma(N) Z^5+\mathcal{O}(Z^7),\] with \[\begin{align} \nonumber \alpha(N) &= -2K\sigma - g \, N\, \ln\!\bigl(eN/2\bigr),\\ \nonumber \beta(N) &= \tfrac{g\,N}{6} - K\sigma,\\ \gamma(N) &= \tfrac{g\,N}{20} - \tfrac{3K}{4}. \end{align}\] Retaining only the leading-order terms, \(F(Z) \approx \alpha(N) Z +\beta(N) Z^3\), recovers the canonical normal form of a pitchfork bifurcation [21]. In the following section, we examine Josephson oscillations between quantum droplets.

5 Droplets interaction↩︎

We employ the variational approximation (VA) to investigate macroscopic quantum tunnelling between quantum droplets. The corresponding Lagrangian density for the system described by Eq. (1 ) is given by: \[\begin{align} \label{eq:lagr} &\mathcal{L} = \frac{i}{2}(uu_t^{\ast} - u^{\ast}u_t) +\frac{1}{2}|\nabla u|^2 +\frac{g}{2}|u|^4\ln(|u|^2) - \nonumber\\ &\frac{g|u|^4}{4}-K(u^{\ast}v + u v^{\ast}) + \frac{i}{2}(vv_t^{\ast} - v^{\ast}v_t)+ \nonumber\\ &\frac{1}{2}|\nabla v|^2 +\frac{g}{2}|v|^4\ln(|v|^2)-\frac{g|v|^4}{4}. \end{align}\tag{13}\] Following the idea of the work [22], let us investigate the Josephson oscillations and self-trapping regimes, using the VA to the QD dynamics. The VA ansatz for the wavefunctions is [23]–[25]: \[\begin{align} \label{eq:ansatz} u=A_1 \exp\left(-\frac{1}{2} (ar)^{2m} + i\phi_1 \right), \nonumber\\ v = A_2 \exp\left(-\frac{1}{2}(ar)^{2m} +i\phi_2\right). \end{align}\tag{14}\] The averaged Lagrangian is given by \(L = \int_{-\infty}^{\infty} \mathcal{L} dxdy\). Substituting the ansatz (14 ) into Eq.(13 ) and integrating, we found the averaged Lagrangian:

\[\begin{align} L &=& \frac{\pi }{a^2} \left\{ \sum_{i=1}^{2}\left[\Gamma(1+M)A_i^2\phi_{i,t}^2 + \frac{1}{2M}a^2A_i^2 \right.\right. \nonumber \\ && \left.\left. - 2^{-2-M}g \Gamma(1+M)((1+M)A_i^4 - 2A_i^4\ln A_i^4)\right] \right. \nonumber\\ && \left. -2K A_1 A_2\cos(\phi_2-\phi_1) \Gamma(1+M) \right\}, \end{align}\]

where \(M=1/m\). The Euler-Lagrange equations give: \[\begin{align} \label{eq:QD95JosepEq} &\theta_t=-\frac{2KZ}{\sqrt{1-Z^2}}\cos{\theta}+\frac{2^{-1-M}a^2gN}{\pi\Gamma(1+M)} -\left[ZM \nonumber \right.\\ &(1+Z)\ln{\left(\frac{N(1+Z)}{2}\right)}+(1-Z)\ln{\left(\frac{N(1-Z)}{2}\right)} \nonumber\\ &\left.-2Z\ln{\left(\frac{a^2}{\pi \Gamma(1+M)}\right)}\right] \nonumber \\ &Z_t=2K\sqrt{1-Z^2}\sin{\theta}, \end{align}\tag{15}\] where \(\theta=\phi_2-\phi_1\) is relative phase, \(Z=(N_2-N_1)/N\) is the atomic imbalance, \(N_1=2\pi\int_{0}^{\infty}|u|^2dr=2^{-2-M}\pi g A_1^4 \Gamma(1+M)/a^2\), \(N_2=2\pi\int_{0}^{\infty}|v|^2dr = 2^{-2-M}\pi g A_2^4 \Gamma(1+M)/a^2\) are the number of atoms in the cores and \(N\) is the total number of atoms. We can write the Hamiltonian function for these new variables as: \[\begin{align} \label{Hamiltonian} & H=2K\sqrt{1-Z^2}\cos{\theta}+\frac{2^{-2-M}a^2gNZ^2}{\pi \Gamma(M)}\left[1+\frac{1}{M}- \right. \nonumber\\ & \left.\frac{2}{M} \ln{\left(\frac{a^2}{\pi \Gamma(1+M)}\right)}\right]-\frac{2^{-2-M}a^2gN}{\pi \Gamma(1+M)}\left[(1-Z)^2\times \right. \nonumber\\ &\left. \ln{\left(\frac{N(1-Z)}{2}\right)}+(1+Z)^2\ln{\left(\frac{N(1+Z)}{2}\right)} \right]. \end{align}\tag{16}\] Using the system of equations (15 ), it is possible to find analytical expressions for Josephson frequencies for small amplitude oscillations. a) zero-phase mode: \[\begin{align} \label{JFzero} \omega_J=2K \left\{1-\frac{2^{-M}a^2gN}{2K\pi M \Gamma(M)}\times \right. \nonumber\\ \left.\left[\frac{M}{2}-\ln{\left(\frac{e a^2N}{2\pi M \Gamma(M)}\right)}\right]\right \}^{1/2}, \end{align}\tag{17}\] b) \(\pi\) -phase mode: \[\begin{align} \label{JFpi} \omega_J=2K \left\{1+\frac{2^{-M}a^2gN}{2K\pi M \Gamma(M)}\times \right. \nonumber\\ \left.\left[\frac{M}{2}-\ln{\left(\frac{e a^2N}{2\pi M \Gamma(M)}\right)}\right]\right \}^{1/2}. \end{align}\tag{18}\]

We can derive the self-trapping condition using Eq. (16 ): \[\begin{align} H(Z_0,\theta_0)>\pm 2K-\frac{2^{-1-M}a^2gN\ln{(N/2)}}{\pi \Gamma(1+M)}, \end{align}\] and can be found \(K_{\mathrm{cr}}\): \[\begin{align} \label{critical95v} &K_{\mathrm{cr}}=\frac{2^{-3-M}Ma^2gN\ln{(N/2)}}{\pi \Gamma(1+M)(\pm 1-\sqrt{1-Z_0^2}\cos{\theta_0}}\left\{\frac{2}{M}\ln{\left(\frac{N}{2}\right)}+\right. \nonumber\\ & Z_0^2-\frac{2Z_0^2}{M}\ln{\left(\frac{a^2}{\sqrt{e}\pi \Gamma(1+M)}\right)}-\frac{1}{M} \left[(1-Z_0)^2 \right. \nonumber \\ & \ln{\left(\frac{N(1-Z_0)}{2}\right)}+(1+Z_0) \ln{\left(\frac{N(1+Z_0)}{2}\right)} \left ] \left\}\right. \right. , \end{align}\tag{19}\] here, the \(\pm\) signs correspond to the zero- and \(\pi\) -phase modes, respectively. In the symmetric case, Eq. (1 ) can be rewritten as \[i\,\partial_{t} u + \frac{1}{2} \nabla^{2} u - g\,|u|^{2} \log\left(|u|^{2}\right) u + K u = 0.\] The linear coupling term can be eliminated by applying the transformation \(u(\mathbf{r},t) = \Psi(\mathbf{r},t) \exp(-i K t)\). Here, the new wave function \(\Psi\) inherits exactly the same modulus, i.e., \(|\Psi| = |u|\), so the nonlinear term remains unchanged. The stationary and dynamical properties of the resulting equation, analysed within the VA using a super-Gaussian ansatz, have been thoroughly studied in Ref. [25]. Therefore, in our subsequent theoretical analysis in this section, as well as in solving Eq. (15 ), we implement the variational results obtained in Ref. [25] to determine the stationary parameters of quantum droplets. It is worth noting that, in the symmetric case, the VA provides highly accurate predictions for the stationary parameters of quantum droplets, in very good agreement with numerical simulations.

Figure 7 (a) shows the dependence of the chemical potential on the number of atoms \(N\), obtained from imaginary-time simulations of Eq. (1 ) with an initial imbalance \(Z_0 = 0.1\) between the components. The red solid line represents the chemical potential \(\mu_1(N_1)\) of the first component \(u_1\), while the blue dashed line represents \(\mu_2(N_2)\) of the second component \(u_2\). In the numerical computations, the individual component norms are defined as \(N_1 = N(1 - Z_0)/2\) and \(N_2 = N(1 + Z_0)/2\), respectively. The bottom black solid line illustrates the total chemical potential \(\mu = \mu_1 + \mu_2\) as a function of the total norm \(N = N_1 + N_2\). According to the Vakhitov-Kolokolov stability criterion, the negative sign of the derivative \(d\mu/dN < 0\) corresponds to dynamically stable quantum droplets. This stability was further confirmed by direct numerical simulations under small perturbations.

Figure 7 (b) presents the dependence of parameters of stationary quantum droplets, namely, amplitudes, widths, and super-Gaussian indices, on the number of particles in each component, for an initial imbalance \(Z_0 = 0.1\). The red solid line and the blue dashed line correspond to the first and second components, respectively. As the number of particles increases, the amplitudes of both components saturate, while their widths begin to grow noticeably only at higher norms. The increase in the super-Gaussian index \(m\) with particle number indicates a transition of the droplet profile from a Gaussian to a flat-top shape. This behaviour highlights the emergence of an incompressible fluid-like nature in the quantum droplets (see also Fig. 8), consistent with the findings of Ref. [2].

Figure 8 displays the density profiles of both components for various values of total norm \(N\), assuming \(Z_0 = 0.1\). For small \(N\), the droplet profiles are approximately Gaussian, whereas for larger \(N\), the profiles develop flat-top structures, further supporting the identification of the droplets as incompressible quantum fluids.

Figure 9 (a) displays the Josephson frequency as a function of the norm for the zero-phase regime. The dashed curve represents the theoretical results obtained from Eq. (17 ), while the points joined by a line represent the numerical simulation results. As shown in the figure, the theoretical and numerical results agree well for a small number of atoms. As the number of atoms increases, the overall trend remains similar, but the actual values begin to differ significantly. At certain values of the number of atoms, the Josephson frequency obtained from Eq. (17 ) becomes a complex number, leading to a discontinuity in the frequency curve. In theory, the part of the Josephson frequency where it is a complex number is represented as zero on the graph for these specific frequency values. According to theoretical results, this gap appears within the norm range of approximately \(2.25\) to \(23.70\), while numerical simulations indicate a broader interval from \(2.25\) to \(86.54\). This range of numbers of atoms corresponds to the self-trapping regime. However, this interval was derived for fixed values of the parameters \(g\) and \(K\). As can be seen from Eq. (17 ), within this interval of \(N\), where self-trapping occurs, the Josephson frequency becomes complex. Nevertheless, by appropriately tuning the parameters \(K\) and \(g\), it is still possible to observe Josephson oscillations in this interval of \(N\). A similar behaviour is found in the interval obtained from the numerical simulations, where Josephson oscillations are also present, provided that the parameters \(K\) and \(g\) are chosen appropriately.

Figure 9 (b) corresponds to the \(\pi\)-phase regime. In this case, the results obtained from the VA agree well with the numerical results only when the number of atoms is sufficiently small (\(N<16\)). Numerical simulations reveal that a self-trapping regime occurs within the range of the number of atoms from \(2.85\) to \(7.22\). When \(N>7.22\), quantum droplets, after a few oscillation periods up to \(5-7\), move away from each other. This separation is attributed to the interaction energy becoming repulsive, which occurs because it is proportional to \(\cos{\theta}\) (see Eq. (20 )), where \(\theta\) is the relative phase of QDs.

For large particle numbers \(N\), the results of the VA and numerical simulations show noticeable discrepancies. The primary reason is that, in deriving the Josephson equations, we assumed identical super-Gaussian parameters for both quantum droplets, i.e., \(m_1=m_2=m\), as well as equal widths, \(w_1=w_2=1/a\), which were considered time-independent, while only the amplitudes were allowed to vary in time. In reality, as the number of atoms increases, the super-Gaussian index \(m\) and the droplet width \(w=1/a\) both undergo significant changes, whereas the amplitude remains almost unaffected. Consequently, the discrepancy between the VA and numerical results becomes more pronounced as \(N\) increases. If one were to incorporate the time dependence of the super-Gaussian parameters and the droplet width into the VA, the derivation of the Josephson equations would become mathematically intractable.

Figure 10 displays the time evolution of the atomic population imbalance \(Z(t)\) and the relative phase \(\theta(t)\). It can be seen from the figure that for small norms, \(N=1\), the results obtained through numerical simulations and the VA are in good agreement.

In Fig. 11, the dynamics of the relative atomic imbalance \(Z(t)\) and the relative phase \(\theta(t)\) are shown for different norms and phase regimes. For \(N=100\), Josephson oscillations are observed (see the middle solid and dash-dotted lines in Fig. 11 (a)). In the \(\pi\)-phase regime, the droplets separate after a few Josephson oscillations, whereas in the zero-phase regime, the long-time dynamics of \(Z(t)\) indicate that the droplets remain bound. The upper red dashed line and the lower gray dotted line in Fig. 11 (a) correspond to \(N=40\) and \(N=5\), respectively, both representing the self-trapping regime. The associated dynamics of \(\theta(t)\) are displayed in panel (b) with matching line styles: in the Josephson oscillation regime, \(\theta(t)\) exhibits harmonic oscillations (left axis), while in the self-trapping regime, the relative phase increases monotonically without bound (right axis). The corresponding density profiles are shown in Figure 12. These observations in the dynamics of \(Z(t)\) and \(\theta(t)\) are further supported by the frequency spectra shown in Fig. 9.

Figure 13 shows the dependence of the critical coupling \(K_{\mathrm{cr}}\) on the norm \(N\) for the zero-phase mode. The points represent the theoretical prediction from Eq. (19 ), while the solid curve denotes the results of numerical simulations. For small norms (\(N \lesssim 20\)), the agreement is reasonable. For larger norms (\(N>20\)), however, a clear discrepancy emerges: the numerically obtained values of \(K_{\mathrm{cr}}\) approach zero, whereas the theoretical prediction becomes negative. Hence, the two approaches agree only qualitatively. In the simulations, \(K_{\mathrm{cr}}\) is determined with a numerical precision of \(10^{-4}\). Notably, dynamics are self-trapped for \(K>K_{\mathrm{cr}}(N)\), while Josephson oscillations occur for \(K<K_{\mathrm{cr}}(N)\).

6 Quantum droplets moving in two channels↩︎

In this section, we analyze moving of the QDs in two channels. There is an effective interaction energy between the QDs, which can be repulsive or attractive depending on the relative phase \(\theta\). To qualitatively estimate this effective energy, we choose the solutions of Eq. (1 ) as simple Gaussian functions:

\[\begin{align} & u=A_1 e^{-((x-x_1)^2+ (y-y_1)^2)/2a_1^2 + i\phi_1} \nonumber \\ & v=A_1 e^{-((x-x_2)^2+ (y-y_2)^2)/2a_2^2 + i\phi_2} \end{align}\] Then the energy of interaction can be calculated from the overlap integral \[E_{int} = K\int\int dx dy (u^{\ast}v + v^{\ast}u)\] Substituting the expressions for the trial functions \(u_i\) and calculating integrals, we obtain \[\label{eq:int} E_{int} = A_1A_2 \frac{\pi}{a_1^2 + a_2^2}e^{-\alpha^2 \Delta R^2}\cos(\theta),\tag{20}\] where \[\alpha^2 = \frac{a_1^2a_2^2}{2(a_1^2 + a_2^2 )}, \Delta R^2 = (x_1-x_2)^2 + (y_1-y_2)^2, \theta =\phi_1 - \phi_2.\] The character of interaction between droplets depends on the value of their relative phase difference \(\theta\). For \(\theta =0\), we have the attractive interaction; for \(\theta =\pi\), the interaction is repulsive.

Although the expression given in Eq. (20 ) is formally valid only for small values of \(N\), it qualitatively captures the behaviour of the interaction energy even for larger \(N\). Figure 14 (a) illustrates the time evolution of the relative displacement of the centre of mass between two droplets initially confined to separate channels. The solid curve corresponds to the zero-phase regime, while the dashed curve represents the \(\pi\)-phase regime. The displacement is shown only along the \(x\)-axis; however, a similar dynamical behaviour is observed along the \(y\)-axis. As evident from the figure, in the zero-phase regime, the droplets undergo a brief oscillation before moving apart. This behaviour is attributed to the evolution of the relative phase \(\theta(t)\). Since \(\theta(t)\) is non-linearly coupled to the atomic imbalance \(Z(t)\) (see Eq.(15 )), variations in \(Z(t)\) lead to the corresponding changes in \(\theta(t)\), as shown in Fig. 14 (b). As \(\theta(t)\) increases, the interaction energy decreases according to Eq. (20 ), which weakens the attractive force between the droplets and causes them to separate. In contrast, in the \(\pi\)-phase regime, the initial interaction is repulsive (see Eq. (20 )), and the droplets immediately repel each other without exhibiting oscillatory motion (see Fig. 14 (a) dashed line). A similar behaviour is observed for the initial relative phase of \(0.25\pi\) and \(0.5\pi\). Figures 14 (a) and (b) correspond to the case where the interaction parameter is set to \(K = 0.01\), under which no Josephson oscillations are observed.

Figures 14 (c) and (d) correspond to the case with \(K = 0.02\) and are presented for the zero-phase mode. In this case, the interaction between the droplets is stronger, and since the relative phase oscillates around zero (see Fig. 14 (d)), the centres of mass of the droplets exhibit oscillations relative to each other over time. As a result, Josephson transitions become apparent (see Fig. 14 (d)). However, because the distance between the droplets oscillates in time, the resulting Josephson oscillations are not perfectly harmonic.

The Andreev-Bashkin effect refers to a nondissipative (entrainment) drag between two superfluids, where a superflow in one component induces a mass current in the other [26], [27]. In our system, the Andreev-Bashkin effect may also occur: if one quantum droplet is set into motion, under certain conditions, it can drag the second droplet, leading to their joint motion. The Andreev-Bashkin effect is illustrated in Figure 15, where the flow of one quantum fluid induces another to flow in the same direction without viscosity. At the initial time (\(t=0\)), the droplet of the second component is given a small kick (phase gradient \(k=0.02\) along the positive \(x\)-axis. During the subsequent evolution, the kicked droplet drags the first-component droplet, which has no initial velocity, so that both move together. The time evolutions of their centres of mass and the relative phase are shown in Fig. 15 (a).

While moving together, the droplets also undergo Josephson tunnelling between components, as seen in Fig. 15 (b). This behaviour is consistent with the presence of the Andreev-Bashkin effect in our system.

7 Conclusion↩︎

We investigate a binary Bose gas in a symmetric, dual-core, pancake-shaped trap. The system is modelled by two linearly coupled two-dimensional Gross-Pitaevskii equations that include Lee-Huang-Yang corrections. We begin with a spatially uniform condensate, derive the Hamiltonian, and obtain the Josephson oscillation frequencies in the zero- and \(\pi\)-phase modes. For each phase, we analyse the phase portrait, distinguishing finite and infinite trajectories and the separatrix that divides them. Finite trajectories correspond to Josephson oscillations, while infinite ones reflect self-localisation. For both the zero- and \(\pi\)-phase modes we derive, in closed form, the critical linear coupling (or, equivalently, the critical initial imbalance) that places the system on the separatrix. We also determine the parameter conditions set by the strength of quantum fluctuation interactions and the linear coupling under which bifurcations occur. In the zero-phase, two pitchfork bifurcations appear in sequence: a supercritical bifurcation at smaller norms followed by a subcritical one at larger norms. The resulting bifurcation diagram exhibits bistability and hysteresis. In the \(\pi\)-phase, we find a single subcritical pitchfork.

We then turn to quantum droplets and develop a variational approach for their coupled dynamics. Stationary droplet parameters are obtained via the variational method and imaginary-time simulations, and we verify their stability in direct numerical simulations. The VA predicted profiles match the numerically computed QDs with high accuracy. We also derive analytic expressions for the Josephson frequencies in both the zero- and \(\pi\)-phase droplet regimes and confirmed these against numerical simulations. For smaller particle numbers, the variational approach results for the Josephson oscillations agree well with numerics. For larger particle numbers, the values of the predicted frequencies can differ by orders of magnitude, but the qualitative trends remain consistent, so we rely primarily on numerical simulations in that regime.

In the \(\pi\) phase for quantum droplet tunnelling, we observe that the two droplets split in opposite directions after roughly \(6-7\) oscillation periods. This splitting is explained by evaluating the droplet-droplet interaction energies. In numerical simulations, we map the critical linear coupling that separates Josephson and self-trapped regimes as the particle number varies. Finally, in the two-core geometry, we examine inter-droplet interactions and find numerical signatures consistent with the Andreev-Bashkin superfluid drag effect.

8 Acknowledgments↩︎

This work has been supported from the State Budget of the Republic of Uzbekistan (Grant No. 2025 year award).

References↩︎