Dynamics of quantum droplets in an external harmonic confinement

Recent theoretical and experimental results show that one-dimensional (1D) weakly interacting atomic Bose-Bose mixtures with repulsive interspecies mean field (MF) interaction are stabilized by attractive quadratic beyond-mean-field (BMF) effects into self-bound quantum droplet (QD) in free space. Here, we construct an exact analytical model to investigate the structure and dynamics of QDs in presence of external harmonic confinement by solving the 1D extended Gross–Pitäevskii equation (eGPE) with temporal variation of MF and BMF interactions. The model provides the analytical form of wavefunction, phase, MF and BMF nonlinearities. The generation of QDs and interesting droplet to soliton transition in presence of regular/expulsive parabolic traps by taking the comparable MF and BMF interactions are illustrated. We derive the phase diagram of the droplet-soliton phase transition between amplitude of MF, BMF interactions and harmonic oscillator frequency. The strength and form of oscillator frequency are identified as key parameter for tuning the compression, fragmentation and transport of droplets. Finally, the stability of the obtained solutions are confirmed from Vakhitov–Kolokolov (VK) criterion and are found stable.

The formation of droplets is one of the characteristic feature of classical liquids which results due to counter balancing of attractive van der Waal's interactions with repulsive interactions generated from high density of its constituents 1 . Recently, in one of the pioneering experiment in the field of Bose-Einstein condensate (BEC) and superfluids, a new type of dilute quantum droplets (QDs) are observed in Bose-Bose mixtures with isotropic contact interactions 2,3 and anisotropic dipolar interactions 4,5 . Different from the classical liquids, these QDs are generated at very low temperatures ( 10 9 order less than classical liquids) in BEC phase which makes them ultradilute quantum liquid with its density 10 8 orders of magnitude lower than liquid helium's 6 . Although, the formation of liquid state in weakly interacting BEC is not possible according to mean-field theory, however, Petrov pointed out that the QDs phase can be stabilized when the dominant contribution to the energy due to the contact mean field (MF) interaction is close to vanishing and becomes comparable to the beyond mean field (BMF) effects Lee-Huang-Yang (LHY) corrections arising from quantum fluctuations 7 . Based on this stabilization proposal, the dilute QDs are experimentally realized in binary Bose-Bose mixtures 8 , dipolar gases 9 and similar stabilization mechanisms are proposed for its existence in Bose-Fermi mixtures 10,11 and dipolar mixtures 12,13 . These developments led to a plethora of works in this field [14][15][16][17][18][19][20][21] . It is relevant to mention that the formation of vortices in droplets 22,23 , multiple droplets 24 , existence of striped states 25 , droplet clusters 26 and supersolid properties in phase coherent droplets 27 , are also predicted in BECs. The precise experimental control and tunability over the QDs parameters: population density of mixture constituents, inter-and intra-component interactions, confinement, temperature and spatial dimension, opens new direction in modelling the other complex manybody systems like liquid helium 6,21 . Further, one possible application of the QDs is considered in the precise designing of matter-wave interferometer 28 .
In this work, we aim to construct an exact analytical model for investigating the generation of QDs in cigar-shaped two component BEC mixture and study its spatio-temporal dynamics in the presence of external harmonic confinement. In order to construct the analytical model, we consider binary Bose-Bose mixture in one dimensional geometry. The choice of binary BECs offers the presence of competing intra-and inter-component interactions which is required for the generation and stabilization of QDs 7 and these interactions are experimentally tunable using Feshbach resonance technique 29 . Additionally, the quasi 1D condensates are quite useful for studying fundamental physics problems and a variety of nonlinear excitations like solitons (bright, dark, grey), compactons, rogue waves, Faraday waves, etc. are theoretically and experimentally explored in it 30 . Further, the choice of 1D geometry offers significant stability advantage of nonlinear excitations in comparison www.nature.com/scientificreports/ to three-dimension geometry and a clean controllable many-body test bed for precise measurement of quantum many-body effects. Although, in recent times, the dynamics of QDs is investigated in 1D geometry in absence of any external confinement i.e. free space. Astrakharchik et. al. identified the two physically different regimes (Gaussian shape and broad flat top plateau) of QDs depending its size in free space using eGPE 31  In the following section, we provide an exact analytical model for investigating the dynamics of QDs in a weakly interacting 1D mass-balanced binary Bose-Bose mixture in the presence of external harmonic trap using 1D eGPE. Here, we have considered the system with temporally varying competing repulsive cubic MF and attractive quadratic BMF interactions. The model and analytical framework for the calculated system variables is explained by finding the wavefunction, oscillator frequency, MF and BMF nonlinearities. It is shown that with suitable choice of the harmonic oscillator frequency, one can realize following experimentally realizable forms of trap configurations: (a) free space; (expulsive harmonic). Moreover, the the consistency conditions governing the oscillator frequency interestingly maps with the linear Schrödinger equation which leads to a significant temporal control over the droplet dynamics of the system. We then investigate the dynamics of condensate under the influence of regular and expulsive parabolic traps. In comparison to free space, the axial compression of the atomic number density is observed in above-mentioned confinement with the variation of oscillator frequency. We illustrate the droplet to soliton transition in the considered system and remarkably identify the strength of oscillator frequency and sign of MF, BMF interaction's amplitude as the key parameters for realizing this transition. For a better insight of the phenomenon of droplet-soliton transition, a phase diagram between the MF, BMF interaction amplitude's and the oscillator frequency is also plotted. Moreover, by choosing the elliptic solution of eGPE and oscillator frequency of trap, we achieve fragmentation, merging and transport of QDs across the potential. The stability of obtained solution is confirmed using Vakhitov-Kolokolov (VK) criterion through an estimation of slope of norm w.r.t. chemical potential.

Results
The model and analytical framework. We begin by considering the binary BEC with equal masses and equal number of atoms in the components under the influence of BMF (LHY corrections for quantum fluctuations) in presence of temporally varying external harmonic confinement. The same magnitude of equilibrium densities of binary mixture components makes the result analysis clearer and easier. In 1D geometry, the system is described by the following equation 14,33 : where ψ 1 ( ψ 2 ) are the wavefunctions of first and second component of binary mixture. Further, we assume that the coupling constants describing the repulsion between the atoms in each components are equal: g ↑↑ = g ↓↓ ≡ g = 2 2 a s (t)/(ma 2 ⊥ ) and g c = g ↑↓ . Here, � s (t) = (g c + 3g)/2 represents the self interaction coefficients whereas � c (t) = (g c − g)/2 is the cross interaction coefficients along with Ŵ(t) = √ mg 3/2 /(π ) 33 . We also consider binary mixture under influence of an external harmonic trapping potential V(x, t) with frequency ω T , i.e., mω 2 T x 2 /2 . a ⊥ is the transverse oscillator length of the trap and a s (t) represent the time dependent interand intra-components atomic scattering lengths which is tunable through Feshbach resonance technique 29 . www.nature.com/scientificreports/ Considering the binary mixture as mutually symmetric spinor components with ψ 1 = ψ 2 = ψ 0 ψ , then the dynamics of the considered system reduces to the time-dependent dimensionless single eGPE 14 : Here, ψ(x, t) can be treated as the condensate wave function of the droplet with mass m and |ψ(x, t)| 2 represents the normalized density of each component. The parameters g 1 (t) = Ŵ(t) and g 2 (t) = � s (t) + � c (t) are the binary mixture components coupling constants representing BMF and MF interactions. In experiments it is possible to tune g 1 (t) and g 2 (t) , both to positive or negative values by changing the sign of atomic scattering length. Here, g 1 (t) and g 2 (t) are non-zero functions. Further, M(t) is the time dependent harmonic oscillator frequency; a constant M(t) implies an oscillator potential which can be confining or expulsive for M > 0 or M < 0 , respectively. Here, the value of scaling parameters x 0 , t 0 , ψ 0 are 31 . Equation (3) is the extended GPE added with external harmonic confinement and for M(t) = 0 , it has been extensively used in the past theoretical studies to describe the structure and dynamics of quantum droplets 7,21 .
In order to obtain analytical solutions of (3) in presence of of the time-dependent harmonic trapping potential and discuss QDs spatio-temporal dynamics in it, we assume the following ansatz solution: where A(t) and ξ(x, t) = γ (t)x + δ are the time-modulated amplitude and travelling coordinate of the condensate, respectively. Here, γ (t) is the positive definite functions of time with δ > 0 . In principle, the γ (t) and −δ/γ (t) represents the inverse of the width and the position of the center of mass of the excitation 34 , respectively. Further, we assume the quadratic form of the phase as: Now, by substitution of the ansatz Eq. (4) and phase form (5) in Eq. (3), we obtain the following consistency conditions: where E is eigenvalue of Eq. (6), G 1 , G 2 denote nonlinearity constant which can have positive or negative depending on scattering length along with: where M(t) is the oscillator frequency of the harmonic confinement. It is important to note here that if γ (t) is constant then M(t) = 0 i.e. no external harmonic trap. The Eq. (9) can be transformed into the well-known linear form of Schrödinger equation by taking γ (t) = 1 a(t) : Equation (10) reveals a non-trivial correlation between the oscillator frequency and the amplitude of condensate wavefunction. This is one of the main results of this paper. By taking advantage of this connection, one can in principle identify QDs corresponding to each solvable quantum-mechanical system. Further, the Eq. (10) can be transformed into the Riccati equation, The merit of these correlations is that it provides freedom to control the dynamics of the QDs in a number of analytically tractable ways as the Schrödinger equation and the Riccati equation can be exactly solved for a variety of M(t). The solution of Eq. (6) can be given as 7,31 : where, µ 0 = −2/9 with E < 0 , G 1 < 0 and G 2 > 0 . Thus, the complete solution of Eq. (3) can be written as: www.nature.com/scientificreports/ Equation (13) represents the localized solution of Eq. (3). Below, we study the dynamics of QD's by considering the exact solution (13) of eGPE in presence of temporally varying repulsive cubic MF and attractive quadratic BMF interactions in harmonic confinement. Based on Eq. (9), we investigate the droplet dynamics for the following experimentally relevant forms of harmonic traps with suitable choice of oscillator frequencies M(t) and γ (t) : where, γ 0 , and M are real positive constants. In principle, depending on the need of the physical situation, one can choose the suitable γ (t) form for controlling the depth and width of the resultant traps. Here, in the following sections, we demonstrate the generation of QD's in above mentioned trap configurations i.e. in absence/presence of regular/expulsive parabolic traps.
Quantum droplets in free space with M(t) = 0. In this section, we investigate the generation of QD in absence of external confinement i.e. free space. In the 1D geometry, the dynamics of QDs is extensively investigated in free space scenario 7,21,31 . Here, from the constructed analytical model, we write down the wavefunction solution of Eq. (3) for V (x, t) = 0 . For this, we consider γ (t) =constant = γ , which ensures that M(t) becomes zero from Eq. (9). Thus, the complete wavefunction solution form of Eq. (13) for δ = 0 becomes: where γ > 0 , with E < 0 , G 1 < 0 and G 2 > 0 . The obtained form (14) of free space wavefunction solution for eGPE (4) is identical to the results reported in the literature with G 1 = −1 7,21,31 . Further, for this form of γ (t) , the BMF and MF interactions becomes: g 1 (t) = (1/2)G 1 γ 3 2 , and g 2 (t) = (1/2)G 2 γ , respectively. Depending upon the sign of G 2 , Eq. (14) leads to investigation of QDs dynamics in three characteristic regimes: (a) G 2 > 0 : flat top density profile i.e. droplet domain, (b) G 2 < 0 : integrable GPE with bright soliton solution with quantum fluctution term approaching to zero, and (c) G 2 is very small: GPE with only quadratic nonlinearity giving rise to the droplet wave function of the Korteweg-de Vries-type 32 .
In Fig. 1a, we plot the variation of condensate density profiles for the wavefunction (14) w.r.t. the changing magnitude of mean field interaction i.e. G 2 with G 2 > 0 . Here, the physical parameter values: γ (t) = 1 , E = −2/9 , G 1 = −1 , along with G 2 = 0.99 (dotted purple), 0.9999 (dotdashed blue), 0.999999 (dashed red) and 0.9999999 (thick black line). It is apparent from the figure that as MF interaction strength G 2 → 1 and becomes comparable to BMF interaction strength G 1 = −1 , the signature feature of QD's i.e. flat top density profile is realized. The variation of condensate density with the MF interaction parameter G 2 is in conformity with the physical situations reported in the literature 32 .
In Fig. 1b, we compare the profile of condensate densities for the regular harmonic trap (dotdashed blue line) with the free space domain (black line) using wavefunction solutions (14) and (15)   www.nature.com/scientificreports/ condensate density having QD profile which with increase in time demonstrate strong compression at t = nπ/2 , where n is a real number, and these condensate atoms coalesces forming large density peaks. Here, other physical parameters are: γ 0 = 1.5 , E = −2/9 , G 1 = −1 , and G 2 = 0.999999 . This behavior of condensate density can be understood by plotting the g 1 (t) (thick line) and g 2 (t) (dashed line) profiles w.r.t. M at t = 1 (Fig. 2c) for the same parameter values. It shows that the magnitude of g 2 (t) and g 1 (t) increases sharply near nπ/2 points due to presence of sec(Mt) term which physically represents a strong repulsive MF and attractive BMF interactions. Further, the magnitude of g 1 (t) is proportional to [γ 0 sec(Mt)] 3 2 leading to compression of condensate atomic density. Analytically, the sign of g 1 (t) and g 2 (t) changes at nπ/2 values i.e. in the [0, π/2] and [3π/2, 2π] domains, g 2 (t) is positive whereas in between [π/2, 3π/2] , this becomes negative whereas g 1 (t) is negative in [0, π/2] and [3π/2, 2π] . The change of sign g 2 (t) results in making MF interactions repulsive to attractive with variation of M. As the balance of MF and BMF interaction is required for realizing QD and as discussed earlier, for g 2 (t) > 0 , the density profile of Eq. (3) will be droplet whereas if g 2 (t) < 0 , then it becomes soliton 7,32,35 . In the time interval [π/2, 3π/2] , the considered system will have attractive MF interaction with continuous loss of condensate atoms due to imaginary BMF term leading to formation of bright solitons 36 . On this basis, we reveal an analytical scenario for the droplet to soliton in regular harmonic trap with the change in sign of ±g 2 (t) in Fig. 2c at nπ/2 points. Based on this, in Fig. 2d, we construct a phase diagram in between γ 0 i.e. amplitude of MF, BMF interaction's and M. For γ 0 > 0 , g 2 (t) is positive in [0, π/2] and [3π/2, 2π] domains i.e. droplet density profile whereas it becomes negative in between [0, π/2] with γ 0 < 0 resulting into soliton formation. This is one of the important result of this article revealing non-trivial correlation between the droplet to soliton phase transition with M and γ 0 . Here, all the parameters are in dimensionless unit with magnitude: γ 0 = 1.5 , E = −2/9 , G 1 = −1 , and G 2 = 0.999999 . Further, in Fig. 2b, we illustrate the impact of the strength of M on single QD profile using (15). It is apparent from the Fig. 2b that with M tending from 0 → 1 leads to increase in height and decrease in the width of QDs profile due to its compression. In principle, by tuning the oscillator frequency of external confinement, one can tune the width and height of QDs. Here, we depicted the profile with physical parameters: Fragmentation and merging of QDs by tuning oscillator frequency: We demonstrate a novel mechanism for the fragmentation and merging of QDs with the tuning of oscillator frequency of regular harmonic trap. We show that by choosing the elliptic solution of Eq. (6) in place of localized solution, one can induce the fragmentation in QDs. In order to demonstrate this effect, we write the solution of Eq. (6) as: Here, q is the modulus parameter for the Jacobi elliptic function cn. For investigating the dynamics of condensate in the regular harmonic confinement, we take the same parameter values of Fig. 2 i.e. γ 0 = 1.5 , E = −2/9 , with γ (t) = γ 0 sec(Mt) . Thus, using this, the complete wavefunction solution of (13) becomes: with δ = 0 , M > 0 with E < 0 , G 1 < 0 and G 2 > 0. Figure 3a and b depicts the impact of magnitude of M on the fragmentation and merging of QDs for M = 0.75 and M = 0.5 with q = 0.9 . Here, we have taken comparable strength of MF and BMF interaction strengths with G 1 = −1 , and G 2 = 0.999999 , respectively so that density profile remains in droplet regime (Fig. 2b). For exploring the condensate dynamics, we considered potential in the range [−1, 1] . It is evident from the figures that, initially at t = 0 , a single QD profile is observed and with increase in time, it gets fragmented and forms the QD train. Further, with increase in time, these fragmented components starts to coalesce and again merges to form a single QD. Importantly, with M changing from 0.5 → 0.75 , the number of fragmented droplets increase. In Fig. 3c,  Transport of droplets across potential: In Fig. 4a, we illustrate an analytical scenario for the oscillation and transport of QDs in regular harmonic confinement in the considered system. In order to realize the transport of QDs in parabolic trap, we have taken δ = 0 in Eq. (13) along with γ (t) = γ 0 sec(Mt) which ensures the form of trap as (1/2)M 2 x 2 . Thus, the resultant wavefunction becomes: with γ 0 > 0 , M > 0 , E < 0 , G 1 < 0 and G 2 > 0.
In Fig. 4a, we plot the condensate density profile for wavefunction solution (18)    www.nature.com/scientificreports/  www.nature.com/scientificreports/ with [ −δ/γ ((t) ] and is oscillating w.r.t. time 37 . For the considered case, −δ/γ ((t) = −(δ/γ 0 ) cos(Mt) i.e. for M = 0.5 , the maximum magnitude of −(δ/γ 0 ) cos(Mt) appears at t = 2nπ . This is validated from Fig. 4a as the direction of oscillation of QDs changes at t = 2nπ locations. For providing a better interpretation of the oscillation of QDs in this system, we plot −(δ/γ 0 ) cos(Mt) w.r.t. M in Fig. 4b  , and g 2 (t) = (1/2)G 2 γ 0 sech(Mt) , respectively. It is apparent from Fig. 5c that g 2 (t) and g 1 (t) have maximum magnitude at M = 0 and t = 1 i.e. strong repulsive and attractive interactions whereas their amplitude decreases with increase in M for γ 0 > 0 . However, g 1 (t) is proportional to [γ 0 sech(Mt)] 3 2 leading to strong compression. Similar, temporal profile of g 1 (t) , g 2 (t) can be obtained by keeping M = 1 and varying time. Fig. 5c clearly justifies the strong compression in condensate density at t = 0 in Fig. 5a. It is important to note here that the profile of Fig. 5c get reversed for γ 0 < 0 which means g 2 (t) becoming attractive with continuous loss of condensate atoms due to imaginary BMF term leading to formation of bright solitons 36 . As discussed above, in this domain, the condensate density becomes soliton. Based on this, we construct a phase diagram in between γ 0 and M in Fig. 5d which shows that droplet to soliton transition in expulsive trap can be observed by changing the sign of MF and BMF interaction's amplitude. Different from regular harmonic confinement, the droplet to soliton transition occurs in expulsive trap by tuning the sign of MF and BMF interactions. It is due to the change in nature of MF, BMF interactions from sec(Mt) (regular harmonic) to sech(Mt) (expulsive harmonic). Figure 5b illustrates the variation of droplet profile w.r.t. changing magnitude of M. With increase in magnitude of M from 0 → 1 leads to decrease in height and increase in the width of QDs profile due to its expansion with changing strength of oscillator frequency. In principle, by tuning the oscillator frequency of external confinement, one can tune the width and height of QDs. Here, we depicted the profile with physical parameters: γ 0 = 1.5 , E = −2/9 , G 1 = −1 , G 2 = 0.999999 and t = 1.

Stability analysis.
Using various forms of γ (t) in Eq. (13), we investigated the generation and spatio-temporal dynamics of droplets in absence/presence of regular/expulsive traps. In this section, we perform the stability analysis of (13) using the VK criterion 41 . The VK criterion is widely used for the stability analysis of nonlinear Schrödinger equation and it is already established in the literature that necessary condition for a solution to be stable is positive slope of dN/dE, where N is the normalization and E is the chemical potential of the system. If the dN/dE of obtained solution is negative then the solution is unstable and for the condition that dN/dE = 0 , the solution remains marginally stable. Now, using Eq. (13) and N = +∞ −∞ |ψ| 2 ∂x , one can estimate the correlation between normalization N and E as: where G 1 = −1 and G 2 G 2 1 ≈ 1 . Equation (20) estimates the magnitude of N in presence of harmonic trap and is equal to the N reported for free space 31,32 . Thus, even in the presence of external trap, the system has a continuous symmetry property and N is conserved in time as per Noether's theorem 42 . In Fig. 6, using Eq. (20), we plot dN/dE w.r.t. E where G 1 = −1 . It is evident from the figure, that the magnitude of dN/dE is positive w.r.t. its variation E which indicates the stable nature of the obtained solution. www.nature.com/scientificreports/  www.nature.com/scientificreports/

Summary
The main results reported in this paper are summarized as follows. We report a large family of exact analytical solution of 1D eGPE for investigating the structure and dynamics of QDs in harmonic confinement for temporally varying MF and BMF interactions. We have considered a weakly interacting 1D mass-balanced binary Bose-Bose mixture with competing repulsive cubic MF and attractive quadratic BMF interactions. The dynamics of QDs investigated w.r.t. the oscillator frequency and the strength of BMF interaction in different experimentally realizable forms of trap configurations: (a) V (x, t) = 1 2 M 2 x 2 ; (regular harmonic); and (b) V (x, t) = − 1 2 M 2 x 2 ; (expulsive harmonic). The constructed solutions for the free space is in agreement with existing results in the current literature 31 . Compression and flat top profile of the atomic number density is illustrated in above-mentioned confinement with the variation of M and G 2 . Such flat top density profile was reported in the literature for QDs in free space domain. We have observed the droplet to soliton transition in regular/expulsive harmonic confinement depending on the magnitude of M and sign of γ 0 . We have produced the phase diagram between the BMF interaction amplitude and the oscillator frequency for a better insight of droplet-soliton transition. Utilizing the elliptic solution of eGPE and oscillator frequency, we demonstrates a novel mechanism for fragmentation, merging and transport of QDs. The stability of obtained wavefunction solutions are confirmed using VK criterion.
As an extension of the present work, it may be interesting to verify the collision of small/large droplets, its impact on fragmentation/merger and study its dynamics in 2D/3D dimensions in presence of harmonic trap.

Methods
We consider the Bose-Bose mixture BEC with equal masses and equal number of atoms in the components under the influence of BMF (LHY corrections for quantum fluctuations) in presence of temporally varying external harmonic confinement. By taking the spinor components of binary mixture as mutually symmetric, the Hamiltonian is taken as i.e. in the form of time-dependent dimensionless 1D eGPE. Using similarity transformation technique, we assume an ansatz solution and connect the 1D eGPE with a solvable PDE. This results into a number of constraints over phase, MF and BMF nonlinearities as PDEs. Next, we solve these equations and obtain the general form of wavefunction solution.

Data availability
All data generated or analysed during this study are included in this published article. It can be reproduced by utilizing the form of wavefunction and considered trap form.