Dynamical formation of quantum droplets in a $^{39}$K mixture

We report on the dynamical formation of self-bound quantum droplets in attractive mixtures of $^{39}$K atoms. Considering the experimental observations of Semeghini et al., Phys. Rev. Lett. 120, 235301 (2018), we perform numerical simulations to understand the relevant processes involved in the formation of a metastable droplet from an out-of-equilibrium mixture. We first analyze the so-called self-evaporation mechanism, where the droplet dissipates energy by releasing atoms, and then we consider the effects of losses due to three-body recombinations and to the balancing of populations in the mixture. We discuss the importance of these three mechanisms in the observed droplet dynamics and their implications for future experiments.


I. INTRODUCTION
Self-bound droplets of ultracold atoms have been recently discovered as a new exotic quantum phase [1][2][3][4]. Although forming in dilute atomic gases, they display properties which are unique in that context but common to different systems such as classical liquids, helium nanodroplets or atomic nuclei. After their existence in attractive bosonic mixtures was theoretically predicted in [1], they have been experimentally observed in dipolar condensates [5][6][7][8], homonuclear mixtures of 39 K [9,10] and recently in a heteronuclear mixture of 41 K and 87 Rb [11].
The first pioneering experiments performed with homonuclear mixtures at ICFO [9,12] and LENS [10,13] have been able to demonstrate the existence of self-bound droplets in these systems and to provide a first characterization of their peculiar features. While these works have been mainly devoted to study the droplets equilibrium properties, the experiment reported in Ref. [10] has also shown that during the formation of the droplet an interesting and complex dynamics takes place. The non-adiabatic preparation of the mixture leads to an initial compression of the atomic cloud and to following oscillations of its size, while the presence of strong threebody losses (3BL) continuously drives the system out of equilibrium. During this evolution, the droplet eventually reaches a metastable state, which was the main focus of the investigation performed in [10]. In this work we concentrate instead on the dynamical evolution observed during the droplet formation, trying to understand the different mechanisms involved and their relative importance. While three-body recombinations are a well-known phenomenon in the field of ultracold atoms, the other two mechanisms playing a role in this evolution are specific of quantum droplets and thus require further attention. The first is related with the need to adjust the populations in the components of the mixture to bal-ance its interaction energy. The second, the so-called self-evaporation, is a more peculiar dissipation mechanism predicted in Ref. [1]. Calculating the excitation spectrum of the droplet, Petrov noticed that, in some specific conditions, the droplet cannot host any discrete excitation, since all the excited states are higher in energy than the particle emission threshold. This suggested the idea that the droplet could be able to dissipate any excess of energy by expelling atoms, from which the term self-evaporation.
In this paper, we consider the specific experimental case of Ref. [10] and use numerical simulations to understand how these different mechanisms come into play in the evolution of the droplet. In Sec. II we summarize the conditions for the existence of self-bound states in a 39 K atomic cloud at zero temperature. In Sec. III we analyze the phenomenon of self-evaporation, considering the ideal case of a mixture without 3BL. We first study the linear regime, where the droplet is prepared with a small initial excitation, compatible with the assumptions of Ref. [1], and then we investigate how the concept of self-evaporation extends to the more realistic case where the mixture is prepared far from equilibrium. In Sec. IV we introduce 3BL, thus fully recovering the experimental conditions of Ref. [10], and we analyze the dynamical evolution of the droplet, identifying which mechanisms play a key role in the droplet dynamics. Finally, conclusions are drawn in Sec. V.

II. SELF-BOUND DROPLETS IN 39 K
The key ingredient for the formation of self-bound atomic clouds is the competition between attractive and repulsive forces which, scaling differently with the atomic density, generate a binding potential [1]. In a mixture of 39 K atoms in the hyperfine states |1, 0 ≡ |1 and |1, −1 ≡ |2 , this situation occurs in a specific range of magnetic fields B, where the intraspecies scattering lengths a 11 and a 22 are positive, while the interspecies a 12 is negative. For B < 56.85 G, the quantity δa ≡ −|a 12 | + √ a 11 a 22 becomes negative, so that the global mean-field (MF) interaction is attractive. In this regime, a repulsive effect is provided by quantum fluctuations, corresponding to the so-called Lee-Huang-Yang (LHY) energy term [14], that stabilizes the system against collapse and gives rise to a self-bound atomic droplet. As derived in Ref. [1], the competition between the LHY energy and the attractive MF term locks the ratio between the equilibrium densities in the two species to with The experiments reported in Refs. [9,10] have verified the existence of these self-bound atomic clouds in the predicted interaction regime at the nominal population ratio of Eq. (1).

III. SELF-EVAPORATION
As discussed in Ref. [1], for small atom numbers (see later on for a more precise definition) the droplet has no collective modes with energy lower than the particle emission threshold. Therefore, no bound collective excitation can be sustained in that regime, so that any perturbation of the equilibrium state would result in a release of particles or in a break-up of the droplet into smaller pieces. In this sense, a quantum droplet represents a self-evaporating object. In this section we use numerical simulations to characterize this phenomenon, looking at the dynamical evolution of an excited droplet. We first consider the linear regime, where the system is prepared with a small initial excitation, and then we extend our analysis to a regime closer to the experimental conditions, where the droplet is prepared in a highly excited state.

A. Linear regime
For small-amplitude excitations, we can describe the system with a single wave function φ(r, t), neglecting any possible relative motion of the mixture components. We use the same formalism of [1], which we briefly summarize here. We introduce the rescaled spatial coordinate ρ = r/ξ, with and the rescaled time τ ≡ t/mξ 2 . The droplet wave function φ(ρ, τ ) evolves according to the time-dependent Gross-Pitaevskii equation where |φ| 2 d 3 ρ = N defines the chemical potentialμ and N is related to the number of atoms in the two atomic species N i by It was predicted in Ref. [1] and confirmed by the experiments in Refs. [9,10] that stable droplets exist only for N > N c ≈ 18.65. The regime of self-evaporation corresponds to 20.1 < N < 94.2 [1].
In the following, we will restrict our analysis to the case of a spherically symmetric system, which reduces Eq. (4) to an effectively one-dimensional equation which depends on the radial coordinate ρ only. This assumption will be maintained throughout this work, for easiness of calculations and conceptual clarity. Note that in this approximation the only possible excitation is that with angular momentum ℓ = 0, i.e., the monopole mode. Higher angular momentum modes are not accounted for, so that the system is self-evaporating up to N ≃ 934, where the monopole mode reenters into the spectrum (see Ref. [1]).
To study the dynamics of an excited droplet, we prepare the system slightly out of equilibrium, by solving the stationary version of Eq. (4) with a finite tolerance [15], so that the initial wave function φ (i) 0 (ρ) corresponds to an excited state with a slightly larger energy than the actual ground state φ 0 (ρ) [16]. We then calculate the evolution of the system by solving Eq. (4) with the above initial condition. To distinguish the droplet from unbound expanding components that may form during the evolution, we define a droplet volume as that contained within a certain bulk radius R d (τ ), defined as the position of the minimum ofn(ρ, τ ) ≡ ρ 2 n(ρ, τ ) (see Appendix A for details) [17].
The evolution of the droplet size σ(τ ), defined as the rms size of the density distribution within the droplet volume, is shown in Fig. 1. We find two different behaviors depending on the value of N . For low values of N [in panel (a)], we observe a damped oscillation with decreasing frequency ω 0 (τ ). In this regime, where all the droplet excitations lie in the continuum, the droplet cannot host a bound monopole mode. This means that when the cloud starts oscillating, part of the atoms move away from it, thus reducing the droplet energy by ∆E = −µ∆N + E kin , where ∆N is the number of atoms leaving the droplet and E kin is their average kinetic energy. For large N instead, when the monopole has reentered into the spectrum [panel (b)], the size σ performs a sinusoidal oscillation with a well defined frequency ω 0 . To characterize the latter regime, we fit σ(τ ) for N > 940 with A cos(ω 0 τ ) to extract the monopole frequency. For lower N instead, we fit the oscillation with separate damped cosine functions of the form A cos(ω 0 τ ) exp (−τ /τ ) on different intervals of length ∆τ = 25, and we extract the corresponding frequency ω 0 and damping rateτ . In Fig. 2a we report the fitted values of ω 0 for different values of N and we compare them with the theoretical predictions of Ref. [1] for ω 0 ( N ) and for the emission threshold −μ. For N > 940 we find a very good agreement between our numerical results and the prediction for ω 0 . In the self-evaporation regime we see that the oscillation frequency decreases in time until it reaches its lower value, set by −μ. In Fig. 2b we plot the extracted values of ω 0 (τ ) andτ (τ ) from the fits of Fig. 1a. In Fig. 2(c,d) we report the variation of the parameter N and of the droplet energy E in the corresponding time We observe that the cloud initially oscillates with large frequencies and fast damping rates, associated with a significant release of atoms from the droplet into an unbound cloud. As the droplet energy decreases due to atom losses, the oscillation slows down and ω 0 eventually saturates at −μ. Close to this threshold the atoms that leave the droplet carry away ≃ −μ and thus have negligible kinetic energy. The damping of these small final oscillations is then extremely slow and the droplet reaches the stationary ground state only asymptotically. However, these residual excitations are extremely small and one can effectively consider that the droplet has reached its equilibrium configuration well before the asymptotic regime is reached.

B. Non-linear regime
In the experiment of Ref. [10], the mixture is prepared far from the equilibrium profile φ 0 of a droplet. Then, to understand the dynamics which takes place during the droplet formation we extend the previous analysis beyond the linear case. Since we cannot exclude a priori a different behavior of the two components of the mixture, here we replace Eq. (4) with a system of two coupled generalized Gross-Pitaevskii (GP) equations [18] including LHY corrections [10,11] (we keep the assumption of spherical symmetry to simplify the discussion) where ψ i (r) are the wave functions of the two species and the corresponding chemical potentials are given by with g ij = 4π 2 a ij /m. The LHY energy is [1] so that Note that a similar model, obtained by means of a local density approximation, has been used also for the description of dipolar quantum droplets [19][20][21].
In our numerical simulations we consider a preparation of the mixture similar to the one implemented in the experiment of Ref. [10]. A single species Bose-Einstein condensate (BEC) of N atoms of 39 K in the hyperfine level |2 is prepared in the ground state of a harmonic trap with trapping frequency ω. At t = 0, N 1 atoms are instantaneously transferred to |1 , with N 2 = N − N 1 atoms remaining in the initial state. The harmonic potential is then turned off to study the evolution of the mixture in free space. Here, we consider a typical set of experimental parameters, with N = N 1 + N 2 = 4 × 10 5 and scattering lengths a 11 ≃ 69.99a 0 , a 12 ≃ −53.37a 0 and a 22 ≃ 34.11a 0 , which correspond to a Feshbach magnetic field B = 56.45 G. For this set of parameters, we have N = 200, which lies in the self-evaporation regime identified above.
The ground state wave function ψ 0 of the initial condensate is obtained from the following stationary GP equation [18] with |ψ 0 | 2 d 3 r ≡ 1. Considering the preparation sequence described above, we simulate the instantaneous transfer of atoms in |1 by assuming that the (normalized) wave functions of the two components at t = 0 are equal to ψ 0 , with the corresponding densities being n i = N i |ψ 0 | 2 . Here we keep the ratio N 1 /N 2 fixed to the nominal ratio in Eq. (1). The following evolution is then obtained by solving Eq. (6). One can easily guess that the droplet is minimally excited when the trapping frequency ω is such that the initial atomic distribution is as close as possible to that of the droplet ground state for that value of N , namely |ψ 0 (r)| 2 ≈ ξ −3 |φ 0 (r/ξ)| 2 . To quantify the difference between the initial profile of the condensate and that of the droplet, one can consider, e.g., the relative deviation between the corresponding energies, ∆[ N , ω] ≡ |E(ψ 0 ) − E(φ 0 )|/E(φ 0 ). For the present case, we find that ∆[ N = 200, ω] is minimized for a trapping frequency ω/2π ≃ 600 Hz.
We then study the dynamics of the system as a function of the distance of the initial state with respect to φ 0 , by varying the trap frequency ω. Similarly to the case discussed in the previous section we observe that, after the mixture is formed and released into free space, the atomic clouds undergo damped oscillations. For ω/2π = 600 Hz the density of the binary mixture smoothly adapts to a droplet profile with N = 193 in a few milliseconds after the release from the trap (see top insets in Fig. 3). In this case, the initial excitation energy is very small and indeed we observe a limited dynamics. Remarkably, a droplet is quickly formed even if the trapping frequency differs significantly from the optimal value, as in the case of ω/2π = 200 Hz (here with N = 159), which highlights the efficiency of the self-evaporation mechanism to dissipate the initial excitation energy.
To characterize the relaxation process towards a stationary droplet, here we indicate with N d i (t) (i = 1, 2) the number of atoms remaining within the droplet volume and we define a running value of N from Eq. (5) as where the second (approximate) equality holds for the current values of the scattering lengths. We indicate with σ(t) the droplet size (evaluated as the average rms of the density distributions of the two atomic components defined as in the previous section, see Appendix A), and with E d (t) the corresponding energy. The evolution of σ(t), E d (t), N d 1 (t)/N d 2 (t) and of the running value of N R (t) are plotted in Fig. 3, for ω/2π = 200 Hz and 600 Hz. We find that the relaxation process occurs via a sudden expulsion of atoms (at about t = 2 ÷ 4 ms, depending on ω), which allows to dissipate most of the initial excitation energy. From Fig. 3b we can estimate an average dissipation rate, by evaluating the decrease in energy occurring in the first part of the evolution, until the deviation from the equilibrium droplet energy becomes very small. We find that it varies from 0.3 to 2.5 MHz/ms, as ω/2π goes from 600 Hz to 200 Hz. Figure 3c shows that the ratio between the atom numbers in the two components always remains very close to the nominal equilibrium value. We recall that a droplet can sustain an excess of particles in one of the two components δN i /N i up to a critical value ∼ |δa|/a ii , beyond which the particles in excess are expelled [1]. Here the deviations of N d 1 (t)/N d 2 (t) are always within this tolerance (shaded area in Fig. 3c).
We can conclude that -in the absence of 3BL -the relaxation dynamics of a binary mixture produced by means of the non-adiabatic experimental protocol of [10]  is dominated by the self-evaporation mechanism, consisting in the dissipation of the initial excitation via the release of particles wave packets emitted from the droplet (see also Fig. 5 in Appendix A).

IV. DYNAMICS OF THE DROPLET FORMATION IN THE PRESENCE OF 3BL
Having now a clear idea of how self-evaporation works, we can discuss the dynamical formation of the droplet in the realistic experimental conditions of Ref. [10], where a significant role is played by 3BL. To do this, we use the same model as in Sec. III B, adding to each equation in (6) a non-unitary term where K iii are the intra-species 3BL rates [22]. Their values are K ′ 111 /3! = 9×10 −40 m 6 /s, K ′ 222 /3! = 1×10 −41 m 6 /s, where the primed values correspond to the loss rates of thermal atoms and the factor 1/3! accounts for the Bose statistics of condensed atoms [23]. Following Ref. [10], we assume that the two hyperfine levels are equally populated initially, N 1 = N 2 = 2 × 10 5 , and we fix the trapping frequency in Eq. (10) to a value of the same order of the geometric average of the experimental frequencies, namely ω = 2π × 200 Hz.
In the present case, owing to the complex dynamics that originates from the presence of 3BL, especially in the first part of the evolution, we use a different strategy to distinguish the droplet from unbound expanding atoms, similar to the protocol implemented in [10]. We integrate the 3D density profiles along one direction and we fit the column-density with a double gaussian function , where the first gaussian corresponds to the central droplet and the second one to the unbound expanding cloud.
In Fig. 4 we show the evolution of the droplet size σ(t), of the atom number in the droplet and in the expanding cloud surrounding it, of the droplet population ratio N 1 /N 2 , and finally of the droplet energy. In Fig. 4e we compare the measured droplet energy with the energy of the ground state for the same atom number E eq (N (t)). In the first stages of the evolution, the dynamics is dominated by losses of atoms from state |1 , which bring the ratio N 1 /N 2 closer to the nominal value of Eq. (1) and thus significantly reduce the droplet energy. These losses come both from 3B recombinations and from the release of atoms from the droplet into the unbound component, due to the imbalance in the interaction energies. In this stage, the system dissipates the exceeding energy pretty quickly: We measure an energy loss rate of about 20 MHz/ms. Notice that this is almost 10 times larger than the self-evaporation cooling rate measured in the previous section for the same initial ω. When the ratio N 1 /N 2 reaches its nominal value anyway, the droplet keeps losing atoms in |1 due to 3BL. To compensate for that, the droplet starts releasing atoms in |2 . This population dynamics, together with the compression of the atomic cloud, causes the two apparent bumps in Fig. 4e. At this point we find the system close to its equilibrium configuration, with the central cloud forming a metastable droplet with an atom number significantly smaller than the initial one. Since the system keeps losing atoms, due to 3BL in |1 and to population re-equilibration in |2 , we see that around t = 12 ms the droplet population drops below the critical value N c [1,9,10], so that the binding mechanism breaks and the cloud starts expanding. From this analysis we can conclude that the energy variations driven by 3BL and the corresponding population balancing occur on timescales much shorter than those typical of the self-evaporation mechanism described above, so that it is hard to determine if the latter occurs at all. Even if that were the case, its effect would be negligible with respect to the two leading loss mechanisms.
We have also investigated whether a different choice of the initial parameters, such as the trapping frequency or the population ratio N 1 /N 2 , would lead to a faster equilibration time and/or to longer droplet lifetimes. We find that the values of ω and N 1 /N 2 discussed above, which correspond to the values chosen in the experiment [10] by optimizing the droplet lifetime, represent indeed a convenient choice. In fact, we observe that a larger value of the initial trapping frequency ω leads to stronger 3BL in the first stages of the evolution, thus reducing even further the droplet lifetime. On the other hand, a ratio N 1 /N 2 closer to its equilibrium value a 22 /a 11 does not help either, because the stronger losses in |1 quickly drive the system away from this equilibrium condition, triggering an earlier release of atoms from state |2 .
To study experimentally the effect of self-evaporation, one would need to find a suitable attractive mixture with significantly lower 3BL. This might be the case of the 41 K-87 Rb mixture of the recent experiment in Ref. [11], where thanks to the different values of the scattering lengths the droplet forms at lower densities, so that 3BL are less effective. Indeed, there the lifetime is expected to be one order of magnitude larger than that of the 39 K mixture considered here, at least. In this case, as discussed in Sec. III B, it is preferable to prepare the initial state at the nominal equilibrium ratio N 1 /N 2 and as close as possible to the droplet density profile, to have short equilibration times and limited overall dynamics (see Fig. 2). Remarkably, if the droplet is prepared in the self-evaporation regime, any unwanted excitation would be then efficiently dissipated via that mechanism.

V. CONCLUSIONS
We have characterized, by means of numerical simulations, the dynamics of self-evaporation of an atomic quantum droplet, considering both the case of small and large initial excitations. We have verified that the regime of its occurrence corresponds to that predicted in Ref. [1], and we have discussed how and on which timescales this mechanism allows the dissipation of energy in the droplet. We have then simulated the same preparation of the droplet realized in Ref. [10], including the effect of 3BL. This has allowed to identify the relevant mechanisms involved in the dynamical formation of a self-bound droplet of 39 K. In this case, the evolution of the system is dominated by the presence of 3BL and by a continuous release of atoms to restore the proper population ratio, whereas the self-evaporation mechanism plays -if anyonly a negligible role. We have finally discussed the optimal strategies for the preparation of droplets close to their ground state, both for the experimental case considered here and for possible future experiments with reduced 3BL. The results reported here, which provide an accurate characterization of the dynamical evolution of excited droplets, are both useful for a deeper understanding of the recent experimental observations of quantum droplets and relevant for future experimental studies, especially regarding the effect of self-evaporation.
The dynamical formation of stationary (or metastable) droplets occurs via a relaxation process during which part of the atoms are expelled from the central volume. Therefore, a criterion for distinguishing between the droplet and the outgoing density waves, is needed. Here we exploit the radial symmetry of the present problem and we consider -along with the proper density distribution n(r, t) in radial coordinates (Fig. 5a,b) -the quantitȳ n(r, t) ≡ r 2 n(r, t) which includes the factor r 2 corre- Plot of the total density n(r, t) = 2 i=1 Ni|ψi(r, t)| 2 (a,c) and of the quantityn(r, t) ≡ r 2 n(r, t), at different evolution times (t = 4, 12 ms -here for the case ω/2π = 200 Hz of the simulations discussed in Sec. III B). The abscissa of the red dot (see text) represents the value of the droplet radius R d (t).
sponding to the Jacobian determinant in spherical coordinates (Fig. 5c,d). Therefore, a linear integration of n(r, t) corresponds to a volume integration of the density n(r, t). This representation also provides a convenient way to visualize the point at which the outgoing density waves "detach" from the inner volume, which we identify with the first minimum ofn(r, t) at the right-hand side of the bulk, marked by a red dot in Fig. 5c,d. The abscissa of this point is the value that we identify with the droplet radius R d (t). In the numerical simulations of Sec. III B for example, the droplet radius is initialized at R d (0) = 10 µm and it then evolves -both continuously or via sudden jumps -following the position of the leftmost minima (outside of the origin). The latter is typically found between 5 µm and 10 µm during the whole evolution considered here (up to 20 ms), corresponding to the fact that the bulk hosts a breathing mode and it does not expand. Contrarily, the density waves which are expelled from the droplet move outward while expanding, as shown in Fig. 5c,d.