Breathing mode in two-dimensional binary self-bound Bose gas droplets

In this work, we present the study of the stationary structures and the breathing mode behavior of a twodimensional self-bound binary Bose droplet. We employ an analytical approach using a variational ansatz with a super-Gaussian trial order parameter and compare it with the numerical solutions of the extended GrossPitaevskii equation. We find that the super-Gaussian is superior to the often used Gaussian ansatz in describing the stationary and dynamical properties of the system. We find that for sufficiently large non-rotating droplets the breathing mode is energetically favourable compared to the self-evaporating process. For small self-bound systems our results differ based on the ansatz. Inducing angular momentum by imprinting multiply quantized vortices at the droplet center, this preference for the breathing mode persists independent of the norm.


I. INTRODUCTION
For bound systems to form without external confinement, a balance between repulsive and attractive forces is required, this principle holds for systems ranging from water droplets to atomic nuclei and metallic clusters. For ultra-cold bosonic atoms the possibility to form self-bound droplets was proposed for binary gases with suitably tuned contact interactions [1,2], where higher-order terms in the total energy density functional may become sizeable. In fact, such Bose droplets were first discovered for dipolar bosonic systems [3][4][5][6][7][8], where the effective interatomic interactions may be adjusted such that self-bound states become stabilized by quantum fluctuations [6,[9][10][11][12][13][14]. The experimental discovery of droplets in dipolar gases was soon followed by their realization in binary Bose gases [15,16], implementing Petrov's original proposal [2]. The quantum fluctuation contributions to the total energy density functional are effectively represented by the Lee-Huang-Yang (LHY) terms [17], extending the usual mean field (MF) Gross-Pitaevskii equation (see [18,19] for a review). Usually, the LHY terms can be neglected as their contribution to the total interaction energy is negligible. However, in the case of a binary condensate where the MF interactions between components are tuned such that they become very small, the LHY contributions may dominate the system's properties [2]. In a binary three-dimensional condensate the LHY terms are positive definite [17] and thus purely repulsive. A three-dimensional binary condensate on the verge of collapse in the MF description can thus be stabilized via the repulsive LHY terms. In lower dimensions, however, the corrections can take on an attractive or repulsive form [20][21][22]. The proposal and realization of self-bound Bose gas droplets sparked recent research efforts to set focus on different aspects such as vorticity embedded in dipolar [23,24] or contact-interacting systems [25][26][27][28][29], supersolid behavior in dipolar and spin-orbit coupled systems [30][31][32][33][34][35][36][37][38] or collective excitations [39][40][41][42][43]. For trapped gases, the elementary excitation modes have been of fundamental interest ever since the celebrated experimental realizations of atomic Bose-Einstein * philipp.sturmer@matfys.lth.se condensates [44][45][46]: Examples range from a hydrodynamic description [47] to the Bogoliubov-de Gennes equations [18,19,48] and to solving the Euler-Lagrange equations [49][50][51][52]. Breathing mode oscillations were discussed in terms of the underlying symmetry properties of a Bose gas in two dimensions [53], and semiclassical methods have been probed to describe low-lying excitations in the limit of large boson numbers [54]. Collective excitations of BECs carrying vortices or vortex lattices have also been studied extensively; to mention only a few examples from the vast literature, see Refs [53][54][55][56][57][58][59][60][61][62]. More recently and in the context of Bose droplets, in Ref. [63] a variational approach for solving the Euler-Lagrange equations was applied to calculate the breathing mode in a harmonically confined three-dimensional Bose gas where the MF interactions are completely cancelled, such that the LHY interaction is effectively the only interaction acting on the gas, creating a so-called LHY fluid. Excitations compete with the self-evaporative process, naturally occurring in self-bound SBose gases due to a non-positive chemical potential µ 0, and are generally only observable if the associated excitation energy is lower than |µ|. However, this requirement can stay unfulfilled for selfbound systems [2] depending on the system's dimensionality, if the excitation energy exceeds |µ|, see Refs. [40,42]. Recent works on collective excitations in self-bound systems have mainly focused on one-dimensional droplets [20,39] or on the cross-over between the soliton and droplet regime [39]. Besides solving the Bogoliubov-de Gennes equations in one [42] and in three dimensions [2,43], a Gaussian trial order parameter is often utilized to describe the properties of the droplet state [39,40,43]. However, due to the self-bound nature and the resulting flat-top profile of the droplet, the classical Gaussian shape proves to be a poor approximation of stationary properties of droplets. Yet the Gaussian shape still represents the dynamics reasonably well in one dimension [39], although overestimating the results in three dimensions [43]. In search for a better approximation of the order parameter, a 1 − tanh(x) function in the Thomas-Fermi limit has been utilized [29], which imposes difficulties in an analytical variational ansatz. To overcome the limitations by these functions, we here study the breathing mode and stationary properties using a super-Gaussian ansatz for droplets in two dimensions, following Refs. [41,64] for one-dimensional systems. The super-Gaussian allows for a more precise description of a flat-top shape than the pure Gaussian, while still allowing for an easier analytical treatment compared to a 1 − tanh(x) function.
In continuation of our previous work [26], in this paper we investigate the lowest-lying collective excitations of self-bound Bose-gas droplets in two dimensions. These collective excitations are generally referred to as (monopole) breathing modes, where the radial extent of the droplet oscillates unidirectionally. The work is structured as follows. In Sec. II we introduce the basic model for a two-dimensional binary self-bound Bose gas droplet. In Sec. III we use a super-Gaussian function as a trial order parameter and solve the Euler-Lagrange equations for complex fields of droplets. The calculation is followed by an analysis and comparison of stationary properties to numerical results self-bound systems in Sec. III A and the breathing frequency in Sec. III B. We continue in Sec. IV by studying droplets carrying angular momentum in the form of a phaseimprinted singly-or multiply-quantized vortex at the droplet center. We analyze the breathing mode frequency and energetically compare the mode with the self-evaporation process.

II. MODEL
Let us now consider a binary Bose-Einstein condensate in two dimensions, and label the two components by the indices 1 and 2. The atoms are chosen to have equal mass and equal intra-component scattering lengths a 11 = a 22 = a, with an inter-component scattering length denoted by a 12 . The number of atoms in each component is set to be equal and consequently one also has equal atom densities n 1 = n 2 = n. The two-component system thus reduces to effectively only one component, where the MF interaction cancels and only the LHY terms beyond mean field remain, and the time-dependent extended Gross-Pitaevskii equation (eGPe) takes the form see Ref. [20]. Here, ψ = ψ(r, t) is the order parameter with r = (x, y). As in Refs. [20,25,26,28,34], the scaling invariances have been used to bring the system into dimensionless form, such that it conveniently only depends on the norm N = d 2 r|ψ(r)| 2 . Given Eq. (1) the energy density E[ψ] of the self-bound system is then given by For the comparision of the variational ansatz discussed in the following, Eq. (1) is solved numerically using the Fourier split-step method in imaginary and real time.

III. SUPER-GAUSSIAN VARIATIONAL ANSATZ
In this section we use a super-Gaussian trial order parameter to calculate the stationary and dynamic properties of a self-bound binary droplet by taking advantage of the Euler-Lagrange equations for a complex field. We then continue by comparing the analytical solution to the numerical solution of the eGPe in Eq. (1). Recent works which made use of a variational analysis of such self-bound states used either a Gaussian ansatz for one or three-dimensional droplets [39,40,43], or a √ 1 − tanh x ansatz in the Thomas-Fermi approximation for two-dimensional droplets [29]. While the Gaussian ansatz for the breathing frequency yields small errors compared to the numerical eGPe solution for one-dimensional problems [39], it has limitations to properly describe the spatial properties of self-bound systems as it lacks the characteristic flat-top shape for high norms. Furthermore it overestimates the breathing mode frequency in three dimensions [43]. On the other hand, a √ 1 − tanh x function accurately describes the spatial properties, but faces limitations to describe the system for small norms [29] and also becomes analytically cumbersome when treated with the Euler-Lagrange equations. In comparison, a super-Gaussian ansatz offers a good middle ground, as it can be analytically handled and describes the low-and high-norm systems properly [41,[64][65][66]. For the following calculation we thus use the super-Gaussian as a trial order parameter. Assuming a circular droplet shape, we write with radial coordinate r, real amplitude A(t), chirp b(t) and width R(t). The real and positive exponent m is determined as a function of the norm N . The dynamics of a complex field, such as the reduced twodimensional two-component Bose gas for some generalized coordinate q i = {A(t), b(t), R(t)}, emerges from the respective Euler-Lagrange equations where the Lagrangian L is calculated by averaging the Lagrangian density L over the full space The Lagrangian density L for the two-dimensional droplet system results from the energy density for a self-bound system in Eq. (2) as Normalizing the wave function gives the number of particles in the condensate, where Γ(m) is the gamma function. We continue to write After evaluation of the integral in Eq. (5) by inserting the trial wavefunction in Eq. (6) and eliminating powers of A in favor of the norm N , the Lagrangian reads Thus, the Euler-Lagrange equations for the variational param- Let us now first determine m as a function of N . We extremize the Lagrangian L with respect to m, which gives the transcendental equation with the analytical exact value m exact For large N we have an approximation m approx m approx ≈ N 2π (ln 4 − 1).
In Fig. 1 we plot m according to Eqs. (11) and (12) as dashed and dotted-dashed lines respectively, as well as the difference ∆m = m exact − m approx (see inset). The difference ∆m approaches 0 as N → ∞, so that we can use the approximate expression in Eq. (12) for droplets with large N . Looking at the minimum of the effective potential U eff , as defined in Eq. (10), gives us an expression for the maximum density n m and the radial extent of the condensate R 0 , such that For large norms the maximum density approaches 1/ √ e as expected for this choice of scaling [25][26][27]34]. The mean radius r 2 depends on the width of the cloud, which simply reduces to the Gaussian case r 2 = R 2 for m = 1. Furthermore, the width of a stationary condensate is The low-lying excitation frequencies are found by assuming a single sinusoidal oscillation in R(t) with frequency ω 0 , such that the breathing frequency ω 0 is given by and one obtains The result in Eq. (17) can be divided into a low and high norm N regime. In the low norm regime, the behavior is dominated by the exponential behavior of n m , which approaches a constant value for high norm, such that the 1/N dependency governs the breathing frequency's behavior. From the time-independent eGPe, Eq. (1), we obtain the chemical potential µ as With N → ∞, µ then approaches −1/(2 √ e) as suggested by Refs. [20,25].

A. Stationary properties
Let us now use the previously obtained expressions for the exact and approximative value of m as well as the Gaussian case (m = 1) and compare them to the numerical solution of the eGPe equation, Eq. (1). In the following, unless otherwise specified, we will use solid lines for the Gaussian case with m = 1, dashed lines and dotted-dashed lines for m according to Eq. (11) and (12), respectively, and numerical results of the eGPe in bullets for norms 1 ≤ N ≤ 1000. We start by comparing the mean radius r 2 in Fig. 2. As expected from our previous discussion about m, the exact and approximate descriptions perfectly overlap for large N . With decreasing N the approximate description starts to deviate and our variational approach is well approximated by the Gaussian form of the order parameter for m = 1. For all values of N the numerical results obtained by solving the eGPe coincide directly with m exact solution. For the case of a large norm, it follows the behavior expected of a growing droplet. For a smaller norm, however, the system's properties are dominated by the kinetic energy. This bears some similarity to the bright soliton-like structures known to occur in one-dimensional systems [67][68][69][70].
As mentioned above, the maximum density n m approaches 1/ √ e and m approx in Eq. (12) gives the correct behavior in the limit of large N , as it can be seen in Fig. 3. With increasing N the Gaussian with m = 1 starts to deviate strongly from the super-Gaussian and numerical solution of the selfbound system. For small N , however, m = 1 approximates the super-Gaussian reasonably well, with only small deviations, while m approx differs significantly as expected. All our numerical values lie within a 2% deviation directly below the super-Gaussian with m exact . Having analyzed the mean radius r 2 and the maximum density n m we now look at the density profile, comparing the exact super-Gaussian, the Gaussian and the numerical solutions of the eGPe in Fig. 4. The dashed lines represent the Gaussian with m = 1, and the full lines the super-Gaussian. As already anticipated from the previous discussion, for small N both solutions coincide with each other. For increasing N the deviations increase and the super-Gaussian solution slowly approximates the characteristic flat-top shape.
Finally, let us look at the chemical potential µ, as specified it in Eq. (18) above. Due to µ 0 in self-bound condensates, an underlying self-evaporation process exists. Thus, in prin- ciple no other mode with energy E should be observable as long as E/|µ| 1. In a self-bound droplet, this may not be always the case, depending on the system's dimensionality, as shown in [2,20,40]. As before we compare the exact and approximate solutions for m, the Gaussian solution with m = 1 and the numerical results in Fig. 5. While the exact and approximate solutions converge towards the theoretical value of 1/(2 √ e) for high N , the Gaussian solution deviates significantly. Again, the numerical solution of the eGPe, Eq. (1) agrees well with the variational approach both in the limits of small and large norms N .
FIG. 5. Chemical potential µ of the droplet based on the variational and numerical solutions, as specified by the legend. The super-Gaussian ansatz correctly predicts the large N limit, while the Gaussian ansatz clearly deviates from it. For small N the approximate approach overestimates the chemical potential, while the Gaussian offers a good fit.

B. The breathing mode in a two-dimensional droplet
We now compare the analytical results for the breathing mode frequency ω 0 in Eq. (17) to the numerical results obtained by solving the eGPe, Eq. (1), for 1 ≤ N ≤ 1000. We excite the breathing mode by introducing a small interaction perturbation, multiplying the interaction term in Eq. (1) with a (dimensionless) factor k = 1.001 which is decreased back to unity linearly in a time interval t 0 = 500. The interaction perturbation k is chosen such that the size modulations during the time propagation are small compared to r 2 . The breathing frequency ω 0 is measured via Fourier-analysis of r 2 . We then calculate ω 0 /|µ| via the analytical chemical potential in Eq. (18) and Fig. 5, and also compare to the numerical solution. Fig. 6 shows ω 0 in comparison to analytical results with varying m, and in comparison to the numerical solution of the eGPe. For large N the approximate and exact solution coincide and give the correct behavior for ω 0 as we can confirm from the comparison to the numerical solution of the eGPe. However, when approaching the low norm regime the numerical results start to deviate strongly from the analytical solutions. This is explained by an emerging beating pattern in the numerical results, such that there are several significant breathing mode frequencies. We assume however in Eq. (16) that the oscillating system only produces one single frequency. Similar to our observation for a two-dimensional droplet, the deviations between the Gaussian ansatz and the numerical results of the eGPe have also been found for a three-dimensional system [43]. However, such deviations between analytical and numerical approaches are absent for the one-dimensional selfbound system [39]. With the chemical potential in Eq. (18) we calculate the ratio ω 0 /|µ| shown in Fig. 7. Surprisingly, all three analytical solutions overlap for high N , while the eGPe solution starts to deviate as ω 0 /|µ| → 1. From the analytical results, we obtain that the breathing mode is observable from N 36, while the numerical solution of the eGPe suggests N 21.

IV. BREATHING MODE IN DROPLETS WITH ANGULAR MOMENTUM
We continue our analysis by studying systems with angular momentum, obtained numerically by starting from initial states of the form ψ 0 (r) = Cr S exp (−αr 2 + iSθ), The size of the droplets and minimum norm N for stability follow the heuristically derived relation in Ref. [25]. The phase in cases with non-zero angular momentum is shown as small insets in each figure.
FIG. 9. Frequency spectrum for droplets with angular momentum as specified in the legend. With increasing L/N the breathing frequency ω0 decreases. The onset of each branch is given by the realtime stability of the system. For the multiply quantized systems, the stability requirement for the minimum required norm N is equal to that found in [25].
where C is a normalization constant, r and θ radial and angular coordinate, α some positive number, and S is the vorticity [25]. Starting the time evolution from the states obtained this way, we then bring the interaction strength back to unity and obtain the breathing mode frequency ω 0 via a Fourieranalysis of r 2 and the chemical potential µ. Here, we choose systems with L/N = 1.0, 2.0, 3.0 as shown in Fig. 8 for N = 850. As stated before, the breathing mode will only be observable if the condition ω 0 < µ is fulfilled. Due to the small oscillations of r 2 when measuring ω 0 we can calculate µ numerically at any given point in the real-time evolution. Figs. 9 and 10 show the breathing frequency ω 0 and the ratio ω 0 /|µ| for the systems as described above in the range 1 ≤ N ≤ 1000. Both quantities, ω 0 and ω 0 /|µ| follow the same pattern, where the droplet without angular momentum corresponds to the highest lying branch followed by the multiply-quantized vortices with lower ω 0 as L/N increases.
With increasing L/N the droplet requires a higher norm N in order to support the size of the vortex core, as otherwise the modified kinetic energy due to angular momentum overcomes the attractive two-dimensional LHY term. We observe this behavior in Fig. 9 and 10 as the onset of the branches for the respective L/N values. Further the onset follows the minimum N requirement for a droplet to support a vortex with L/N as in Ref. [25]. As the numerical solution is the same as in Fig. 7, the breathing mode continues to be observable from N 21 while the analytical solution in Fig. 7 becomes observable for N 36. This is in stark contrast to the one-dimensional case where the breathing mode is observable at any norm N [42] and the quasi one-dimensional case [40] where it is observable for most configurations.

V. CONCLUSIONS
We have studied the breathing mode for symmetric twodimensional self-bound binary Bose-gas droplets with components in the non-rotating case as well as for systems with a multiply-charged vortex imprinted at the center. For the breathing mode of the droplet with L/N = 0 we employed a variational super-Gaussian ansatz [41] and compared it to the widely used Gaussian ansatz [39,40,43] as well as the numerical solution of the eGPe. We found that the super-Gaussian supersedes the Gaussian in predictability of stationary and dynamical properties of the self-bound system. However, for lower norms there is a great discrepancy between the numerical solution of the eGPe and the analytical solutions for the breathing mode, such that the ratio ω 0 /|µ| barely exceeds 1 when the analytical solution suggests that ω 0 /|µ| > 1. Unlike two-and three-dimensional droplet [43] systems, such discrepancy between the numerical and analytical solutions is not obtained for the one-dimensional system [39]. This suggests that additional considerations should be taken into account for the real-time evolution of the higher dimensional weakly interacting system within the numerical eGPe framework, and can be subject for future studies. We furthermore found that the breathing mode in self-bound droplets without angular momentum can be observed from N 36, for which ω 0 /|µ| < 1. This is similar to the quasi one-dimensional and three-dimensional case [2,40,43]. However this behavior is different from that of a one-dimensional system, for which the breathing mode is observable at any N [39]. We continued our investigation by introducing angular momentum in the shape of a multi-charged vortex imprinted at the droplet center. For this metastable system [26] we find that with increasing L/N the breathing frequency decreases, and droplets with angular momentum have ω 0 /|µ| < 1. Furthermore, the stability requirement in N for droplets to support a vortex with L/N is the same as in [25]. The current work may be continued by further investigating the discrepancy between the numerically obtained breathing frequencies and the analytical solution. The super-Gaussian may be used as a trial order parameter in calculating the ideal trap-opening mechanism, while retaining the droplet state, commonly referred to as a shortcut to adiabaticity [71,72].