Quantum Droplets in Imbalanced Atomic Mixtures

Quantum droplets are a quantum analogue to classical fluid droplets in that they are self-bound and display liquid-like properties -- such as incompressibility and surface tension -- though their stability is the result of quantum fluctuations. One of the major systems for observing quantum droplets is two-component Bose gases. Two-component droplets are typically considered to be balanced, having a fixed ratio between the densities of the two component. This work goes beyond the fixed density ratio by investigating spherical droplets in imbalanced mixtures. With increasing imbalance, the droplet is able to lower its energy up to a limit, at which point the droplet becomes saturated with the atoms of the majority component and any further atoms added to this component cannot bind to the droplet. Analysing the breathing mode dynamics of imbalanced droplets indicates that the droplet can emit particles, as in balanced mixtures, but the imbalance leads to an intricate superposition of multiple simultaneously decaying collective oscillations.


I. INTRODUCTION
Quantum gases are highly controllable systems due in part to their diluteness, yielding a highly versatile platform for exploring quantum many-body physics [1].For a system of bosons, a mean-field model for quantum gases is Gross-Pitaevskii (GP) theory [2,3], which captures much of the physics of these systems, whilst neglecting beyond-mean-field effects such as quantum fluctuations [4,5].For a single-component homogeneous Bose gas stability is governed by the atomic interactions, with the gas becoming unstable to collapse for attractive interactions.This has been demonstrated experimentally by using a Feshbach resonance to tune the interactions from repulsive to attractive [6,7].For a two-component system, the stability depends on both the inter-and intraspecies contact interactions, which can likewise be tuned via Feshbach resonances [8][9][10].For repulsive intraspecies interactions -to ensure stability of each individual component -and attractive interspecies interactions, meanfield theory once more predicts an unstable gas.However, this attractive collapse can be stabilised at high densities by quantum fluctuations [11] forming a self-bound quantum droplet; this has inspired work in so called Lee-Huang-Yang (LHY) fluids [12,13], named for the correction describing quantum fluctuations to first-order.By tuning mean-field interactions between the components to vanish, the interactions of LHY fluids are described by quantum fluctuations alone [12,13].
In free space quantum droplets exist in equilibrium with the vacuum and form a seemingly counter-intuitive dilute liquid-like state [11,14].This adds to the variety of properties that quantum gases can be used to investigate, as the majority of experiments are inherently within the gas phase -a property widely exploited in time-of-flight imaging [15] -whereas quantum droplets open up the field to explore the properties of quantum liquids, such as surface tension [11,16,17] and incompressibility [18,19], within controllable experimental conditions.
A defining property of three-dimensional (3D) twocomponent quantum droplets is self-evaporation [11].In certain regimes of the droplet's phase diagram, excitations will cause the droplet to shed atoms in order to relax to a system of lower energy.This occurs when the energies of these excitations exceed −µ, the particle emission threshold, where µ is the chemical potential of the droplet.Furthermore, the lowest energy monopole mode -the breathing mode -decays across the largest proportion of the droplet's phase diagram relative to other collective modes [11,35].Self-evaporation is a nonintuitive and remarkable property that is not exhibited by dipolar quantum droplets [18], or for the breathing mode of a 1D two-component droplet [36], and has not been experimentally observed, but it is crucial in understanding the dynamics of these objects.
The formation of a two-component quantum droplet also has an interesting property: density balancing.A key result from the pioneering work of Petrov is the en-ergetic favourability for the two component densities to maintain a fixed ratio n 2 /n 1 = const.[11] where n i is the number density of the ith component.By pairing this assumption with negligible spin modes [11] -i.e., assuming only in-phase density oscillations -the mixture can be modelled via a single macroscopic wavefunction.The majority of the literature has focused on such balanced droplets, with theoretical studies of imbalanced systems limited to dipolar mixtures [37,38] and lowdimensions [39,40].This work is a systematic study of the ground states and breathing modes of 3D spherical imbalanced quantum droplets in homonuclear mixtures.Adding atoms to one component of the mixture yields a lower energy configuration than the associated balanced droplet.This forms a droplet with a density imbalance in the core, and this imbalance can be increased to a limit at which point any further atoms cannot bind to the droplet.By investigating imbalanced droplets in the free space limit, this work explores how the density profiles, chemical potentials and breathing mode dynamics of the droplet are modified by the presence of a population imbalance.
This work begins with a discussion of the underlying theory in Section II.This theory is then applied in Sections III and IV to first isolate imbalanced droplet ground states and then to propagate these states in time, subject to an initial perturbation, to analyse the droplet breathing modes.Finally the main results and potential future research avenues are discussed in Section V.

II. THE MODEL
A zero-temperature mixture of two weakly-interacting, dilute Bose gases can be described by the energy functional [11,41] in which m i are the atomic masses of the i-th component, g ii and g 12 are the effective intra-and interspecies interaction strengths respectively, and are related to the intraand interspecies scattering lengths by g ii = 4π 2 a ii /m i and g 12 = 2π 2 a 12 (1/m 1 + 1/m 2 ).The first two terms of Equation ( 1) are the kinetic energy contributions; the next three terms describe the two-body interactions within and between the components.The final term is the LHY correction which, to first-order, describes the effects of quantum fluctuations on the condensate [42].
For a Bose-Bose mixture, the LHY correction takes the form [11] where f is a function, defined in [16,43], with arguments z = m 1 /m 2 , u = g 2 12 /(g 11 g 22 ) and x = (g 22 n 2 )/(g 11 n 1 ).The function f (z, u, x) means that the LHY term takes a relatively simple form in homonuclear mixtures (m 2 = m 1 ) [11], but takes a far more complex form in heteronuclear mixtures (m 2 = m 1 ) [11,16].The function argument, u, reduces to one by assuming that the mixtures lies at the critical point of attractive stability, i.e., g 2 12 = g 11 g 22 =⇒ u = 1, removing the issue of complex contributions resulting from an unstable phonon mode [11,44,45].It should be noted that this approximation is made only in this derivation and does not imply any parameter choice in subsequent analyses.Crucial to the formation of quantum droplets, quantum fluctuations stabilise the attractive mixture against collapse, and furthermore are ubiquitous in nature though often play a limited role in the physics of many quantum gas experiments as discussed in Section I.
The two component densities can each be related to a macroscopic wavefunction or order parameter, , can then be used to derive the coupled extended GP equations, ( The computational complexity of these LHY terms can be minimised by assuming a homonuclear mixture, giving the equal-mass coupled extended GP equations [11] i ∂Ψ Finally, the dimensional scalings r = ξr, t = τ t and Ψi result in the dimensionless, equal-mass coupled extended GP equations, is the equilibrium density of component-1 for the balanced mixture.The expression of the equilibrium density is calculated in a homogeneous infinite system under the criterion of a vanishing pressure -i.e., the droplet in equilibrium with the vacuum -and takes the form [11] The density scalings ρ i correspond to rescaled normalisation constants, Ñi = N i /(ρ i ξ 3 ), in which N i is the normalisation constant of the i-th component wavefunction.In subsequent sections the dimensionless forms of N i , Ψ i , r and t are used, with tildes omitted for clarity.By assuming a constant density ratio, n 2 /n 1 = g 11 /g 22 , the two component wavefunctions can be expressed in terms of a single wavefunction, neglecting any out-of-phase motion between the components [11,46,47].Equations (2) reduce to a single equation, with the system described by a single parameter, an effective atom number, Ñ , given by [11] in which N here is the total atom number N = N 1 + N 2 .Within this work, balanced and imbalanced droplets are both modelled by Equations ( 3), though it should be noted that, for a balanced droplet, the dimensionless parameters (N 1 , N 2 , α, β, η) can be recast to Ñ .In the density-locked model, a given set of interaction strengths, g ii and g 12 , correspond to a fixed population number ratio, N 2 /N 1 = g 11 /g 22 , imposing that the atoms can only reside bound to the droplet.By breaking the assumption of density-locking, it is possible to imbalance population numbers such that N 2 /N 1 = g 11 /g 22 .The next section explores the effect of this imbalance on the density structure and energy of 3D spherical droplet ground states.

III. GROUND STATES
This work considers spherically-symmetric droplets.Density is hence assumed to be a function of radius only, reducing the computational problem to an effective 1D system with the kinetic term becoming ]/r.Further assumptions within this work include: a homonuclear mixture and balanced intraspecies scattering lengths (a 11 = a 22 =⇒ β = 1), to simplify the problem.Thus, the only differences between the components can arise from an imposed atom number imbalance of N 1 = N 2 + δN 1 .

A. Flat-top Droplet Limit
A simplified ansatz for the density of large droplets, as used in Ref. [11], is a flat-top, step function in which kinetic energy contributions are neglected.The twocomponents of the imbalanced droplet are modelled as majority and minority components with N + δN and N atoms, respectively.Substituting this ansatz into the two-component energy functional gives the equilibrium energy of the flat-top droplet, E eqbm , as a function of the imbalance, δN (see the Appendix for further details, including a schematic of the ansatz).The results of this energy calculation are given in Figure 1 with the total equilibrium energy given in the main plot, and the inset showing the equilibrium energy per particle.
Firstly, the main plot of Figure 1 shows that there is a minimum at δN = 0 in the total energy, indicating the droplet can lower its energy by absorbing particles into one component.However, the equilibrium energy per particle given in the inset of Figure 1 exhibits a minimum at δN = 0 demonstrating that whilst the total energy of the droplet can be lowered by single-component absorption, there is a corresponding increase in the energy per particle resulting in a less stably bound droplet.Secondly, the Appendix shows that the minimum in the equilibrium energy occurs at limit.The dependence of this quantity on the scattering lengths is consistent with the prediction of Ref. [11] at linear order [48].For the parameters of Figure 1 δn ≈ 0.458, as highlighted by the orange vertical line.
At this limit the droplet becomes saturated with the majority component.Hence if further majority component atoms are available, they cannot be absorbed into the droplet and will form an unbound gas around the droplet (as described in Ref. [11]).There are thus three droplet states: (1) balanced droplets; (2) bound, imbalanced droplets; and (3) saturated, imbalanced droplets (which can be immersed in a halo of unbound atoms).
The decrease in energy of an imbalanced droplet is consistent with the predictions of Ref. [11] and stems from the basic argument that droplets always seek to lower their energy by absorbing atoms.For example, in the density-locked model, atoms can only be added to both components to ensure that N 2 /N 1 = g 22 /g 11 is preserved, and by adding atoms to both components, a lower energy state is recovered.One way to conceptualise this is due to droplets having negative chemical potentials, meaning that it is always energetically favourable to absorb more particles into both components.However, the subtlety in the energy calculation given in the Appendix is that the energy of a droplet can also be lowered by absorbing particles into one component if there is an excess of this component available, up to the saturation limit discussed above.
By using the simple flat-top density ansatz, key insights are gained into effects of a population imbalance on the density structure and energy of a two-component droplet.However, for smaller droplet sizes kinetic energy cannot be neglected; thus, to extend this analysis to general, spherical droplets, the coupled GP equations are solved numerically.

B. Numerical Solutions
To find ground state solutions, the coupled GP equations are propagated numerically in imaginary time until the energy of the mixture is deemed adequately converged.The numerical scheme is a 4 th -order Runge-Kutta method, with centred finite-difference methods for the spatial derivatives.Neumann boundary conditions ( ∂Ψ i /∂r = 0) are applied at the centre of the computational box as the density in the droplet core is approximately constant.The boundary at r = L r , where L r is the radial computational box size, has either Neumann or Dirichlet [Ψ i (r = L r ) = 0] boundary conditions, which are discussed further below with reference to Figure 2.
By enforcing the Dirichlet boundary condition at r = L r , Figure 2(a), (b) and (c) show the three different ground state solutions (discussed above) by varying δN 1 , i.e., by varying the imbalance of the mixture.Figure 2(a) presents an example of a balanced droplet (N 1 = N 2 =⇒ δN 1 = 0) with an inset of the density difference ∆n(r) = n 1 (r) − n 2 (r), showing that the two component densities are identical.By increasing the population imbalance to δN 1 ≈ 450, as shown in Figure 2(b), the two component densities start to split within the droplet core, with both an increase and decrease in the majority and minority components' central densities, respectively.Hence, the two components are no longer identical but notably the density differences in the inset shows that this imbalanced droplet has no unbound atoms around the droplet.This implies that this state is a bound, imbalanced droplet.Driving the imbalance higher to δN 1 ≈ 5631, Figure 2(c) demonstrates a more pronounced density splitting within the droplet core, but crucially also exhibits a non-zero gas density outside of the droplet [highlighted in the inset of (c)].This ground state hence corresponds to a saturated, imbalanced droplet at the centre of the box with an unbound gas cloud outside of the droplet core.
For a set droplet size -i.e., fixed α, η and N 2 -the only free parameter is δN 1 and thus Figure 2(d) shows two measures of the droplet ground states for varying δN 1 .The two measures are the central density difference, ∆n(r = 0) = n 1 (r = 0) − n 2 (r = 0), and the majority component density, n 1 , at a fixed radius outside of the droplet.The chosen fixed radius is r = 0.75L r , which is chosen to be relatively far from the droplet surface but not too close to r = L r , to avoid the density drop resulting from the Dirichlet boundary condition.These measures are used to illustrate the transition between bound, imbalanced droplets, and saturated, imbalanced droplets.Figure 2(d) shows that for increasing δN 1 from zero, there is an immediate density splitting within the droplet corresponding to the approximately linear increase in ∆n(r = 0).It should be noted that in this regime, n 1 (r = 0.75L r ) stays fixed at zero, meaning that this relatively small δN 1 regime corresponds to bound, imbalanced droplets.However, increasing δN 1 further eventually leads to the formation of a shoulder  .At δN1 = 0 the system is balanced, with the leftmost vertical grey line corresponding to (a).Increasing δN1 leads to an approximately linear increase in ∆n(r = 0), there are however no atoms outside of the droplet as n1(r = 0.75Lr), a measure of density outside of the droplet, is zero.This regime corresponds to a bound, imbalanced droplet with the central vertical grey line corresponding to (b).The linear increase in ∆n(r = 0) eventually plateaus to a regime in which the central densities in the droplet cannot imbalanced further and the imbalanced droplet saturates with the density outside of the droplet increasing, resulting in a non-zero density of unbound atoms.The inset shows a measure of the difference in central densities in the limit of a saturated imbalanced droplet, ∆n sat (r = 0), as a function of droplet size given by the effective atom number Ñ in Equation (4), using the parameters for the δN1 = 0 case.The orange, dashed, horizontal line corresponds to δn ≈ 0.458, with the droplet size used in the main plot highlighted by the grey, vertical line.(e) Neumann boundary conditions are applied at both boundaries of the computational box, with Lr = 1024.The main plot shows the difference in central densities, recreating the results from (d), showing that the behaviour is not a function of the boundary conditions.It likewise demonstrates that the saturated droplet reaches a lower energy, E[Ψ1, Ψ2], state than the balanced droplet.The inset shows an example imbalance of δN1 ≈ 22524, with varying box size Lr.In the limit of large a computational box any background gas will be effectively zero density, thus converging to a saturated droplet in free space.
in ∆n(r = 0) which is the saturation of the imbalanced droplet, i.e., the droplet is reaching a limit of how many majority component atoms can be absorbed.Beyond this shoulder ∆n(r = 0) then approaches a constant value [labelled here as ∆n sat (r = 0)].Once in the saturated droplet regime, there is a corresponding increase in n 1 (r = 0.75L r ), indicating that the excess atoms reside outside of the droplet.The phase diagram in Figure 2(d) corresponds to the three ground state densities in (a), (b) and (c) given by the vertical grey lines.Furthermore, the inset of Figure 2(d) illustrates how the central density splitting of the saturated droplet varies with droplet size, given by Ñ , the effective atom number of the balanced droplet in Equation (4) (from Ref. [11]).To vary droplet size, the fixed parameters are α and η, thus N 2 (and hence N 1 ) is varied for the different droplet sizes, with Ñ calculated for the balanced droplet (δN 1 = 0) parameters.The central density difference of the saturated droplet increases with droplet size and asymptotically approaches a constant value of δn ≈ 0.458 [given by the orange, dashed, horizontal line in the inset of (d)] in the limit of an large, flat-top droplet as discussed above.
Figure 2(c) and (d) show that in the limit of the sat-urated droplet, there exists the unbound gas.This work seeks to explore imbalanced droplets, in the absence of this unbound gas, to probe the physics of imbalancing the droplet in free space.Thus, rather than imposing a Dirichlet boundary condition at r = L r , instead the same Neumann boundary condition used at r = 0, is used at r = L r , also.This gives a numerical approximation to a free space system.Figure 2(e) uses equivalent parameters as in (d), with a slightly increased δN 1 range, and exhibits the same central density difference behaviour as with the Dirichlet boundary conditions at r = L r .From this point on, Neumann boundary conditions are applied at both boundaries.
Figure 2(e) also shows the energy of the mixture which uses the dimensionless form of Equation ( 1), by defining Ẽ = E/ where = [4/(3π 2 α)][ 2 /(mξ 2 )], in which tildes are subsequently neglected.The energy decreases with imbalance -as predicted by the analytic calculated presented in the Appendix -until again reaching the saturation limit.Note again that a saturated, imbalanced droplet cannot absorb further majority component atoms and thus the energy of the droplet does not vary with δN 1 .However, this is an effect of the large computational boxes used, in that the density of the unbound atoms is negligibly small such that once in this saturation limit the system is essentially identical with changing δN 1 .In smaller computational boxes -or equivalently, exploring much larger imbalancesthe density of the background gas becomes significant [e.g. Figure 2(c)], which is explored in the inset of (e).By fixing δN 1 ≈ 22524, i.e., an imbalanced droplet in the saturation limit, the size of the computational box is varied.In small boxes, the density of the unbound gas is non-zero and thus has a significant, positive energy contribution.The existence of the background gas likewise modifies the internal structure of the droplet as can be seen by the increased value of ∆n(r = 0).In the limit of large L r (or equivalently small 1/L r ) the energy contribution of the background approaches zero as the density of the gas approaches zero.Thus in the limit of a large box, all of the ground states converge to the saturated droplet in free space, as can be seen from the convergence of both the energy and ∆n(r = 0) in the inset of Figure 2(e).
Further evidence of the energy saturation is given by the two component chemical potentials in Figure 3(a).For increasing δN 1 the majority component chemical potential (the dark orange curve) converges towards zero, implying there is no energetic favourability for the droplet to absorb more majority component atoms.Additionally, if more atoms of the minority component were made available to the droplet, then these atoms will likewise be absorbed into the droplet, reducing the energy of the droplet even further.
In summary, droplets will always seek to absorb atoms to reduce the system energy.If there is an imbalance of atoms from N 2 /N 1 = g 22 /g 11 then the droplet will absorb more of the majority component but this reaches a limit in which the droplet is saturated with the majority component atoms.It should be noted that in the case of a non-zero density unbound gas, this limiting behaviour will be modified by the positive energy contribution of the majority component density tails, similar to submerged droplets observed in 1D imbalanced droplets [39,40].These results have been presented both from an analytic energy calculation, and from numerical simulations, the latter of which is to be used next to explore the dynamical stability of these imbalanced droplets across their parameter space.

IV. BREATHING MODE
Self-evaporation is a key property of 3D twocomponent quantum droplets [11].Recalling that the density-locked droplet can be described by a single parameter, Ñ , there are three main regimes of interest for the droplet collective modes: 1) all modes exceed the particle emission threshold and hence are evaporated ( Ñ 94.2); 2) the monopole mode evaporates but other non-zero angular momentum modes are stable (94.2 Ñ 933.7); 3) the monopole is stabilised ( Ñ 933.7) [11,35].This first phase is the principal argument behind the self-evaporative property of a twocomponent droplet, since a finite-temperature droplet is inherently excited by contributions from, for example, the non-condensate components.Thus, these excitations are evaporated and the droplet relaxes to a lower energy configuration, akin to self-cooling [11].
In this work, the system is restricted to spherically symmetric droplets, meaning that the only observable mode is the breathing mode.The restriction of spherical symmetry is applied to reduce computational cost.This is due to particle shedding of self-evaporative droplets resulting in reflections from the computational box boundaries which become a substantial issue for long-time dynamics.To avoid this issue, very large box sizes -L r = 8192, approximately 500 times the droplet sizes considered here -are used to observe the dynamics of the droplets without interference from reflected particles.If simulating a general 3D droplet, these box sizes would quickly become infeasible.Note that the use of large computational boxes is crucial to the focus of this paper, namely imbalanced droplets in free space, i.e., without the external effects of unbound atoms.By excluding any non-zero angular momentum modes, the spherical symmetry yields two regimes: a decaying and a stable breathing mode, thus this work analyses the stability of breathing modes in the presence of an imbalance.
The breathing mode frequencies of the balanced, density-locked droplet are a function of the single parameter, Ñ [11].However, for an imbalanced system, a further degree of freedom is introduced, the size of the imbalance.As shown in the previous section, the central density differences, chemical potentials and energies, as functions of imbalance, reach the saturation limit and any further imbalance has no significant effect on the internal structure of the droplet.In order to observe these collective modes, the ground state solutions found via imaginary time propagation, shown in Section III, are then evolved in real time via Equations (3).To trigger a breathing mode in the droplet, a harmonic phase is imprinted on the ground state wavefunction, i.e., e i r 2 where is small (here = 10 −5 ) [4,49].This phase is always imprinted onto the minority component, though breathing modes are an in phase oscillation so the subsequent dynamics are not dependent on this choice of component.

A. Self-Evaporative Regime
The first modes to consider are within the selfevaporative regime ( Ñ 933.7).Breathing modes of spherical, balanced droplets have been quite extensively investigated [11,35,50].Dynamically, within this regime, the droplet begins to oscillate at a frequency exceeding the particle emission threshold, −µ, and rapidly decays with the decay rate asymptotically tending to zero and the oscillation frequency tending to the particle emission threshold [50].The initial rapid decay is due to a high dissipation of energy through particle emission, though in the long-time limit this corresponds to limited particle emissions at the energy of the chemical potential.This asymptotically decaying behaviour is likewise recovered in the imbalanced case though with three modes instead of one due to the imbalance yielding two chemical potentials.Figure 3(a) describes the three super-imposed modes: the early-time, rapidly decaying mode -see purple dashed lines in Figure 3(a) denoting the extracted equal frequencies of the two components -that varies marginally with imbalance, which can be considered as the intrinsic droplet breathing mode; at late times there are two further modes replacing the early-time mode.The late-time modes arise from the splitting of the chemical potentials which diverge with increasing imbalance before again reaching the saturated droplet [see orange lines in Figure 3(a)].
To visualise the single early mode and two late modes, Figure 3(b) and (c) represent a measure of the droplet central density, ni (t) = n i (r = 0, t) − n i (r = 0) t , with insets showing the associated power spectra |F [n i ]| 2 in which F [•] denotes the power spectrum rescaled by the mean, and all negative frequencies set to zero purely for data visualisation.The early-time mode is a highamplitude, rapidly decaying mode that is given in Figure 3(b), with an inset of the associated power spectrum highlighting the frequency (vertical dashed line) corresponding to the purple dashed lines in Figure 3(a).As this mode decays, it is then replaced by the two latetime modes corresponding to the split chemical potentials.The late-time dynamics is thus a superposition of a higher and lower frequency mode, as can be seen in Figure 3(c), once more with the associated power spec- 4257, N2 ≈ 17027, α ≈ 0.00657 and η ≈ −1.11.At δN1 = 0, the chemical potentials are equal but diverge as the imbalance is increased eventually reaching a saturation point where the chemical potential of the majority component is approximately zero.The minority component chemical potential likewise saturates, to an increased value.(b) A highlighted simulation of δN1 ≈ 1447 [corresponding to the bold vertical line in (a)] at early times, with an inset of the associated power spectrum.This indicates that at early times, there is a single rapidly decaying mode.(c) The same δN1 ≈ 1447 simulation at late times in which there is there is a superposition of two modes corresponding to the two chemical potentials, given by the vertical dashed lines in the inset power spectrum of (c).
tra showing peaks at the two chemical potentials (given by the vertical dashed lines).It should be noted that evidence of these late-time modes can even be seen at early times, such as the shorter peak in the inset of Figure 3(b), which roughly corresponds to the chemical potential of the majority component.
In summary, at early times the droplet oscillates with an unstable, high-amplitude mode that decays rapidly due to energy dissipation from particle emission.In the long-time limit however, the particle emission is considerably reduced, with particles only emitted at energies of the two chemical potentials.Hence, these late-time modes decay at much slower rates than the initial mode.This is analogous to the density-locked mixture [50], in which the mode frequency asymptotically converges to the particle emission threshold, with a vanishing decay rate.Physically this vanishing decay rate is the result of late-time emitted particles of energy ≈ −µ, which thus have a negligible effect on the kinetic energy of the droplet.Equivalently, for an imbalanced droplet, there are now two chemical potentials which each have associated particle emission and hence associated residual longlived breathing modes [50].

B. Non-Self-Evaporative Regime
The second set of breathing modes are within the nonself-evaporative regime ( Ñ 933.7).For a balanced droplet in this regime the breathing mode frequency is lower than the particle emission threshold resulting in a stable oscillation [11,50].In the self-evaporative regime, the dynamics of the imbalanced droplet is highly reminiscent of the balanced case but with further modes corresponding to the split chemical potentials, whereas the non-self-evaporative region of the droplet phase has some greater subtleties.
Figure 4(a) shows that, as in the self-evaporative case, the chemical potentials diverge with increasing imbalance until reaching the saturated, imbalanced droplet.The dynamics are again dominated by a high-amplitude mode in both components -given by the central purple dashed lines of Figure 4(a) -that is relatively constant with the changing imbalance.This mode can again be considered the intrinsic droplet breathing mode, as described in [51].In the balanced case, the stability of the breathing mode for droplets of this size is due to the particle emission threshold exceeding the mode frequency.With increasing imbalance however, as the chemical potentials split, the majority component chemical potential [i.e., the lower chemical potential branch in Figure 4(a)] eventually crosses over the frequency of the stable highamplitude mode, at which point this mode will begin to decay.This mode crossing implies that the imbalanced non-self-evaporative regime can instead be split into two regions: 1) a stable breathing mode; 2) a decaying breathing mode.This behaviour is highlighted in the upper inset of Figure 4(a), focusing on smaller imbalances.Decay of the high-amplitude, intrinsic mode occurs when the frequency exceeds the negative of the majority component chemical potential, validating the idea that stability of the mode is entirely dependent on the mode frequency exceeding the droplet's particle emission threshold.The lower left inset of Figure 4(a) shows that the mode is stable, for an imbalance of δN 1 ≈ 45, as the frequency lies beneath both chemical potentials.However, the lower right inset shows that this mode becomes unstable, for an imbalance of δN 1 ≈ 180, and decays as the frequency exceeds one of the chemical potential branches.The critical imbalance between stability and instability of this mode is δN 1 ≈ 124, for this specific droplet.Thus for a sufficiently small imbalance, it is possible to have an imbalanced droplet with a stable breathing mode.
Figure 4(b) and (c) focus on the early and late times of an unstable breathing mode in this regime.At early times, the initial high-amplitude mode dominates the system and oscillates at approximately the frequency of the balanced case.This is highlighted in the power spectrum shown in the inset of Figure 4(b) with the frequency given by the vertical dashed line.As found in the selfevaporative regime, the energy of the droplet is dissipated through atom shedding which then leads to the evaporation of the initial mode.The oscillations at late times -shown in Figure 4(c) -are instead dominated by two other modes, with frequencies corresponding to the chemical potentials of each component given by the two outer vertical dashed lines.There are still residual oscillations from the decaying initial mode which hence explains the interference seen in the oscillations of Figure 4(c) and the central peak shown in the power spectrum.Thus, in both the self-evaporative and non-selfevaporative regimes, providing there is decay of the initial high-amplitude mode, there are two main regions of the dynamics: early times -where the dynamics is dominated principally by a high-amplitude, intrinsic mode that is related to the balanced droplet; late timeswhere there is a superposition of two modes corresponding to each chemical potential.All of these modes are decaying but over different timescales due to the rate of atom shedding.The initial mode decays relatively rapidly due to a high dissipation of energy from the particle emission.However, at late times the particle emission is considerably reduced, with slower emission of particles.This slower emission of particles, at energies of the two chemical potentials, is negligible relative to the droplet kinetic energy and hence the dynamics decay asymptotically [50].
Finally, it should be noted that the breathing modes of imbalanced droplets have a larger range of instability than balanced droplets, which have an unstable breathing mode in the self-evaporative regime only.This is consistent with the discussions in both Section III and the Appendix ,of how the energy, and energy per particle, varies with imbalance.The key point from these discussions is that increasing imbalance corresponds to an increasing energy per particle and thus a less stably bound droplet.Hence, it would be expected that a less stably bound object would be more susceptible to particle shedding when perturbed, as is shown in the main results of this section.
V. DISCUSSION AND CONCLUSIONS By solving the spherically symmetric, coupled generalised GP equations, this work has systematically investigated the spherical ground states and breathing modes of imbalanced droplets across the size of the droplet and imbalance, in free space.Droplets can lower their energy by absorption of atoms, either by absorbing a symmetric number of atoms into each component, under the condition N 2 /N 1 = g 22 /g 11 , or, as has been investigated here, if there exists an asymmetric quantity of atoms available, the droplet will absorb the available atoms forming a majority and minority component.This modifies the core structure of the droplet by splitting the two central component densities, yielding a bound, imbalanced droplet.If further majority component atoms become available then the imbalanced droplet saturates and thus any excess majority components will not bind to the droplet and reside outside of the droplet as an unbound gas.To investigate the density profiles, chemical potentials and breathing dynamics of these imbalanced droplets, without the external effects of the unbound gas, large computational boxes were used such that the background gas density is effectively zero.
One of the main experimental probes to justify the observation of a quantum droplet is measuring a constant width of the atomic cloud, after switching off all traps, i.e., the object is self-bound in free space [25,27,28].In this scenario any unbound gas would be lost, though the imbalanced droplet could still remain stable in free space as demonstrated in this work.The observation of the three breathing mode frequencies could be experimentally achieved in homonuclear mixtures, for which the number of atoms in each component can be tuned by radio-frequency pulses on the cloud [26], and there have been recent experimental observations of imbalanced heteronuclear mixtures [52].However these system will likewise be affected by three-and higher-body losses.
An immediate next question of this work is to explore how the application of trapping potentials effects imbalanced droplets.Isotropic harmonic traps could be applied to the spherically symmetric system considered here, to investigate the modification of both the ground states and breathing mode dynamics.The breathing mode dynamics of trapped, spherically symmetric droplets has been investigated in Ref. [51] and so there are natural extensions of these results for imbalanced droplets.Further ideas include investigating the spherically symmetric ground states and dynamics of a heteronuclear mixture, as the differing kinetic energy contributions of the two components may lead to novel dynamics.This is computationally challenging however, due to the form of the two-component LHY correction of a heteronuclear mixture [11,16].There are approximations for this term, though they are least accurate at the droplet centres [43].This presents an issue for accurate analysis of imbalanced ground state solutions and dynamics, especially for the superimposed, time-dependent collective excitations shown here.
With the continuing development of quantum droplet experiments -for example through the realisation of new mixtures [53]-it will become increasingly important to understand how these mixtures can minimise their energy, from self-evaporating excitations to modifying their internal structure.
The data presented in this paper are available [54].
The density profile of a droplet in the limit of large Ñ approaches a flat-top, step function, i.e., a constant density in the droplet core and a steep drop to zero density at the surface [11].This approximate form is a useful model for infinite sized droplets as it allows for kinetic energy contributions to be neglected.Figure 5 shows a schematic of the step function model of a large imbalanced droplet, used in this derivation.The minority component is modelled as having volume V = 4  3 πR 3 , whilst the majority component has volume V + δV = 4  3 π(R + δR) 3 .The step functions hence correspond to the central densities n and n + δn for the minority and majority components, respectively, with the two components normalised to N and N + δN in which δN ≥ 0, imposing the population imbalance.The dimensionless form of the energy functional given in Equation ( 1), with β ≡ a 22 /a 11 = 1, can be written in terms of these constant imbalanced densities giving, The schematic in Figure 5 states nothing about the sign of δR, i.e., the volume of the majority component, V + δV , can be larger or smaller than the minority component volume, V .However, it can be shown that the only physical system is that where δV = 0, which is demonstrated in the following calculation.
Beginning with the δV > 0 case, Equation (A.1) can be written in terms of the component volumes and population numbers as, Recall from Section II that δa = a 12 + √ a 11 a 22 , and that defining δa < 0 implies dominantly attractive interactions, which can be written as η < −1 in the dimensionless parameters used here.Thus, with δN small, the mean-field terms becomes dominantly positive whilst the LHY contributions are small due to α 1.Hence, the limit of δV → 0 + yields, ∂E/∂(δV ) ≥ 0.
The next case to check is δV < 0, for which Equation (A.1) takes the form, This expression is again minimised with respect to the difference in component volumes, Hence, the limit of δV → 0 − results in, ∂E/∂(δV ) ≤ 0. Thus, the conditions ∂E ∂(δV ) δV →0 + ≥ 0, and, ∂E ∂(δV ) δV →0 − ≤ 0, imply that δV = 0 is either a smooth minimum or a cusp.This is physically realistic in the limit of large droplets, as this implies that the radii of the two components are equal, i.e., no single-component atoms reside outside of the core.However, in the finite droplet limit this is not the case as kinetic energy contributions cannot be neglected, as is highlighted in the ground state profiles in Figure 2 (b) and (c).The imbalanced droplet hence appears to slowly approach the flat-top density profile, relative to the rate at which the large, balanced droplet approaches the flat-top density limit.
Taking the case of δV = 0, the energy can be written as, .
A factor of (1 + δN/2N ) can be factored out of this expression, allowing for the above equation to be written in powers of (N/V ), .
Section II discusses the key property that droplets exist in equilibrium with the vacuum, i.e., they have zero pressure, ∂E/∂V = 0.With N constant, this expression for the zero pressure droplet can be rewritten as ∂(E/2N )/∂(N/V ) = 0, yielding, giving an expression for the equilibrium density.This results in the equilibrium energy, given by, or equivalently as the equilibrium energy per particle, For fixed α, η and N , Equations A.2 and A.3 are plotted as functions of δN in Figure 1.The main plot in Figure 1 shows Equation (A.2) varying with δN .By differentiating Equation (A.2), the minimum is located at For the parameters of Figure 1 the minimum is at δN ≈ 703, which using the expression, corresponds to δn ≈ 0.458, i.e., the orange, dashed, horizontal line in the inset of Figure 2(d).

FIG. 1 .
FIG. 1. Energy of a flat-top density ansatz as a function of imbalance, δN , with parameters N2 ≈ 22524, α ≈ 0.0152, β = 1.0 and η = −1.2.The equilibrium energy, E eqbm , is shown in the main plot with the orange, vertical line indicating the location of the minimum at δN ≈ 703.The inset shows the equilibrium energy per particle, [E/(2N + δN )] eqbm , with the only minimum appearing at the origin.

)FIG. 2 .
FIG.2.Balanced and imbalanced ground state droplets, with fixed parameters: N2 ≈ 22524, α ≈ 0.0152, β = 1.0 and η = −1.2.(a), (b) and (c) Ground state density profiles of: a balanced droplet (δN1 = 0); a bound, imbalanced droplet (δN1 ≈ 450); and a saturated, imbalanced droplet with an unbound halo of majority component atoms (δN1 ≈ 5631), respectively.The inset shows the difference in component densities, ∆n(r).(d) The difference in central densities as a function of imbalanced majority component atoms, δN1, with box size of Lr = 128.At δN1 = 0 the system is balanced, with the leftmost vertical grey line corresponding to (a).Increasing δN1 leads to an approximately linear increase in ∆n(r = 0), there are however no atoms outside of the droplet as n1(r = 0.75Lr), a measure of density outside of the droplet, is zero.This regime corresponds to a bound, imbalanced droplet with the central vertical grey line corresponding to (b).The linear increase in ∆n(r = 0) eventually plateaus to a regime in which the central densities in the droplet cannot imbalanced further and the imbalanced droplet saturates with the density outside of the droplet increasing, resulting in a non-zero density of unbound atoms.The inset shows a measure of the difference in central densities in the limit of a saturated imbalanced droplet, ∆n sat (r = 0), as a function of droplet size given by the effective atom number Ñ in Equation (4), using the parameters for the δN1 = 0 case.The orange, dashed, horizontal line corresponds to δn ≈ 0.458, with the droplet size used in the main plot highlighted by the grey, vertical line.(e) Neumann boundary conditions are applied at both boundaries of the computational box, with Lr = 1024.The main plot shows the difference in central densities, recreating the results from (d), showing that the behaviour is not a function of the boundary conditions.It likewise demonstrates that the saturated droplet reaches a lower energy, E[Ψ1, Ψ2], state than the balanced droplet.The inset shows an example imbalance of δN1 ≈ 22524, with varying box size Lr.In the limit of large a computational box any background gas will be effectively zero density, thus converging to a saturated droplet in free space.

2 FIG. 3 .
FIG.3.Chemical potentials and breathing modes as a function of majority component imbalance, δN1, in the selfevaporative regime (equivalent to a balanced droplet, N1 = N2, with Ñ ≈ 649).(a) The chemical potentials (light orange -minority component, dark orange -majority component) and early-time breathing mode frequencies (light and dark purple dashed lines) as a function of imbalance.These results correspond to the range 0 ≤ δN1 4257, N2 ≈ 17027, α ≈ 0.00657 and η ≈ −1.11.At δN1 = 0, the chemical potentials are equal but diverge as the imbalance is increased eventually reaching a saturation point where the chemical potential of the majority component is approximately zero.The minority component chemical potential likewise saturates, to an increased value.(b) A highlighted simulation of δN1 ≈ 1447 [corresponding to the bold vertical line in (a)] at early times, with an inset of the associated power spectrum.This indicates that at early times, there is a single rapidly decaying mode.(c) The same δN1 ≈ 1447 simulation at late times in which there is there is a superposition of two modes corresponding to the two chemical potentials, given by the vertical dashed lines in the inset power spectrum of (c).

2 FIG. 4 .
FIG. 4. Chemical potentials and breathing modes as a function of majority component imbalance, δN1, in the nonself-evaporative regime (equivalent to a balanced droplet, N1 = N2, with Ñ ≈ 1502), using the same parameters as in Figure 2 with 0 ≤ δN1 5631.(a) The chemical potentials (light orange -minority component, dark orange -majority component) and early-time breathing mode frequencies (light and dark purple dashed lines) as a function of the imbalance.The top inset shows the majority component chemical potential crossing the breathing mode frequencies.This is shown in the lower insets in which the mode is either stable (left) or is decayed (right).(b) A highlighted simulation with δN1 ≈ 1014 [corresponding to the bold vertical line in (a)] at early times, with the associated power spectrum (inset).The early time dynamics are dominated by a single damped breathing mode.(c) The same δN1 ≈ 1014 simulation as in (b), at late times showing a superposition of three modes, with the lowest and highest energy modes corresponding to the majority and minority component chemical potentials, respectively.The third mode [the central peak in the inset of (c)] is the early-time mode still decaying.

FIG. 5 .
FIG.5.Flat-top density profiles as an ansatz for each component of an imbalanced droplet.The orange component is the higher-density, majority component and the purple component is the lower-density, minority component.