Manifestation of relative phase in dynamics of two interacting Bose-Bose droplets

We study coherent dynamics of two interacting Bose-Bose droplets by means of the extended Gross-Pitaevskii equation. The relative motion of the droplets couples to the phases of their components. The dynamics can be understood in terms of the evolution of zero-energy modes recovering symmetries spontaneously broken by the mean-field solution. These are translational symmetry and two U(1) symmetries, associated with the phases of the droplets' two components. A phase-dependent interaction potential and double Josephson-junction equations are introduced to explain the observed variety of different scenarios of collision. We show that the evolution of the droplets is a macroscopic manifestation of the hidden dynamics of their phases. The occurrence of nondissipative drag between the two supercurrents (Andreev-Bashkin effect) is mentioned.


I. INTRODUCTION
Quantum droplets are self-bound objects formed by ultracold atoms. Despite having densities about eight orders of magnitude smaller than air they behave like liquids. Droplets were first observed in dipolar gases in systems of dysprosium [1][2][3] or erbium atoms [4]. Their binding mechanism occurred to be the same as predicted by Petrov [5] for two-component Bose-Bose mixtures. Self-bound systems of ultracold atoms are formed when the mean-field energy of the gas almost vanishes and quantum fluctuations become important. The Lee-Huang-Yang [6][7][8][9] contribution to the energy of the system constitutes an essential ingredient of this stabilizing mechanism. Quantum droplets as considered by Petrov [5] were obtained in a mixture of two hyperfine states of 39 K [10][11][12][13]. Experiments confirmed the main predictions of the theory: nonspreading density profiles and values of equilibrium densities.
The lifetime of droplets is limited by three-body losses, but in heteronuclear mixtures of 87 Rb and 41 K [14] it might exceed 100 ms [15], which is much longer than the lifetime of two-component homonuclear mixtures. Therefore experiments with heteronuclear droplets [14] pave a way toward long-lived binary droplets which in turn opens the possibility to study dynamical situations, such as collisions.
In the case of colliding classical droplets, two scenarios are possible, depending on the relative values of kinetic and surface energies: coalescence or splitting [16][17][18]. This was also observed in binary collisions of quantum droplets composed of two hyperfine states of 39 K [19]. However, the dynamics of these processes may be richer. The superfluid character of the colliding objects introduces additional degrees of freedom: the relative phases of their components. The coupling of two superfluid systems invokes an analogy to the Josephson effect. Coherent Josephson oscillations are at the heart of recently observed supersolid behavior in weakly coupled droplet systems in 1D [20][21][22], as well as in two-dimensional geometries [23,24]. A Josephson-junction based approach was used to describe out-of-equilibrium dynamics of a supersolid [25]. For a review of recent advances in the field of quantum droplets see [26][27][28][29][30].
Finally, we want to mention a totally different system, where effects very similar to those studied here were observed. These are coherent oscillations of neutrons during a collision of two different nuclei. This peculiar result was reported recently [31,32]. A minute Josephson oscillating current of Cooper paired neutrons was found in analysis of data from collisions of tin and nickel nuclei 116 Sn + 60 Ni at energies around 150 MeV, barely sufficient to overcome Coulomb repulsion. As was suggested theoretically [33][34][35][36], the phase difference of a pairing field of two initially independent nuclei is responsible for this effect.
In the present paper we study low-energy collisions of two interacting Bose-Bose droplets taking into account the coherent exchange of atoms between them. Our approach is based on numerical integration of the extended Gross-Pitaevskii (eGP) equations supported by analysis of equations of motion for zero-energy (Goldstone) modes of the system. These two-component Josephson-junction equations allow for deeper understanding of the observed dynamics. We predict a host of different possible scenarios during the droplets' approach. This includes in particular a coherent transfer of atoms in the form of direct or alternating Josephson currents, which may result in total evaporation of one of the droplets, or their attraction, merging, or repulsion. We observe dynamics akin to a droplet subject to a ponderomotive force, but also the Andreev-Bashkin effect, i.e., entrainment between the two superfluids.
The paper is organized as follows. In Sec. II we outline the theory of two-component Bose-Bose droplets using eGP equations. Next in Sec. III we discuss the preparation of the initial state. We compare the fully coherent state, where the relative phase between the two droplets is well defined, with two Fock states, where the relative phase is not controlled. In Sec. IV we find an expression for the interaction potential between two separate quantum droplets overlapping with their tails. In Sec. V we use a formalism based on the so-called zero-energy modes Hamiltonian to derive Josephson-junction (JJ) equations for the moving two-component junction. The equations allow us to find effective simplified dynamics of a collision of interacting droplets in terms of their relative positions and velocities, coupled to their relative phases and Josephson currents. In Sec. VI we present results of time-dependent numerical simulations of collisions for small and large droplets assuming different initial phases of the droplets. These simulations are compared to predictions of the Josephson-junction model which allows for a clear physical picture of the observed dynamics. We conclude with Sec. VII.

II. EXTENDED GROSS-PITAEVSKII EQUATIONS
A two-component ultracold Bose-Bose droplet may be described quite accurately by a mean-field energy functional. The energy density of a droplet is where we assume that atomic masses of both components are equal to m. n 1 and n 2 are the atomic densities related to the corresponding wave functions n i = |ψ i | 2 . We assume intraspecies repulsion of strength g ii > 0 and interspecies attraction proportional to g 12 < 0, where g ij = 4π 2 a ij /m and a ij are the s-wave scattering lengths. Moreover attraction slightly dominates over repulsion, δg = g 12 + √ g 11 g 22 < 0. The energy functional accounts for quantum fluctuations given by the Lee-Huang-Yang (LHY) term [6]. For negative values of δg < 0 this contains a small imaginary part. Second-order Beliaev theory for mixtures allows us to account for higher order terms [37,38] which may cure this problem. However, the Beliaev approach for nonuniform two-component mixtures is extremely challenging. Alternatively, a phenomenological approach based on thermodynamic relations was suggested [39], but sound velocities are still imaginary for very small densities na 3

11
(δg/g) 2 in this framework. Another direction is to include some presumably missing small contributions to the system energy. It was shown in [40] that adding a pairing energy removes the imaginary LHY component. Despite all these different efforts no consistent ultimate solution to the problem of imaginary LHY energy is present yet. For the purpose of the present analysis we simply neglect the imaginary term, as is commonly done in theoretical studies of quantum droplets. This omission is justified close to the instability threshold.
The ground state is obtained by minimization of the following energy functional: where the interaction energy density equals E = (|Ψ 1 | 2 , |Ψ 2 | 2 )/µ 0 . The two chemical potentials µ i = µ i (N 1 , N 2 ) are the only free parameters. Dynamical extended Gross-Pitaevskii equations (eGP) consistent with the above energy functional have the form: The eGP equations quite accurately describe large Bose-Bose droplets [10][11][12][13]. In order to get quantitative agreement for small droplets, more sophisticated approaches are needed [42]. A beyond local density approximation to the LHY energy is desirable then. In our description we do not account for three-body losses, dN/dt −K|Ψ i | 6 . They may be accounted for by including an appropriate imaginary term on the right hand side of Eq. (3). As shown in the experiment [19] the effect of this term was important for qualitative agreement between the experimental results and a theory based on the eGP equation. In our case this term would play a similar role. At present, the life-time of homonuclear droplets (∼ 10ms) is shorter than the complete process, i.e. preparation of a suitable initial state and collision. We do not include this effect since our aim is to study which collision scenarios are possible once three-body losses are overcome. We believe that theoretical and experimental progress in the field will allow for such long-lived droplets in the future. Heteronuclear droplets in lower dimension are a possible remedy, but other unexpected solutions can not be ruled out.

III. INITIAL STATE
Collisions of two droplets were studied experimentally in [19]. The initial state was prepared by forming two droplets in a double well potential followed by removal of the barrier separating the two wells. No relative phase was imprinted onto the droplets and the results obtained were consistent with the assumption that their phases are identical. Therefore many aspects of the collisions were similar to classical liquid droplets. For low velocities the two droplets merged while for larger velocities they separated after the collision. The critical velocity was found to depend on the droplet size. The effect of size in collisions of 1D droplets was discussed in [43].
Here we consider dynamics of interacting droplets having initially nonzero relative phases and moving in opposite directions. Such an arrangement can be experimentally achieved in a way similar to those of experiment [19]. After separation however, a desired phases should be optically imprinted onto the droplets [44].
The preparation of an initial state with well controlled initial phases is a challenging task. In the ideal situation, splitting of a single droplet into two symmetric smaller ones ought to be performed adiabatically in order to avoid excitations, which in the case of a droplet would lead to evaporation of particles. Special attention should be given to the symmetric placing of the splitting barrier. This was successfully achieved in experiments with a quantum gas placed in a double well potential and in Josephson-junction studies [45][46][47][48]. Here, in principle, the same procedure might be implemented In the following we assume that the initial state is a superposition of two-component waves moving towards each other (compare Eq. (8)): where i = 1, 2 enumerates components, and φ ini i,L(R) = Φ i,L(R) (r)e ±i(pr−∆φ/2) . Such an initial state dovetails the N -body "superfluid" state: with every atom being in the superposition of the left and right droplet. The phase between the right and left droplet is well controlled, i.e. does not vary from one realization of the system to the other. However, the number of atoms fluctuates. It follows from Eq. (5) that for large N each droplet is in a coherent state and has on average N/2 atoms with a dispersion equal to √ N /2. The one particle density matrix has only one nonzero eigenvalue and a corresponding eigenvector, Eq. (4). This is a situation ideally suited for eGP equations.
On the other hand one may consider a situation where the two droplets are prepared independently. The Nbody state is then a product of two Fock states: where P stands for permutations of the N atom positions. From the point of view of one-particle measurements, the system is fragmented, with two eigenvalues equal to 1/2 each. The two eigenvectors of the onebody density matrix correspond to two wavefunctions: the right and the left one, |φ ini i,L(R) (r)|. The multiconfiguration time-dependent Hartree approach could be a right tool to tackle such a problem [49][50][51][52][53].
The fragmentation of one-particle density is the result of a lack of a well defined relative phase, which strongly fluctuates from one realization of the system to the other. However, in every single realization (for repeated measurements) it has some (uncontrolled) value which is fixed in course of a measurement. This issue was discussed in detail in the context of interference of two Fock states [54,55]. In a real cold-atom experiment the system is monitored by a CCD camera where a shadow of the atomic cloud is captured at some given instant of time. In such a single shot many atoms, ∼ d 1, interact with light and are monitored simultaneously. Therefore all observables ought to be averaged with respect to a d-body density matrix which in case of the product of two equally populated Fock-states Eq. (6) was found in [55]: where the limit N d 1 was assumed. Thus even if the system seems to be fragmented, its dynamics in a single realization may be safely described assuming a superposition of the left and right droplet states, Eq. (7) with some arbitrary phase φ. The phase will differ from one realization to the other but is fixed for all detected atoms in a single measurement. Therefore it is justified to use eGP equations also in the case of independently prepared droplets. Some information about the actual relative phase can be eventually deduced a posteriori from the observed scattering dynamics.

IV. INTERACTION POTENTIAL
At large distances the eGP equations simplify to Coefficients A i and λ i depend on the number of particles of the two species in the droplet. To find the interaction potential we follow the approach presented in [56,57]. We assume that the two droplets are separated by a distance significantly larger then their radii. The exponentially vanishing tails of their wavefunctions overlap and contribute to the interaction energy. The total wavefunctions Ψ Sc i of the system are assumed to be sums of two stationary droplet solutions, a left and right one, having different phases: and separated by a distance R = |r R − r L |: where B i are chosen to ensure normalization of the wavefunction Ψ Sc i to the value of (N 0 i ) L + (N 0 i ) R . The interaction energy of the left and right droplet may be defined as the difference between the energy of two overlapping and two infinitely separated droplets The potential depends not only on distance U i (R) = A 2 i ( 4π R )e −λiR , but also on the phase difference between the droplets, and may be attractive or repulsive. The characteristic range of the interaction is 1/λ i .
The coefficients A i are found from the continuity condition, where R 0 is the radius of a droplet having N 0 i atoms, i . This is an approximate expression valid only for very large droplets, where the bulk density may be approximated by n 1(2) = |Ψ 1(2) | 2 = 1, 1 √ a . In the general case, the coefficients A i are to be determined by fitting to a numerical solution. This is the approach of the present study.
Equation (9) is valid for two droplets with equal number of atoms. For a more general treatment, we compare this analytic result with numerical calculations of scenarios in which there is a small population imbalance between the droplets, to obtain the following formula for the spatial part of the interaction potential: The exponents λ L,R i characterize the exponential tails of the wavefunctions Ψ L,R i .

V. TWO-COMPONENT JOSEPHSON JUNCTION EQUATIONS
Binary droplet collisions may be described by a set of time-dependent 3D partial differential eGP equations -Eqs. (3). Initially the wavefunction is a sum of two stationary solutions Ψ Sc i (r). If the initial kinetic energy of translational motion of the droplets is low we may assume that they move as a whole without internal excitations, preserving their shape. As long as their relative separation is larger then the diameter 2R 0 of a droplet, i.e. if |r L (t) − r R (t)| > 2R 0 we may assume that no Bogoliubov quasiparticles are excited. The exponential density tails overlap forming a weak link. The phase difference between the left and right parts of the two wavefunctions will trigger a coherent flow of atoms between the droplets.
These are the Josephson-junction-like oscillations of particle number and relative phase.
Instead of solving the full set of eGP equations, Eq. (3), the oscillations may be described adequately by only considering the modes involved in the dynamics -these are the zero-energy or Goldstone modes [58,59]. The energy of a free droplet does not depend on a particular choice of phases of its wavefunctions nor on its position in space. The mean field solutions break these continuous symmetries. In consequence, zero-energy modes which recover the broken symmetries appear in the excitation spectrum. Dynamics of zero modes in a two droplet system is given by the Hamiltonian [59]: In above p α = −i ∂ ∂rα are kinetic momenta of left and right droplet, α = L, R, while r α are positions of droplets' centres. Kinetic momenta of hard and soft modes, p α H and p α S , are where: and the plus (minus) sign refers to the soft (hard) mode respectively while δN α i : are deviations of particle numbers N α i from their equilibrium values. The pairs (r α , p α ) and (φ α i , δN α i ) are three sets of canonically conjugate variables. Coefficients M α are masses of the two droplets, M α = (N 0 1 ) α + (N 0 2 ) α while M α H,S are "masses" of the two phase-modes: The upper sign "+" corresponds to the mass of the soft mode while the lower sign "-" represents the hard mode. As shown in [59] the values of M S and M H differ significantly, 1/M α H ∼ 1 1/|M α S | ∝ |δa|/ √ a 11 a 22 . This inequality justifies addressing the corresponding excitations as hard and soft modes. In addition, the mass of the soft mode is negative and changes its sign only for very large droplets. H may be treated as a quantum Hamiltonian [58,[60][61][62], however for the present purpose we restrict ourselves to a classical treatment. This means neglecting quantum fluctuations of φ α i and δN α i . The equations of motion generated by H ensure conservation of particle number of each species and the center of mass momentum. It is convenient to use relative coordinates, R = r R − r L and P = p R − p L . Furthermore we assume that initial orbital angular momentum is equal to zero and the droplets move towards each other along the x-axis, so we may omit vector notation. The equations of motion for the relevant quantities are: where M = M L M R M L +M R is the reduced mass and ∆µ i = µ R i − µ L i . Dots denote time derivatives while primes denote spatial derivatives. These are the Josephson junction (JJ) equations for a two-component junction. Relative phases are coupled to currents Eq. (18) and to the relative motion Eq. (16). The two Josephson currents δṄ R 1 , δṄ R 2 are mutually coupled via the relative phases ∆φ i . This coupling signifies the Andreev-Bashkin (AB) effect [63][64][65][66].

VI. COHERENT DYNAMICS
In our simulations, the centers of droplets are separated by the distance |r R − r L | = r(0) 15. Due to their finite size, the distance between their surfaces is smaller, r(0) − 2R 0 ≈ 5 − 10. This is comparable to the range of the potential d = 1/ √ 2µ 1 1/ √ 2µ 2 ≈ 1.2 − 2.2. Initial droplets' velocities are zero ensuring low energy collision. Droplets' densities do not depart significantly from the equilibrium values. Any deviations result in a weak atomic evaporation. Fragmentation of droplets as observed in [19] occurs only at higher collision energy. In this regime more sophisticated approaches, like multiconfigurational time-dependent Hartree method [53], might be beneficial.
We begin our discussion by considering the case where both droplets have the same number of atoms and the phases of both components of a droplet are identical, ∆φ ≡ ∆φ 1 = ∆φ 2 . In such an arrangement the system can be described by a single wavefunction.
Our interpretation is based on the JJ equations. If the droplets are identical and initially |∆φ 1 ± ∆φ 2 | < π, they attract mutually and begin to move towards each other. The DC Josepshon current of atoms δṄ i ∼ sin(∆φ) flows from one droplet to the other in direction of the phase gradient. The droplets may merge forming an excited droplet. If the initial phases are such that interaction is repulsive, the droplets repel and move away from each other. The JJ current decreases in time as the coupling between droplets gets weaker. This simplistic description already shows that the motion of droplets may strongly depend on the relative phases of their wavefunctions.
However, the scenario given above is still oversimplified. Attraction between droplets can change to a repulsion. If the direct (DC) Josephson current triggered by a phase difference is large, one droplet may grow at the expense of the other and a difference between their chemical potentials will develop. This condition supports an alternating Josephson current (AC) rather then a DC one. The character of JJ dynamics switches to AC mode when the phase difference grows linearly in time, ∆φ ∼ ωt. The slow motion of droplets and fast phase dynamics happen on two different time scales. Averaging over fast oscillations according to P.L. Kapitza's method [67] allows to separate a fast micromotion from a slow relative movement governed by a repulsive ponderomotive potential This is a mechanism which leads an initially attractive potential eventually to act repulsively and allows the droplets to escape. Due to lack of confinement the distance between the droplets may change, thus the coupling strength varies in time. The Josephson current is a transient effect, and occurs only when the droplets are close to each other. This is a substantial difference from the case of two trapped BECs [47].
The Kapitza mechanism is illustrated in the top panel of Fig. 1, showing the interaction of two droplets with no initial phase difference, but with different atom numbers (and therefore different chemical potentials), N Another situation in which an initial DC current turns into an AC one is illustrated in Fig.1, middle panel, where we show a collision of two large, initially attractive droplets, N R 1 = N L 1 = 203(×2097) and ∆φ = π/4. As the droplets approach each other, the "down" one grows as it is continuously supplied by the DC current, while the "up" one becomes significantly smaller, until the AC mode takes over and the Kapitza mechanism comes into play. The smaller droplet is repelled, but before escaping it is recaptured due to the large amplitude of quadruple oscillations of the excited left droplet. This is an example of a very spectacular scenario where one of the droplets steals atoms from the other "at a distance" and finally devours its smaller companion.
When one of the phase gradients has opposite sign, as in Fig.1, bottom panel (∆φ 1 = −∆φ 2 = π/4), the droplets attract each other, and two equal but opposite Josephson currents are initiated. Such a scenario cannot be described by the single-component eGP equation.
The two counter-flows oscillate around zero while keeping opposite directions. These oscillations have a tendency to separate both species and accumulate them in opposite droplets. Excess atoms gather mostly in low density regions where some deviation from the equilibrium proportion does not destabilize the droplet [41]. None of the droplets is depleted and eventually they merge, similarly as if they had zero initial phase. For comparison see Fig.2, top panel, where merging of two small droplets, N R 1 = N L 1 = 20.3(×2097) is illustrated. Here however, no Josephson current nor phase dynamics take place. The resulting droplet is not as excited as the one discussed previously.
In the middle panel we show the Andreev-Bashkin effect, i.e. entertainment between the two superfluids. We simulate two small droplets N R 1 = N L 1 = 20.3(×2097) Other parameters, dimension of axes and colour coding is the same as in Fig.1. Video links: [71] (top), [72] (middle), [73] (bottom) with vanishing initial phase gradient ∆φ 1 = 0 in one component and ∆φ 2 = π/6 in the second one. If the currents were independent, only one species of atoms would flow from one droplet to the other. Yet both species begin to travel together and the related supercurrents are of the DC type. Eventually the "up" droplet loses so many atoms that it becomes unstable and evaporates during the approach. The minimal number of atoms supporting a stable droplet is about 18.65 × N 0 . The right droplet disappeared altogether while approaching its partner because both currents were acting in the same direction. The JJ equations do not account for evaporation and thus are not accurate at the final stage of such dynamics. The inset (middle panel, right) shows the entertainment between the two supercurrents δṄ R 1 − δṄ L 1 -black line, and δṄ R 2 − δṄ L 2 -red line. Note that the "black" current starts from zero but soon follows the "red" current. We stress that such a spectacular disappearance of one droplet in the presence of another is not unique and does not necessarily signify the AB effect. Very similar behaviour (not illustrated here) may be observed if both initial phase gradients are equal (e.g. ∆φ 1 = ∆φ 2 = π/10), i.e. if the two supercurrents act together from the very beginning. Therefore experimental proof of the AB effect requires precise control over the droplets' phases. In turn, if both initial phase gradients are equal to ∆φ = π/6, the JJ currents are larger than in the aforementioned case (Fig. 2, bottom panel), droplets initially approach but the AC sets on while the droplets are still far away from each other. Thus repulsion takes over before the smaller droplet disappears.
Experimental realizations of phase dependent collisions crucially depend on the sensitivity of the entire process with respect to the stability of the relative phase difference. To investigate this problem we randomly perturb the initial phases of the droplets. After preparing the initial state with a well defined phase difference, we add to both components some local phase perturbations: where δφ i,L(R) (r) is a random variable uniformly distributed in the interval [−δ, δ]. In the example we show below, we chose the initial phase difference ∆φ 0 = 3π/4 and δ = 0.1×∆φ 0 . The perturbation is thus quite strong, at the level of ±10%. The kinetic energy of the system after this perturbation is so large that the total energy of the droplets becomes positive, signifying thermal instability. Not only soft, but also hard modes are excited this way, and the phase difference between the two components varies from point to point. These phase fluctuations convert into density fluctuations at early stages of the dynamics, on a time scale of about 0.05-0.1 ms. In Fig.(3) relative density fluctuations of the right droplet (normalized to the non-perturbed density) are plotted. Fluctuations of the order of ±2% are clearly visible. As densities depart locally from their equilibrium values, particles are emitted, which ultimately ceases both density and phase fluctuation. The phases stabilise at some values not significantly different from the initial setting (see right panel of Fig(3)). The internal dynamics of the droplets settle down, and the droplets continue their mutual approach until they eventually escape in opposite directions, see Fig(4). The right panel shows dynamics of unperturbed droplets, and in the left panel a collision of perturbed droplets is illustrated. Both processes look very similar, however careful inspection suggests that the repulsion of perturbed droplets starts marginally later than in the analogous undisturbed setting. tions of phases of droplets is a very interesting feature which deserves detailed study in the future.

VII. CONCLUSIONS
In conclusion, we showed that the dynamics of interacting droplets and their ultimate fate depend crucially on the relative phases of their wavefunctions. Two liquid quantum droplets, which constitute identical macroscopic objects, can be made to merge, repel or evaporate only by manipulating their quantum phases. Thus the processes studied in this work are macroscopic manifestations of the quantum nature of ultracold droplets. The interaction potential derived here as well as the twocomponent Josephson-junction equations may prove useful in studying the Andreev-Bashkin effect or modelling arrays of coupled droplets in a supersolid-like arrangement. Our Josephson-junction equations do not have any free parameters -the chemical potentials and their derivatives as well as the long-range behaviour of the droplets' wavefunctions were found from stationary solutions of the extended Gross-Pitaevskii equations. The collisions discussed in this work enfold over a time span of more than 100 ms. The experimental verification of our predictions depends crucially on the realization of long-lived quantum droplets. One of the possible ways towards experiments are heteronuclear droplets in lower dimensions. On the other hand, one component Josephson junction dynamics is also possible in the case of dipolar droplets with sufficient lifetime. In fact Josephson-like oscillation of phases and atom number in droplets forming a supersolid were reported [21]. Collisions of dipolar droplets remain to be studied.
We believe that experimental advances in creating long-living droplets will sooner or later allow for verification of our predictions in their full extent.