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

P. Stürmer<EMAIL_ADDRESS>Mathematical Physics and NanoLund,
Lund University, Box 118, 22100 Lund, Sweden  M. Nilsson Tengstrand
Mathematical Physics and NanoLund, Lund University, Box 118, 22100 Lund,
Sweden  R. Sachdeva  Mathematical Physics and NanoLund, Lund University, Box
118, 22100 Lund, Sweden  S. M. Reimann  Mathematical Physics and NanoLund,
Lund University, Box 118, 22100 Lund, Sweden

###### Abstract

In this work, we present the study of the stationary structures and the
breathing mode behavior of a two-dimensional 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 Gross-Pitaevskii 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 Bulgac (2002);
Petrov (2015), 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 Kadau  _et al._ (2016); Schmitt  _et al._ (2016);
Ferrier-Barbut  _et al._ (2016); Chomaz  _et al._ (2016); Wenzel  _et al._
(2017); Ferrier-Barbut  _et al._ (2018), where the effective interatomic
interactions may be adjusted such that self-bound states become stabilized by
quantum fluctuations Saito (2016); Wächtler and Santos (2016a); Baillie  _et
al._ (2016); Bisset  _et al._ (2016); Wächtler and Santos (2016b); Chomaz  _et
al._ (2016); Baillie and Blakie (2018). The experimental discovery of droplets
in dipolar gases was soon followed by their realization in binary Bose gases
Semeghini  _et al._ (2018); Cabrera  _et al._ (2018), implementing Petrov’s
original proposal Petrov (2015). The quantum fluctuation contributions to the
total energy density functional are effectively represented by the Lee-Huang-
Yang (LHY) terms Lee  _et al._ (1957), extending the usual mean field (MF)
Gross-Pitaevskii equation (see Pethick and Smith (2008); Pitaevskii and
Stringari (2016) 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 Petrov (2015). In a binary three-dimensional
condensate the LHY terms are positive definite Lee  _et al._ (1957) 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 Petrov and Astrakharchik (2016); Zin  _et al._ (2018); Ilg
_et al._ (2018).
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 Cidrim  _et al._ (2018); Lee  _et al._ (2018) or contact-
interacting systems Li  _et al._ (2018); Nilsson Tengstrand  _et al._ (2019);
Kartashov  _et al._ (2018, 2019); Examilioti and Kavoulakis (2020), supersolid
behavior in dipolar and spin-orbit coupled systems Ancilotto  _et al._ (2020);
Böttcher  _et al._ (2019); Chomaz  _et al._ (2019); Nilsson Tengstrand  _et
al._ (2020); Sachdeva  _et al._ (2020); Roccuzzo and Ancilotto (2019);
Roccuzzo  _et al._ (2020); Tanzi  _et al._ (2019a, b) or collective
excitations Astrakharchik and Malomed (2018); Cappellaro  _et al._ (2018);
Otajonov  _et al._ (2019); Tylutki  _et al._ (2020); Hu and Liu (2020). For
trapped gases, the elementary excitation modes have been of fundamental
interest ever since the celebrated experimental realizations of atomic Bose-
Einstein condensates Anderson  _et al._ (1995); Bradley  _et al._ (1995);
Davis  _et al._ (1995): Examples range from a hydrodynamic description
Stringari (1996) to the Bogoliubov-de Gennes equations Edwards  _et al._
(1996); Pethick and Smith (2008); Pitaevskii and Stringari (2016) and to
solving the Euler-Lagrange equations Pérez-García  _et al._ (1997, 1996);
Salasnich (2000); Castin and Dum (1996). Breathing mode oscillations were
discussed in terms of the underlying symmetry properties of a Bose gas in two
dimensions Pitaevskii and Rosch (1997), and semiclassical methods have been
probed to describe low-lying excitations in the limit of large boson numbers
Sinha (1997). 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 Pitaevskii and Rosch (1997); Sinha (1997);
Svidzinsky and Fetter (1998); Linn and Fetter (2000); Mizushima  _et al._
(2004); Watanabe (2006); Simula  _et al._ (2008); Simula (2013); Teles  _et
al._ (2013a, b). More recently and in the context of Bose droplets, in Ref.
Jørgensen  _et al._ (2018) 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 $\mu\leqslant
0$, and are generally only observable if the associated excitation energy is
lower than $|\mu|$. However, this requirement can stay unfulfilled for self-
bound systems Petrov (2015) depending on the system’s dimensionality, if the
excitation energy exceeds $|\mu|$, see Refs. Tylutki  _et al._ (2020, 2020);
Cappellaro  _et al._ (2018).
Recent works on collective excitations in self-bound systems have mainly
focused on one-dimensional droplets Astrakharchik and Malomed (2018); Petrov
and Astrakharchik (2016) or on the cross-over between the soliton and droplet
regime Astrakharchik and Malomed (2018). Besides solving the Bogoliubov-de
Gennes equations in one Tylutki  _et al._ (2020) and in three dimensions
Petrov (2015); Hu and Liu (2020), a Gaussian trial order parameter is often
utilized to describe the properties of the droplet state Astrakharchik and
Malomed (2018); Cappellaro  _et al._ (2018); Hu and Liu (2020). 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 Astrakharchik and Malomed (2018), although
overestimating the results in three dimensions Hu and Liu (2020). In search
for a better approximation of the order parameter, a $\sqrt{1-\tanh(x)}$
function in the Thomas-Fermi limit has been utilized Examilioti and Kavoulakis
(2020), 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. Otajonov  _et al._ (2019); Baizakov  _et al._
(2011) 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 $\sqrt{1-\tanh(x)}$
function.

In continuation of our previous work Nilsson Tengstrand  _et al._ (2019), 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.1
and the breathing frequency in Sec. III.2. We continue in Sec. IV by studying
droplets carrying angular momentum in the form of a phase-imprinted 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

$\displaystyle\begin{split}i\frac{\partial\psi}{\partial
t}=&\Big{(}-\frac{1}{2}\nabla^{2}+|\psi|^{2}\ln|\psi|^{2}\Big{)}\psi~{},\end{split}$
(1)

see Ref. Petrov and Astrakharchik (2016). Here, $\psi=\psi(\mathbf{r},t)$ is
the order parameter with $\mathbf{r}=(x,y)$. As in Refs. Petrov and
Astrakharchik (2016); Nilsson Tengstrand  _et al._ (2019); Li  _et al._
(2018); Kartashov  _et al._ (2019); Sachdeva  _et al._ (2020), the scaling
invariances have been used to bring the system into dimensionless form, such
that it conveniently only depends on the norm $N=\int
d^{2}\mathbf{r}|\psi(\mathbf{r})|^{2}$. Given Eq. (1) the energy density
$\mathcal{E}[\psi]$ of the self-bound system is then given by

$\displaystyle\mathcal{E}[\psi]=\frac{1}{2}|\nabla\psi|^{2}+\frac{1}{2}|\psi|^{4}\ln{\frac{|\psi|^{2}}{\sqrt{e}}}.$
(2)

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 Astrakharchik and Malomed (2018); Cappellaro  _et al._ (2018); Hu and
Liu (2020), or a $\sqrt{1-\tanh{x}}$ ansatz in the Thomas-Fermi approximation
for two-dimensional droplets Examilioti and Kavoulakis (2020). While the
Gaussian ansatz for the breathing frequency yields small errors compared to
the numerical eGPe solution for one-dimensional problems Astrakharchik and
Malomed (2018), 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 Hu and Liu (2020). On the other hand, a $\sqrt{1-\tanh{x}}$
function accurately describes the spatial properties, but faces limitations to
describe the system for small norms Examilioti and Kavoulakis (2020) 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 Baizakov  _et al._ (2011); Otajonov  _et al._ (2019); Karlsson
(1992); Tsoy  _et al._ (2006). For the following calculation we thus use the
super-Gaussian as a trial order parameter. Assuming a circular droplet shape,
we write

$\displaystyle\psi(r,t)=A(t)\exp\left(ib(t)r^{2}-\frac{1}{2}\left(\frac{r}{R(t)}\right)^{2m}\right)~{}$
(3)

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 two-dimensional two-
component Bose gas for some generalized coordinate
$q_{i}=\\{A(t),b(t),R(t)\\}$, emerges from the respective Euler-Lagrange
equations

$\displaystyle\frac{\partial\mathfrak{L}}{\partial
q_{i}}=\frac{d}{dt}\frac{\partial\mathfrak{L}}{\partial\dot{q_{i}}},$ (4)

where the Lagrangian $\mathfrak{L}$ is calculated by averaging the Lagrangian
density $\mathcal{L}$ over the full space

$\displaystyle\mathfrak{L}=\langle\mathcal{L}\rangle=\int
d^{2}\mathbf{r}\,\mathcal{L}\,.$ (5)

The Lagrangian density $\mathcal{L}$ for the two-dimensional droplet system
results from the energy density for a self-bound system in Eq. (2) as

$\displaystyle\begin{split}\mathcal{L}[\psi,\psi^{*}]=\frac{i}{2}\left(\psi^{*}\frac{\partial\psi}{\partial
t}-\psi\frac{\partial\psi^{*}}{\partial
t}\right)-\mathcal{E}[\psi].\end{split}$ (6)

Normalizing the wave function gives the number of particles in the condensate,

$\displaystyle N=\int
d^{2}\mathbf{r}\,|\psi|^{2}=\frac{\pi}{m}\Gamma\left(\frac{1}{m}\right)A^{2}R^{2},$
(7)

where $\Gamma(m)$ is the gamma function. We continue to write
$\Gamma\left(\frac{1}{m}\right)=\Gamma_{1}$,
$\Gamma\left(\frac{2}{m}\right)=\Gamma_{2}$ and
$\Gamma\left(\frac{2}{m}\right)/\Gamma\left(\frac{1}{m}\right)=\Gamma$. 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

$\displaystyle\begin{split}\frac{\mathfrak{L}}{N}=&-\dot{b}R^{2}\Gamma-2b^{2}R^{2}\Gamma-\frac{m^{2}}{2R^{2}\Gamma_{1}}-\\\
&\frac{Nm}{2\pi R^{2}}\frac{2^{-1/m}}{\Gamma_{1}}\ln\left(\frac{Nm}{\pi
R^{2}}\frac{1}{\Gamma_{1}}\right)+\frac{N}{4\pi
R^{2}}\frac{2^{-1/m}}{\Gamma_{1}}\left(m+1\right).\end{split}$ (8)

Thus, the Euler-Lagrange equations for the variational parameters $b(t)$ and
$R(t)$ become

$\displaystyle b$ $\displaystyle=\frac{1}{2}\frac{\dot{R}}{R},$ (9)

and

$\displaystyle\begin{split}\ddot{R}=&\frac{1}{\Gamma_{2}}\Bigg{[}\frac{m^{2}}{R^{3}}+\frac{N}{2\pi
R^{3}}2^{-1/m}(m-1)+\\\ &\frac{Nm}{\pi R^{3}}2^{-1/m}\ln\left(\frac{Nm}{\pi
R^{2}}\frac{1}{\Gamma_{1}}\right)\Bigg{]}=-\frac{dU_{\textrm{eff}}}{dR}.\end{split}$
(10)

Let us now first determine $m$ as a function of $N$. We extremize the
Lagrangian $\mathcal{L}$ with respect to $m$, which gives the transcendental
equation with the analytical exact value $m_{\textrm{exact}}$

$\displaystyle
0=\frac{m_{\textrm{exact}}}{2}(\ln{2}-m_{\textrm{exact}})+\frac{N}{4\pi}2^{-1/m_{\textrm{exact}}}(2\ln{2}-1).{}$
(11)

For large $N$ we have an approximation $m_{\textrm{approx}}$

$\displaystyle m_{\textrm{approx}}\approx\sqrt{\frac{N}{2\pi}(\ln{4}-1)}.$
(12)

In Fig. 1 we plot $m$ according to Eqs. (11) and (12) as dashed and dotted-
dashed lines respectively, as well as the difference $\Delta
m=m_{\textrm{exact}}-m_{\textrm{approx}}$ (see inset). The difference $\Delta
m$ approaches $0$ as $N\rightarrow\infty$, so that we can use the approximate
expression in Eq. (12) for droplets with large $N$.

Figure 1: Optimal value for the exponent $m$ of the super-Gaussian ansatz with
increasing norm $N$. Solved numerically for Eq. (11) _(dashed line)_ and
approximated for large $N$ according to Eq. (12) _(dotted-dashed line)_. The
inset shows the difference $\Delta m=m_{\text{exact}}-m_{\text{approx}}$. With
increasing $N$ the error $\Delta m$ approaches zero, suggesting that the
analytical approximation is suitable for flat-top droplets.

Looking at the minimum of the effective potential $U_{\text{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

$\displaystyle
A^{2}|_{t=0}=n_{m}=\exp\left[\frac{1}{2m}(1-m)-m\pi\frac{2^{1/m}}{N}\right].$
(13)

For large norms the maximum density approaches $1/\sqrt{e}$ as expected for
this choice of scaling Nilsson Tengstrand  _et al._ (2019); Li  _et al._
(2018); Kartashov  _et al._ (2018); Sachdeva  _et al._ (2020). The mean radius
$\sqrt{\langle r^{2}\rangle}$ depends on the width of the cloud,

$\displaystyle\langle r^{2}\rangle=\Gamma R^{2},$ (14)

which simply reduces to the Gaussian case $\langle r^{2}\rangle=R^{2}$ for
$m=1$. Furthermore, the width of a stationary condensate is

$\displaystyle R_{0}^{2}=\frac{N}{\pi n_{m}}\frac{m}{\Gamma_{1}}.$ (15)

The low-lying excitation frequencies are found by assuming a single sinusoidal
oscillation in $R(t)$ with frequency $\omega_{0}$, such that the breathing
frequency $\omega_{0}$ is given by

$\displaystyle\frac{d^{2}U_{\textrm{eff}}}{dR^{2}}=\omega_{0}^{2}\,,$ (16)

and one obtains

$\displaystyle\omega^{2}_{0}=n_{m}^{2}\frac{\Gamma_{1}}{\Gamma}\frac{2\pi}{Nm}2^{-1/m}\,.$
(17)

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
$\mu$ as

$\displaystyle\mu=-\frac{n_{m}}{2}\left(\frac{m\pi}{N}+2^{-1/m}\right)$ (18)

With $N\rightarrow\infty$, $\mu$ then approaches $-1/(2\sqrt{e})$ as suggested
by Refs. Petrov and Astrakharchik (2016); Li  _et al._ (2018).

###  III.1 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\leq N\leq 1000$.
We start by comparing the mean radius $\sqrt{\langle r^{2}\rangle}$ 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_{\textrm{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 Cheiney  _et al._ (2018); Carr and Castin (2002);
Pérez-García  _et al._ (1998); Salasnich  _et al._ (2002).

Figure 2: Mean radius $\sqrt{\langle r^{2}\rangle}$ of the self-bound droplet
based on different approaches, as indicated by the legend. While the Gaussian
trial order parameter ($m=1$) is a good approximation for small $N$, it fails
to describe the mean radius accurately for larger norm.

As mentioned above, the maximum density $n_{m}$ approaches $1/\sqrt{e}$ and
$m_{\textrm{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 self-bound system. For small $N$, however, $m=1$ approximates
the super-Gaussian reasonably well, with only small deviations, while
$m_{\textrm{approx}}$ differs significantly as expected. All our numerical
values lie within a $2\%$ deviation directly below the super-Gaussian with
$m_{\textrm{exact}}$.

Figure 3: Maximum density $n_{m}$ of the droplet based on the variational and
numerical solutions, as indicated in the legend. As proposed in Refs. Petrov
and Astrakharchik (2016); Li  _et al._ (2018) $n_{m}$ approaches $1/\sqrt{e}$
with increasing $N$ for the super-Gaussian trial order parameter in its exact
and approximate form. With increasing norm $N$ the Gaussian ($m=1$) deviates
significantly. The numerical solutions of the eGPe, Eq. (1), are in good
agreement with the variational results for $m_{\mathrm{exact}}$ and
$m_{\mathrm{approx}}$.

Having analyzed the mean radius $\sqrt{\langle r^{2}\rangle}$ 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.

Figure 4: Density distributions $n(R)$ based on the Gaussian trial order
parameter (dashed), the super-Gaussian trial order parameter with
$m_{\text{exact}}$ (dotted-dashed) and a numerical solution of the eGPe (full
lines) for different $N$. With increasing $N$ the super-Gaussian density
distribution approaches the flat-top limit.

Finally, let us look at the chemical potential $\mu$, as specified it in Eq.
(18) above. Due to $\mu\leqslant 0$ in self-bound condensates, an underlying
self-evaporation process exists. Thus, in principle no other mode with energy
$E$ should be observable as long as $E/|\mu|\geqslant 1$. In a self-bound
droplet, this may not be always the case, depending on the system’s
dimensionality, as shown in Petrov (2015); Petrov and Astrakharchik (2016);
Cappellaro  _et al._ (2018). 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\sqrt{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$.

###  III.2 The breathing mode in a two-dimensional droplet

We now compare the analytical results for the breathing mode frequency
$\omega_{0}$ in Eq. (17) to the numerical results obtained by solving the
eGPe, Eq. (1), for $1\leq N\leq 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 $\sqrt{\langle r^{2}\rangle}$. The breathing frequency
$\omega_{0}$ is measured via Fourier-analysis of $\sqrt{\langle
r^{2}\rangle}$. We then calculate $\omega_{0}/|\mu|$ via the analytical
chemical potential in Eq. (18) and Fig. 5, and also compare to the numerical
solution.

Figure 5: Chemical potential $\mu$ 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. Figure 6: Breathing
frequency $\omega_{0}$ of the non-rotating droplet. For large $N$ the
numerical solutions of the eGPe, Eq. (1), coincide with the super-Gaussian
ansatz calculation, while for smaller $N$ the assumption of a single breathing
frequency in Eq. (16) no longer holds, due to the emergence of a beating
pattern.

Fig. 6 shows $\omega_{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
$\omega_{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 Hu and Liu (2020). However, such
deviations between analytical and numerical approaches are absent for the one-
dimensional self-bound system Astrakharchik and Malomed (2018).

Figure 7: $\omega_{0}/|\mu|$ for droplets based on the different approaches,
as specified in the legend. All three analytical solutions coincide for the
high norm regime, and are energetically favorable over the self-evaporation
for $N\gtrsim 36$.

With the chemical potential in Eq. (18) we calculate the ratio
$\omega_{0}/|\mu|$ shown in Fig. 7. Surprisingly, all three analytical
solutions overlap for high $N$, while the eGPe solution starts to deviate as
$\omega_{0}/|\mu|\rightarrow 1$. From the analytical results, we obtain that
the breathing mode is observable from $N\gtrsim 36$, while the numerical
solution of the eGPe suggests $N\gtrsim 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

$\displaystyle\psi_{0}(\mathbf{r})=Cr^{S}\exp{(-\alpha r^{2}+iS\theta)},$ (19)

where $C$ is a normalization constant, $r$ and $\theta$ radial and angular
coordinate, $\alpha$ some positive number, and $S$ is the vorticity Li  _et
al._ (2018). 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 $\omega_{0}$ via a Fourier-analysis of $\sqrt{\langle
r^{2}\rangle}$ and the chemical potential $\mu$. Here, we choose systems with
$L/N=1.0,\,2.0,\,3.0$ as shown in Fig. 8 for $N=850$.

Figure 8: Density distributions for self-bound droplets at $N=850$ with
angular momentum $L/N=1.0,$ $2.0,$ $3.0$. The size of the droplets and minimum
norm $N$ for stability follow the heuristically derived relation in Ref. Li
_et al._ (2018). The phase in cases with non-zero angular momentum is shown as
small insets in each figure. Figure 9: Frequency spectrum for droplets with
angular momentum as specified in the legend. With increasing $L/N$ the
breathing frequency $\omega_{0}$ decreases. The onset of each branch is given
by the real-time stability of the system. For the multiply quantized systems,
the stability requirement for the minimum required norm $N$ is equal to that
found in Li  _et al._ (2018). Figure 10: Availability of the breathing mode
for droplets with angular momentum as given in the figure. The breathing mode
exists if $\omega_{0}<|\mu|$ and is energetically favorable for a rotating
system as long as the stability requirement is fulfilled and occurs for
$L/N=0$ if $N\gtrsim 21$.

As stated before, the breathing mode will only be observable if the condition
$\omega_{0}<\mu$ is fulfilled. Due to the small oscillations of $\sqrt{\langle
r^{2}\rangle}$ when measuring $\omega_{0}$ we can calculate $\mu$ numerically
at any given point in the real-time evolution. Figs. 9 and 10 show the
breathing frequency $\omega_{0}$ and the ratio $\omega_{0}/|\mu|$ for the
systems as described above in the range $1\leq N\leq 1000$. Both quantities,
$\omega_{0}$ and $\omega_{0}/|\mu|$ follow the same pattern, where the droplet
without angular momentum corresponds to the highest lying branch followed by
the multiply-quantized vortices with lower $\omega_{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. Li  _et al._ (2018).
As the numerical solution is the same as in Fig. 7, the breathing mode
continues to be observable from $N\gtrsim 21$ while the analytical solution in
Fig. 7 becomes observable for $N\gtrsim 36$. This is in stark contrast to the
one-dimensional case where the breathing mode is observable at any norm $N$
Tylutki  _et al._ (2020) and the quasi one-dimensional case Cappellaro  _et
al._ (2018) where it is observable for most configurations.

##  V Conclusions

We have studied the breathing mode for symmetric two-dimensional 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 Otajonov  _et al._ (2019) and compared it to the widely used
Gaussian ansatz Astrakharchik and Malomed (2018); Cappellaro  _et al._ (2018);
Hu and Liu (2020) 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
$\omega_{0}/|\mu|$ barely exceeds $1$ when the analytical solution suggests
that $\omega_{0}/|\mu|>1$. Unlike two- and three-dimensional droplet Hu and
Liu (2020) systems, such discrepancy between the numerical and analytical
solutions is not obtained for the one-dimensional system Astrakharchik and
Malomed (2018). 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\gtrsim 36$, for
which $\omega_{0}/|\mu|<1$. This is similar to the quasi one-dimensional and
three-dimensional case Petrov (2015); Cappellaro  _et al._ (2018); Hu and Liu
(2020). However this behavior is different from that of a one-dimensional
system, for which the breathing mode is observable at any $N$ Astrakharchik
and Malomed (2018).
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 Nilsson Tengstrand  _et al._ (2019) we find that with increasing $L/N$
the breathing frequency decreases, and droplets with angular momentum have
$\omega_{0}/|\mu|<1$. Furthermore, the stability requirement in $N$ for
droplets to support a vortex with $L/N$ is the same as in Li  _et al._ (2018).
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 Xu  _et al._ (2020);
Guéry-Odelin  _et al._ (2019).

###### Acknowledgements.

Acknowledgements. Discussions with J. Bengtsson and A. Idini are gratefully
acknowledged. This work was financially supported by the Knut and Alice
Wallenberg Foundation, the Swedish Research Council and NanoLund.

## References

  * Bulgac (2002) A. Bulgac, Phys. Rev. Lett. 89, 050402 (2002).
  * Petrov (2015) D. S. Petrov, Phys. Rev. Lett. 115, 155302 (2015).
  * Kadau  _et al._ (2016) H. Kadau, M. Schmitt, M. Wenzel, C. Wink, T. Maier, I. Ferrier-Barbut, and T. Pfau, Nature 530, 194 (2016).
  * Schmitt  _et al._ (2016) M. Schmitt, M. Wenzel, F. Böttcher, I. Ferrier-Barbut, and T. Pfau, Nature 539, 259 (2016).
  * Ferrier-Barbut  _et al._ (2016) I. Ferrier-Barbut, H. Kadau, M. Schmitt, M. Wenzel, and T. Pfau, Phys. Rev. Lett. 116, 215301 (2016).
  * Chomaz  _et al._ (2016) L. Chomaz, S. Baier, D. Petter, M. J. Mark, F. Wächtler, L. Santos, and F. Ferlaino, Phys. Rev. X 6, 041039 (2016).
  * Wenzel  _et al._ (2017) M. Wenzel, F. Böttcher, T. Langen, I. Ferrier-Barbut, and T. Pfau, Phys. Rev. A 96, 053630 (2017).
  * Ferrier-Barbut  _et al._ (2018) I. Ferrier-Barbut, M. Wenzel, F. Böttcher, T. Langen, M. Isoard, S. Stringari, and T. Pfau, Phys. Rev. Lett. 120, 160402 (2018).
  * Saito (2016) H. Saito, Journal of the Physical Society of Japan 85, 053001 (2016).
  * Wächtler and Santos (2016a) F. Wächtler and L. Santos, Phys. Rev. A 93, 061603 (2016a).
  * Baillie  _et al._ (2016) D. Baillie, R. M. Wilson, R. N. Bisset, and P. B. Blakie, Phys. Rev. A 94, 021602 (2016).
  * Bisset  _et al._ (2016) R. N. Bisset, R. M. Wilson, D. Baillie, and P. B. Blakie, Phys. Rev. A 94, 033619 (2016).
  * Wächtler and Santos (2016b) F. Wächtler and L. Santos, Phys. Rev. A 94, 043618 (2016b).
  * Baillie and Blakie (2018) D. Baillie and P. B. Blakie, Phys. Rev. Lett. 121, 195301 (2018).
  * Semeghini  _et al._ (2018) G. Semeghini, G. Ferioli, L. Masi, C. Mazzinghi, L. Wolswijk, F. Minardi, M. Modugno, G. Modugno, M. Inguscio, and M. Fattori, Phys. Rev. Lett. 120, 235301 (2018).
  * Cabrera  _et al._ (2018) C. R. Cabrera, L. Tanzi, J. Sanz, B. Naylor, P. Thomas, P. Cheiney, and L. Tarruell, Science 359, 301 (2018).
  * Lee  _et al._ (1957) T. D. Lee, K. Huang, and C. N. Yang, Phys. Rev. 106, 1135 (1957).
  * Pethick and Smith (2008) C. J. Pethick and H. Smith,  _Bose–Einstein Condensation in Dilute Gases_, 2nd ed. (Cambridge University Press, 2008).
  * Pitaevskii and Stringari (2016) L. Pitaevskii and S. Stringari,  _Bose-Einstein condensation and superfluidity_, Vol. 164 (Oxford University Press, 2016).
  * Petrov and Astrakharchik (2016) D. S. Petrov and G. E. Astrakharchik, Phys. Rev. Lett. 117, 100401 (2016).
  * Zin  _et al._ (2018) P. Zin, M. Pylak, T. Wasak, M. Gajda, and Z. Idziaszek, Phys. Rev. A 98, 051603 (2018).
  * Ilg  _et al._ (2018) T. Ilg, J. Kumlin, L. Santos, D. S. Petrov, and H. P. Büchler, Phys. Rev. A 98, 051604 (2018).
  * Cidrim  _et al._ (2018) A. Cidrim, F. E. A. dos Santos, E. A. L. Henn, and T. Macrì, Phys. Rev. A 98, 023618 (2018).
  * Lee  _et al._ (2018) A.-C. Lee, D. Baillie, R. N. Bisset, and P. B. Blakie, Phys. Rev. A 98, 063620 (2018).
  * Li  _et al._ (2018) Y. Li, Z. Chen, Z. Luo, C. Huang, H. Tan, W. Pang, and B. A. Malomed, Phys. Rev. A 98, 063602 (2018).
  * Nilsson Tengstrand  _et al._ (2019) M. Nilsson Tengstrand, P. Stürmer, E. O. Karabulut, and S. M. Reimann, Phys. Rev. Lett. 123, 160405 (2019).
  * Kartashov  _et al._ (2018) Y. V. Kartashov, B. A. Malomed, L. Tarruell, and L. Torner, Phys. Rev. A 98, 013612 (2018).
  * Kartashov  _et al._ (2019) Y. V. Kartashov, B. A. Malomed, and L. Torner, Phys. Rev. Lett. 122, 193902 (2019).
  * Examilioti and Kavoulakis (2020) P. Examilioti and G. Kavoulakis, (2020), arXiv:2005.13874 .
  * Ancilotto  _et al._ (2020) F. Ancilotto, M. Barranco, M. Pi, and L. Reatto,  (2020), arXiv:2002.05934 .
  * Böttcher  _et al._ (2019) F. Böttcher, J.-N. Schmidt, M. Wenzel, J. Hertkorn, M. Guo, T. Langen, and T. Pfau, Phys. Rev. X 9, 011051 (2019).
  * Chomaz  _et al._ (2019) L. Chomaz, D. Petter, P. Ilzhöfer, G. Natale, A. Trautmann, C. Politi, G. Durastante, R. M. W. van Bijnen, A. Patscheider, M. Sohmen, M. J. Mark, and F. Ferlaino, Phys. Rev. X 9, 021012 (2019).
  * Nilsson Tengstrand  _et al._ (2020) M. Nilsson Tengstrand, D. Boholm, R. Sachdeva, J. Bengtsson, and S. M. Reimann, (2020), arXiv:2003.06408 .
  * Sachdeva  _et al._ (2020) R. Sachdeva, M. N. Tengstrand, and S. M. Reimann, Phys. Rev. A 102, 043304 (2020).
  * Roccuzzo and Ancilotto (2019) S. M. Roccuzzo and F. Ancilotto, Phys. Rev. A 99, 041601 (2019).
  * Roccuzzo  _et al._ (2020) S. M. Roccuzzo, A. Gallemí, A. Recati, and S. Stringari, Phys. Rev. Lett. 124, 045702 (2020).
  * Tanzi  _et al._ (2019a) L. Tanzi, E. Lucioni, F. Famà, J. Catani, A. Fioretti, C. Gabbanini, R. N. Bisset, L. Santos, and G. Modugno, Phys. Rev. Lett. 122, 130405 (2019a).
  * Tanzi  _et al._ (2019b) L. Tanzi, J. Maloberti, G. Biagioni, A. Fioretti, C. Gabbanini, and G. Modugno,  (2019b), arXiv:1912.01910 .
  * Astrakharchik and Malomed (2018) G. E. Astrakharchik and B. A. Malomed, Phys. Rev. A 98, 013631 (2018).
  * Cappellaro  _et al._ (2018) A. Cappellaro, T. Macrì, and L. Salasnich, Phys. Rev. A 97, 053623 (2018).
  * Otajonov  _et al._ (2019) S. R. Otajonov, E. N. Tsoy, and F. K. Abdullaev, Physics Letters A 383, 125980 (2019).
  * Tylutki  _et al._ (2020) M. Tylutki, G. E. Astrakharchik, B. A. Malomed, and D. S. Petrov, Phys. Rev. A 101, 051601 (2020).
  * Hu and Liu (2020) H. Hu and X.-J. Liu, Phys. Rev. A 102, 053303 (2020).
  * Anderson  _et al._ (1995) M. H. Anderson, J. R. Ensher, M. R. Matthews, C. E. Wieman, and E. A. Cornell, Science 269, 198 (1995).
  * Bradley  _et al._ (1995) C. C. Bradley, C. A. Sackett, J. J. Tollett, and R. G. Hulet, Phys. Rev. Lett. 75, 1687 (1995).
  * Davis  _et al._ (1995) K. B. Davis, M. O. Mewes, M. R. Andrews, N. J. van Druten, D. S. Durfee, D. M. Kurn, and W. Ketterle, Phys. Rev. Lett. 75, 3969 (1995).
  * Stringari (1996) S. Stringari, Phys. Rev. Lett. 77, 2360 (1996).
  * Edwards  _et al._ (1996) M. Edwards, P. A. Ruprecht, K. Burnett, R. J. Dodd, and C. W. Clark, Phys. Rev. Lett. 77, 1671 (1996).
  * Pérez-García  _et al._ (1997) V. M. Pérez-García, H. Michinel, J. I. Cirac, M. Lewenstein, and P. Zoller, Phys. Rev. A 56, 1424 (1997).
  * Pérez-García  _et al._ (1996) V. M. Pérez-García, H. Michinel, J. I. Cirac, M. Lewenstein, and P. Zoller, Phys. Rev. Lett. 77, 5320 (1996).
  * Salasnich (2000) L. Salasnich, International Journal of Modern Physics B 14, 1 (2000).
  * Castin and Dum (1996) Y. Castin and R. Dum, Phys. Rev. Lett. 77, 5315 (1996).
  * Pitaevskii and Rosch (1997) L. P. Pitaevskii and A. Rosch, Phys. Rev. A 55, R853 (1997).
  * Sinha (1997) S. Sinha, Phys. Rev. A 55, 4325 (1997).
  * Svidzinsky and Fetter (1998) A. A. Svidzinsky and A. L. Fetter, Phys. Rev. A 58, 3168 (1998).
  * Linn and Fetter (2000) M. Linn and A. L. Fetter, Phys. Rev. A 61, 063603 (2000).
  * Mizushima  _et al._ (2004) T. Mizushima, Y. Kawaguchi, K. Machida, T. Ohmi, T. Isoshima, and M. M. Salomaa, Phys. Rev. Lett. 92, 060407 (2004).
  * Watanabe (2006) G. Watanabe, Phys. Rev. A 73, 013616 (2006).
  * Simula  _et al._ (2008) T. P. Simula, T. Mizushima, and K. Machida, Phys. Rev. A 78, 053604 (2008).
  * Simula (2013) T. Simula, Phys. Rev. A 87, 023630 (2013).
  * Teles  _et al._ (2013a) R. P. Teles, F. E. A. dos Santos, M. A. Caracanhas, and V. S. Bagnato, Phys. Rev. A 87, 033622 (2013a).
  * Teles  _et al._ (2013b) R. P. Teles, V. S. Bagnato, and F. E. A. dos Santos, Phys. Rev. A 88, 053613 (2013b).
  * Jørgensen  _et al._ (2018) N. B. Jørgensen, G. M. Bruun, and J. J. Arlt, Phys. Rev. Lett. 121, 173403 (2018).
  * Baizakov  _et al._ (2011) B. B. Baizakov, A. Bouketir, A. Messikh, A. Benseghir, and B. A. Pumarov, International Journal of Modern Physics B 25, 2427 (2011).
  * Karlsson (1992) M. Karlsson, Phys. Rev. A 46, 2726 (1992).
  * Tsoy  _et al._ (2006) E. N. Tsoy, A. Ankiewicz, and N. Akhmediev, Phys. Rev. E 73, 036621 (2006).
  * Cheiney  _et al._ (2018) P. Cheiney, C. R. Cabrera, J. Sanz, B. Naylor, L. Tanzi, and L. Tarruell, Phys. Rev. Lett. 120, 135301 (2018).
  * Carr and Castin (2002) L. D. Carr and Y. Castin, Phys. Rev. A 66, 063602 (2002).
  * Pérez-García  _et al._ (1998) V. M. Pérez-García, H. Michinel, and H. Herrero, Phys. Rev. A 57, 3837 (1998).
  * Salasnich  _et al._ (2002) L. Salasnich, A. Parola, and L. Reatto, Phys. Rev. A 66, 043603 (2002).
  * Xu  _et al._ (2020) T.-N. Xu, J. Li, T. Busch, X. Chen, and T. Fogarty, Phys. Rev. Research 2, 023125 (2020).
  * Guéry-Odelin  _et al._ (2019) D. Guéry-Odelin, A. Ruschhaupt, A. Kiely, E. Torrontegui, S. Martínez-Garaot, and J. G. Muga, Rev. Mod. Phys. 91, 045001 (2019).