Self-evaporation dynamics of quantum droplets in a 41K-87Rb mixture

We theoretically investigate the self-evaporation dynamics of quantum droplets in a 41K-87Rb mixture, in free-space. The dynamical formation of the droplet and the effects related to the presence of three-body losses are analyzed by means of numerical simulations. We identify a regime of parameters allowing for the observation of the droplet self-evaporation in a feasible experimental setup.

Besides being self bound, one of the most peculiar characteristics of these quantum objects is the fact that they can be self-evaporating.Namely, in certain regimes the droplet cannot sustain any collective mode, such that any initial excitation is completely dissipated (that is, evaporated) until the system relaxes into a droplet with a lower number of atoms.This remarkable feature, originally predicted in Ref. [1] for bosonic mixtures, has been so far elusive to the experimental detection.Indeed, the first pioneering experiments performed with bosonic mixtures [10,11] have shown that the homonuclear mixtures of 39 K suffer from strong three-body losses that continuously drive the system out of equilibrium, and eventually lead to the depletion of the droplet.This behavior was then confirmed by the theoretical analysis reported in Ref. [28], where the complex dynamics taking place during the droplet formation is thoroughly described.There, it was shown that the evolution of the system is indeed dominated by the presence of three-body losses and by a continuous release of atoms to restore the proper population ratio, whereas the self-evaporation mechanism plays only a negligible role, if any.
In this respect, the recent experiment in Ref. [14] represents a promising setup in which the self-evaporation mechanism could be properly investigated.As a matter of fact, for a 41 K-87 Rb mixture the regime of parameters for which self-bound droplets form is such that three-body losses are expected to be significantly suppressed.Indeed, in such a mixture droplets form at lower densities n ∼ (δg/g) 2 a −3 , and this allows for a larger ratio τ life /τ ∼ n −1 , with τ life and τ being respectively the lifetime (limited by three-body losses) and the characteristic time scale of the droplet dynamics [14].Thus, the fact that the 41 K-87 Rb mixture is characterized by larger scattering lengths with respect to the 39 K mixture employed in Refs.[10,11] allows for lower densities and, therefore, longer lifetimes.Actually, in that experiment no appreciable effect of three-body losses was observed on the timescale of several tens of milliseconds.
Motivated by the previous discussion, in this paper we analyze the self-evaporation mechanism and the role of three-body losses on the collective excitations of a 41 K-87 Rb quantum droplet.For the sake of conceptual clarity, and in view of the fact that quantum droplets at equilibrium are spherically symmetric objects, we shall restrict the analysis to the case of a spherically symmetric system.Though this assumption limits the study to the monopole collective mode, it is however sufficient to draw the general behavior of the system and to discuss the role of three-body losses.
The paper is organized as follows.In Section 2 we review the general formalism for describing quantum droplets in heteronuclear bosonic mixtures.Then, in Section 3 we discuss the preparation of the initial state, namely a droplet compressed by a harmonic trapping.The dynamics of the droplet after the release of the trap is then studied in Section 4 in the unitary case (no losses), and in the presence of three-body losses.Finally, in Section 5 we draw the conclusions.

Self-bound droplets
We consider a binary condensate of 41 K and 87 Rb atoms in the |F = 1, m F = 1 state, as considered in the recent experiment in Ref. [14].The two components will be indicated as 1 and 2, respectively.This system is described by the following Gross-Pitaevskii (GP) energy functional, including both the mean field term and the Lee-Huang-Yang (LHY) correction accounting for quantum fluctuations in the local density approximation: [29] g ij n i (r)n j (r)dr + E LHY (n 1 (r), n 2 (r))dr , (1) where m i are the atomic masses, V i (r) the external potentials, and n i (r) = |ψ i (r)| 2 the densities of the two components (i = 1, 2).The LHY correction reads [1] with κ = 8m 3/2 1 /(15π 2 h3 ) and f (z, u, x) > 0 being a dimensionless function of the parameters z ≡ m 2 /m 1 , u ≡ g 2  12 /(g 11 g 22 ), and x ≡ g 22 n 2 /(g 11 n 1 ) [1,29].The mixture is completely characterized in terms of the intraspecies g ii = 4πh 2 a i /m i (i = 1, 2), and interspecies g 12 = 2πh 2 a 12 /m 12 coupling constants, where m 12 = m 1 m 2 /(m 1 + m 2 ) is the reduced mass.The values of the homonuclear scattering lengths are fixed to a 11 = 62a 0 1 , a 22 = 100.4a0 , whereas the heteronuclear scattering length a 12 is considered here as a free parameter, that can be tuned by means of Feshbach resonances [14].The onset of the MF collapse regime corresponds to δg = g 12 + √ g 1 g 2 = 0, at a c In free space (V i ≡ 0), the equilibrium density of a droplet is obtained by requiring the vanishing of the total pressure [31], which yields [1] Following [1,29], we consider this function at the mean-field collapse u = 1, f (z, 1, x).We note that the actual expression for f can be fitted very accurately with the same functional form of the homonuclear case [19] For a finite number of atoms the droplet has a finite size, and it can be effectively described by a single wave function that satisfies the following dimensionless equation with |φ 0 | 2 dr = 1, and where Here ξ represents the length scale of the droplet [1,25] i the equilibrium density of each component in the uniform case.We also remind that the wave functions of the two condensates forming the droplet are given by ψ i = n (0) i φ 0 .

Preparation of the initial state
Following Ref. [14], we initially prepare the atomic mixture in the ground state of an optical dipole trap.In order to simplify the discussion, and motivated by the fact that quantum droplets at equilibrium are spherically symmetric objects, we assume i) that the two condensates are initially prepared in a spherically symmetric potential U d j (r) (j = 1, 2), and ii) that the differential vertical gravitational sag (due to the different masses of the two atomic species) can be exactly compensated.With these assumptions, which will be maintained throughout this work for easiness of calculations and conceptual clarity, both the ground state and the dynamical evolution can be obtained by solving spherically symmetric equations.We remark that this approach restricts the analysis to the excitations of the monopole mode only.Indeed, the excitation of surface modes (with angular momentum = 0), which arise when the condensate is prepared in an asymmetric trap, is ruled out by the assumption of spherical symmetry.However, the presence of these modes is not expected to modify qualitatively the picture resulting from the following discussion.
The ground state of the system is obtained by minimizing the energy functional in Eq. ( 1) by means of a steepest-descent algorithm [32].Formally, this corresponds to solve the following set of stationary Gross-Pitaevskii (GP) equations for two wave functions ψ j [33] where ∇ 2 r represents the radial component of the Laplacian: The local chemical potentials µ i (n 1 , n 2 ) include both the usual mean-field term and the Lee-Huang-Yang (LHY) correction [1], namely with We recall that the optical dipole potential experienced by the two atomic species can be written as U d j (r) = U 0j I(r), where U 0j is species-dependent (j = 1, 2) and it depends on the atomic species polarizability and the wavelength of the optical potential, and I(r) represents the intensity of the laser beam.For the atomic mixture 41 K-87 Rb and a wavelength λ = 1064 nm [14] one has α ≡ U 02 /U 01 1.15.Then, in the harmonic approximation, we have with ω 1 = ω 2 m 2 /(αm 1 ).In the following, the value of the trap frequency will be considered as a free parameter, to be tuned in order to produce the desired amount of initial excitations.For low enough ω j , the droplet will be prepared close to the equilibrium configuration in free space, whereas a tight confinement obviously corresponds to the excitation of a compressional mode.
The other free parameters of the system are the inter-species scattering length a 12 , that we will consider in the range a 12 ∈ [−95, −73.6]a 0 , below to the mean-field collapse threshold, and the number of atoms of the two species.As for the latter, we shall vary N 2 in the range [1.5, 15] × 10 4 , keeping N 1 = N 2 g 22 /g 11 , so that the atom numbers match the nominal equilibrium ratio in Eq. (4).

Dynamics
Once the initial state is prepared, the optical dipole trap is switched off and the system is let to evolve.The droplet dynamics is studied here by solving the following set of GP equations [28,34,35] in the presence of three-body losses, which are included by adding to the energy functional in Eq. ( 1) a dissipative term −(i/2)hK 3 n 1 (r, t)n 2 (r, t) 2 d 3 r (with K 3 = 7 × 10 −41 m 6 /s, see Ref. [14]) accounting for the dominant recombination channel, i.e.K-Rb-Rb.Since the droplet is prepared in a compressed configuration (owing to the presence of the trap), the droplet will start oscillating.According to discussion in Ref. [1], in the absence of three-body losses the dynamics is expected to be characterized either by sinusoidal oscillations of the droplet width, where the monopole mode exists, or by damped oscillations, in the so-called self-evaporation regime.The latter represents one of the remarkable properties of quantum droplets [1,28], and it takes place in a certain window of N, where the excitation spectrum of the droplet lies entirely in the continuum.The nominal phase diagram for our system (obtained from the predictions Ref. [1]) is shown in Fig. 1 as a function of the interspecies scattering length a 12 and of the (initial) number of atoms N 2 .The specific combination of parameters used in the figure, ( N − N c ) 1/4 , is introduced after Ref. [1].The monopole mode is expected to be stable for N > 933.7, and to evaporate for lower values of N. In the window 94.2 < N < 933.7 other modes with = 0 (not included in the present discussion) may appear, whereas below N < 94.2 neither the monopole nor the surface modes ( = 0) can be sustained, such that the droplet is expected to evaporate any initial excitation [1].Below N c , where a droplet cannot be formed, the mixture forms a so-called LHY fluid: in this regime the MF interactions almost cancel out (δg/ √ g 11 g 22 ≈ 0), and the system is governed only by quantum fluctuations [18][19][20].

Estimate of the droplet lifetime
Before analyzing the detailed dynamical behavior of our system it is convenient to discuss the role of three-body losses on the droplet lifetime.To this aim, we shall consider a droplet in free space, at equilibrium (at t = 0).Then, assuming that the shape of the droplet is weakly affected by the atom losses (which is justified at the early stages at the evolution, at least), we can write where N 1 (0)/N 2 (0) = g 22 /g 11 ≡ γ , and with ρ(r) representing the droplet density profile normalized to unity [see Eqs. ( 4) and (7)].With this in mind, the evolution of N i (t) can be obtained from the previous Eqs.( 13) by neglecting the kinetic and chemical potential terms [which concur in determining the shape ρ(r)], left multiplying each equation by ψ * i , and integrating over the volume (see also Ref. [36]), yielding where we have defined Ni (τ) ≡ N i (τ)/N i (0) and τ ≡ tK 3 N 2 2 (0) ρ 3 d 3 r.Notice that the factor of 2 in Eq. ( 16) corresponds to the fact that here we are considering the dominant recombination channel K-Rb-Rb (two atoms of Rb are lost for each atom of K).The above equations have an approximate solution of the form which is very accurate, indeed.In the present case (γ 0.8736), from a fit of the exact numerical solution of Eqs. ( 15) and ( 16) we find τ 1 0.70, α 1 0.88, β 0.43, τ 2 0.71, α 2 0.86.In particular, τ 2 corresponds to the half-life of the i = 2 component (regardless of the value of α 2 ), in which losses are dominant.Therefore, it represents a characteristic time through which we can measure the impact of three-body losses on the lifetime of the droplet.In our simple model the half-life is therefore t 1/2 τ 2 / K 3 N 2 2 (0) ρ 3 d 3 r that, besides the explicit dependence on N 2 (0), also depends implicitly on a 12 through the density distribution ρ(r).The behavior of t 1/2 as a function of a 12 and N 2 (0) is shown in Fig. 2.There we compare the prediction of the above analytical model with the actual values obtained from the solution of the GP equations in (13).The qualitative agreement is remarkable.

Damped monopole oscillations
Given the above picture, in the following we shall investigate the self evaporation dynamics for a 12 = −82a 0 , where the droplet half-life is larger than 50 ms [see Fig.  18); (b) values extracted from actual decay of N 2 (t) obtained from the solution of the GP equations in (13).The color scale is saturated at t 1/2 = 100 ms (dark blue).The dashed-dotted line corresponds to t 1/2 = 50 ms.The two circles at a 12 = −82a 0 refer to the parameter configurations considered in the GP simulations of Sec.4.2.
monopole mode is expected to be stable (see Fig. 1) the droplet is affected by severe three-body losses (see Fig. 2), which make the detection of this mode unfeasible.
In order to distinguish between the atoms remaining in the droplet and those that evaporate, we define the droplet volume as that contained within a certain bulk radius.In the present case, it can be conveniently fixed to R d = 8 µm 2 .Accordingly, with N d i (t) we indicate the number of atoms of each species (i = 1, 2) within the droplet volume, at time t.Then, we define a running value of N as [28] N R (t) ≡ 1 where k 4.41 × 10 −3 for the current values of the scattering lengths.The evolution of the system is shown in in Fig. 3, where we plot the rms width σ(t) of the droplet, the running value of , and the fraction of evaporated atoms for each species N ev i (t)/N i , with and without three-body losses.In both cases the mixture is initially prepared in the ground state of a dipole trap of frequency ω 2 = 2π × 50 Hz (we recall that ω 1 = ω 2 m 2 /(αm 1 ), see Sec. 3; in the present case ω 1 = 2π × 67.9 Hz).
Let us first focus on the clean case, in the absence of three-body losses.As expected, we find that in both the investigated cases the droplet width performs damped sinusoidal oscillations (corresponding to a damped monopole mode), and the system eventually relaxes to an equilibrium configuration [see Fig. 3(a), along with the top panels], corresponding to a smaller, stationary value of N R (t), see panels (b).This decrease in the number of particles is a consequence of the self-evaporation mechanism [1,28], which takes place within the first 15 ms.Looking at panel (d), where we plot the fraction of atoms of each species that are lost by self-evaporation, it is clear that the values of N ev i (t)/N i soon reach their asymptotic value (modulo small fluctuations).In addition, also the ratio between the atom numbers in the two components remains close to the equilibrium value, see panels (c), with deviations below the reequilibration threshold 2 The numerical simulations of the GP equations are performed on a computational box that is at least one order of magnitude larger than R d .(shaded area in the figure).Indeed, we recall that a droplet can sustain an excess of particles in one of the two components δN d i /N d i up to a critical value ∼ δg/ √ g 11 g 22 [1] (≈ 11% in the present case), beyond which particles in excess are expelled.
Let us now discus how the presence of three-body losses affects the above picture.From panels (b) it is evident that losses produce a continuous drain of particles from the droplet.Nevertheless, after 100 ms of evolution the value of N R (t) is still well above the critical value N c for the existence of a droplet.Indeed, after several tens of milliseconds the mixture still forms a self-bound droplet (see e.g. the density profile at 50 ms) which keeps undergoing damped monopole oscillations, as shown in panels (a).Notably, the initial stage of the evolution is still dominated by self-evaporation, see panel (d), and then both the frequency and the amplitude of the oscillations become larger than in the clean case (without losses), because of the atoms that leave the bulk and populate the tails.Actually, there are two opposite effects taking place: on the one hand the droplet shrinks because of the loss of atoms due to three-body recombination, on the other hand there is an outward flow of K atoms that evaporate from the droplet, see panel (d).This species selective evaporation is due to the fact that the major loss of atoms affects the Rb component (see previous section), so that a progressively increasing fraction of K atoms cannot stay bound inside the droplet, and it is let free to expand.We remark that this is obviously a non-equilibrium process, as it is evident from the fact that N R (t) does not relaxes to a stationary value and that the ratio N d 1 (t)/N d 2 (t) soon run over the reequilibration threshold, anyway.All this is responsible for the different 'asymptotic' behavior of σ(t) in the two cases shown in panels (a).For N 2 = 4 × 10 4 (left), the two effects approximately compensate each other, and the width oscillates close to the nominal equilibrium value without losses.Instead, for N 2 = 1.5 × 10 5 (right), which is characterized by a shorter lifetime, the droplet width decreases progressively (though keeping oscillating).It is worth to remark that the actual behavior is also sensitive to the choice of the droplet radius R d .

Conclusions
We have theoretically investigated the self-evaporation dynamics of quantum droplets in a 41 K-87 Rb mixture in feasible experimental setups, including the effects of three-body losses.The mixture is prepared in the ground state of a spherically symmetric harmonic trap, that is then released thus letting the system evolve in free space.The subsequent dynamics, characterized by the excitation of the monopole breathing mode, has been analyzed by solving the coupled Gross-Pitaevskii equations for the two components.For the estimated values of three-body losses (K 3 = 7 × 10 −41 m 6 /s for the dominant recombination channel K-Rb-Rb [14]), we find that by tuning the interspecies scattering length a 12 , the lifetime of the system can be easily adjusted to be of the order of, or larger than, 100 ms.This makes 41 K-87 Rb droplets much more robust than those realized with two hyperfine states of 39 K, whose lifetime is limited to the order of 10 ms [10,11,28].Such long lifetimes permits to follow the droplet dynamics for several tens of milliseconds, without any appreciable loss of resolution.In this scenario, we have found that the initial stage of the evolution is dominated by the self-evaporation mechanism even in the presence of three-body losses, and that the latter induce an interesting non equilibrium dynamics at later times.These findings make the experimental investigation of collective modes of self-bound droplets in 41 K-87 Rb mixtures very promising.

Figure 1 .
Figure 1.Heat map of ( N − c ) as a function of a and N 2 , with N c = 18.65 being the critical value for the existence of a droplet [1].The dotted-dashed lines, corresponding to N = 933.7,94.2, 18.65 (from left to right), represent the boundaries between the different regimes (see text).The vertical red dashed-dotted line at a 12 = −73.6 corresponds to the onset of the MF collapse.The two circles at a 12 = −82a 0 refer to the parameter configurations considered in Sec. 4.

Figure 2 .
Figure 2. Heat map of t 1/2 as a of a 12 and N 2 (0).(a) Prediction of the analytical model, see Eq. (18); (b) values extracted from actual decay of N 2 (t) obtained from the solution of the GP equations in(13).The color scale is saturated at t 1/2 = 100 ms (dark blue).The dashed-dotted line corresponds to t 1/2 = 50 ms.The two circles at a 12 = −82a 0 refer to the parameter configurations considered in the GP simulations of Sec.4.2.

Figure 3 .
Figure 3. Evolution of (a) the droplet width σ(t), (b) the running value of N R (t) (see text), (c) the ratio N d 1 (t)/N d 2 (t), and (d) the fraction of evaporated atoms for each species N ev i (t)/N i , after the release from a trap of frequency ω 2 = 2π × 50 Hz, with and without three-body losses.The insets in the top row show the total density of the binary mixture, n(r, t) = ∑ 2 i=1 N i |ψ i (r, t)| 2 , at different evolution times, corresponding to the red circles in (a).The horizontal line in (b) represent the nominal equilibrium value N 1 /N 2 = g 22 /g 11 0.873, and the dashed area the corresponding tolerance (see text).Here a 12 = −82a 0 , N 2 = 4 × 10 4 (left), N 2 = 1.5 × 10 5 (right).