The cross-over from Townes solitons to droplets in a 2D Bose mixture

When two Bose-Einstein condensates -- labelled 1 and 2 -- overlap spatially, the equilibrium state of the system depends on the miscibility criterion for the two fluids. Here, we theoretically focus on the non-miscible regime in two spatial dimensions and explore the properties of the localized wave packet formed by the minority component 2 when immersed in an infinite bath formed by component 1. We address the zero-temperature regime and describe the two-fluid system by coupled classical field equations. We show that such a wave packet exists only for an atom number $N_2$ above a threshold value corresponding to the Townes soliton state. We identify the regimes where this localized state can be described by an effective single-field equation up to the droplet case, where component 2 behaves like an incompressible fluid. We study the near-equilibrium dynamics of the coupled fluids, which reveals specific parameter ranges for the existence of localized excitation modes.


I. INTRODUCTION
Mixtures of quantum fluids display novel phenomenology as compared to the single-component case, through the emergence of collective degrees-of-freedom [1,2].In this perspective, ultracold atomic gases have opened a new path for the investigation of such many-body problems, especially thanks to the precise control of interactions [3][4][5][6][7][8][9].In the case of miscible Bose mixtures with two spin components, the dispersion relations of density and spin linear excitations have been studied experimentally [10,11], as well as the nonlinear excitations known as magnetic solitons [12,13].The superfluid character of the spin degree of freedom was also demonstrated by observing undamped spin-dipole oscillations [14] and by moving a magnetic obstacle [15] in such a mixture.Moreover, it has been shown that a coherent coupling between the two components -with or without momentum transfer -can modify the bare dispersion relations in a controlled manner [11], induce spin-orbit coupling to produce a variety of quantum phases [16], and trigger dynamical instabilities [17].
Even when each isolated component is stable, an instability can occur in a binary mixture when the interaction between the two components is attractive and set above a certain threshold.Close to this threshold, the balance between the dominant energy contributions can lead to subtle phases.For instance, it was predicted in Ref. [18] that a mixture of repulsive quantum gases in three dimensions (3D) with finely-tuned mutual attraction may lead to self-bound states stabilized by quantum fluctuations.This novel state of matter, known as a quantum droplet, was realized experimentally in [19,20], and the link between these 3D droplets and 1D solitons was clarified in [21].Interestingly, the role of quantum fluctuations is known to be enhanced in low dimensions, resulting in quantum droplet states with properties distinct from the 3D case, both at equilibrium [22] and close to equilibrium [23].
In a different range of parameters, mutually repulsive fluids may experience phase separation, similarly to solutions of Helium 3 in Helium 4 at low temperatures [24] or, more prosaically, in a combination of oil and water.This demixing dynamical instability was also characterized using Bose mixtures [25,26].In the regime of strong population imbalance, it was realized early that the dynamics of immiscible mixtures may mimic that of a single closed equation for the minority component [27,28].Recently, this mapping was leveraged for the deterministic realization of Townes solitons in a 2D Bose mixture [29], by making judicious use of the almost coincidence of the various interaction strengths (see also [30] for another preparation protocol of the Townes soliton with matter waves and [31] for a general review).This approach was also put forward for the realization of other exotic nonlinear excitations, such as Peregrine solitons [32] or darkbright soliton trains [33].
In this work, we consider a two-component Bose gas and study localized wavepackets of one (minority) component surrounded by a 2D bath of atoms in the other component, see Fig. 1.In section II, we focus on the stationary states of the system.We show that when the atom number in the minority component increases above a threshold value N T , a cross-over transition occurs from a steady-state with a solitonic character (the Townes soliton) to a droplet-like state.Then, in section III, we explore the excitation spectrum of these localized states throughout the crossover.In particular, we show the existence of a given range of atom numbers (1.45 N/N T 3.5) where no localized excitation exist.Finally, we discuss in section IV some possible extensions of this work.

A. The two-component system
We consider the ground state of a 2D Bose mixture made of two components, labelled |1 and |2 with mass m 1,2 , and with short-ranged interactions.We use a classical field description of the two components with the two order parameters ψ 1,2 (r, t).Both intraspecies and interspecies interactions are assumed to be repulsive, g ij > 0 where i, j = 1, 2. In practice, a planar gas is obtained using a strong confinement along z.We consider here the quasi-2D situation where the thickness z is larger than the scattering lengths a ij .In such a situation, the classical field approach is valid when (n 3D a 3 ij ) 1/2 1, where n 3D is the maximal 3D density of the gas.This hypothesis of negligible beyond-mean-field contributions was well satisfied in the experiment reported in [29].
The evolution of ψ 1,2 is given by the set of coupled nonlinear Schrödinger equations (NLSEs) where The steady-state of the two-component system is obtained by solving numerically the set of equations (1) for ψ i (r, t) = e −iµit φ i (r), where µ 1,2 are the chemical potentials of each component.Away from the region where the minority component is localized, the density of the bath brings the energy scale µ 1 = g 11 n ∞ and the corresponding length scale (healing length) ξ = 1/ √ 2m When the bath energy scale µ 1 largely exceeds all energy scales governing the minority component dynamics, it is possible to derive a closed equation for the minority component only.This situation is realized when the density of the minority component n 2 is everywhere much smaller than the bath density n 1 , corresponding to a weak depletion of the bath.Under these conditions, atoms in the minority component get dressed by the bath, which induces an additional effective interaction between atoms in state |2 [34].Using Bogoliubov's approach in 2D and for m 2 m 1 , we show in appendix A that this mediated interaction is described by the potential where K 0 is the zeroth-order modified Bessel function of the first kind, with asymptotic behavior K 0 (r) ∼ e −r / √ r when r → +∞.This expression is analogous to the Yukawa potential which arises in a 3D geometry [1,35,36].Remarkably, these mediated interactions are always attractive -whatever the sign of g 12 -and their range is given by the bath healing length ξ.The same conclusion holds when the masses m 1,2 are comparable, although the mediated interaction has a more complicated structure in this case (see appendix A and [37,38]).When any characteristic length of the minority component is much larger than ξ, we can adopt a zero-range description for the mediated interactions.The effective coupling strength, obtained by summing the bare interaction coupling strength g and the mediated one, is then independent of the ratio m 2 /m 1 and is given by (see Appendix A) In this limit, the minority component time evolution can thus be approximated -at least for short times -by the following single-component NLSE (up to a constant energy contribution) We deduce that, in this weak-depletion regime, the existence of stationary localized states for component |2  1)) with g12 = 1.01 g with g ≡ √ g11g22, the green lines correspond to the weak-depletion model of Eq. ( 5) and the red lines refer to the effective single-component model of Eq. ( 7).All profiles are calculated for m1 = m2 and g11 = g22.
requires effective attractive interactions, i.e. g e < 0. This last condition is equivalent to the criterion for the immiscibility g 12 > g ≡ √ g 11 g 22 of the binary mixture.The weakly-depleted state is thus a precursor to an actual phase separation situation.
Eq. ( 4) is known to host the so-called Townes soliton [31,39,40].Mathematically, this soliton is the unique radially symmetric, real and node-less solution of the stationary version of Eq. ( 4).It exists only when the atom number in the minority component N equals the critical value N T ≡ G T /|g e | with G T 5.85.When this condition is satisfied, the Townes soliton can be formed with any size, a direct consequence of the scale invariance of Eq. ( 4) which does not feature any explicit length scale [41].Formally, the soliton size is set by an effective chemical potential µ < 0 associated to the stationary solution ψ 2 (r, t) = e −iµt φ 2 (r) of Eq. ( 4) with µ = µ 2 − g 12 n ∞ .The upper limit µ = 0 is obtained in the case of small depletion, i.e. n 2 n 1 ≈ n ∞ everywhere.The validity of the single-component description of Eq. ( 4) was demonstrated experimentally in Ref. [29] in this limit.

C. The weak-depletion limit
The scale invariance of Eq. ( 4) results from the assumption that the size of the minority component is large compared to the bath healing length ξ, which provides the range of the mediated interaction.The first-order correction to this assumption adds a weak nonlocal nonlinearity to Eq. ( 4), which explicitly breaks scale invari-ance and leads to the modified NLSE (see Appendix B) with β = −(g 12 /g 11 ) 2 /(4m 1 n ∞ ).This correction remains small in front of the two dominant terms as long as m 2 |β|n 2 1 and ξ s where ξ s = (g 12 /g 11 )/ 2m 1 |g e |n ∞ represents the interpenetration length -or "spin" healing length -of the immiscible mixture.
Eq. ( 5) was studied extensively in Ref. [42].The steady state associated with this equation results from the balance between three energetic contributions: (i) the kinetic energy per particle 1/m 2 2 , (ii) the main part of the interaction energy −N/(N T m 2 2 ) that balances kinetic energy irrespective of for N = N T , and (iii) the correction −βN/ 4 originating from the extra term in Eq. ( 5) in comparison with Eq. ( 4).For 0 . The extension of these states thus becomes very large when → 0, i.e.N → N + T , and their density profile approaches the Townes soliton solution of Eq.( 4).More quantitatively, it is shown in Ref. [42] that Localized states of Eq. ( 5) exist for any N > N T , by contrast to Eq. ( 4) that requires N = N T .When N is close to N T , we recover with this simple approach the phase diagram of Fig. 2(a) as well as the density profiles calculated numerically using Eq.(1) (Fig. 2, case α).For larger N/N T (Fig. 2, cases β and γ), the validity condition ξ s breaks down, the solution of Eq. ( 5) is no-tably different from the result derived from Eq.( 1), and is therefore not relevant for our problem.

D. The strong-depletion limit
When the atom number of the minority component N becomes much larger than N T , the central density of this component grows to n2 = n ∞ g 11 /g 22 and the bath is locally fully depleted, see Fig. 2(b), case γ, and Fig. 2(c).In this phase-separated regime, the pressures g ii n 2 i /2 in the two components are equal [1,2], and component |2 fills approximately uniformly a disk of radius R such that N πR 2 n2 .It thus forms an effective droplet similar to an incompressible fluid of fixed density, although this density is not intrinsic but imposed by the surrounding medium.In the general case, no simple approach is available to describe this regime and one should solve Eqs.(1).Nevertheless, we show in the next paragraph that, close to the SU(2)-symmetry point, the system's equilibrium state can still be described by a single-component equation.

E. The vicinity of SU(2) symmetry
We assume in this paragraph equal masses m 1 = m 2 ≡ m.The interactions are said to be SU(2) symmetric when all interaction parameters g ij are equal.Close to this point, i.e. when g 12 → g + , the stationary state of the mixture in the N > N T case can be determined by solving the single-component effective equation [29] (7) The data in Figs.2(b)-2(c) have been calculated for this regime of nearby coupling constants (g 12 = 1.01g).They show that the predictions derived from Eq. ( 7) accurately describe the equilibrium profiles, from the weak to the full depletion regime.
The comparison of the validity ranges for Eq. ( 7) and for the stationary version of Eq. ( 5) is instructive: both equations coincide when |ψ 2 | 2 n ∞ and g 12 ≈ g (Fig. 2(b), case α).Beyond this common validity domain, Eq. ( 5) allows one to address the case where g 12 differs notably from g, whereas Eq. ( 7) is valid for arbitrary depletions, hence arbitrary values of N/N T .Finally we note that one should refrain from using in the general case a time-dependent version of Eq. ( 7), which would be obtained by replacing the left-hand-side µψ 2 by i∂ t ψ 2 .Indeed, in the strong-depletion regime, there is no hierarchy between the time scales for the minority component and for the bath.Therefore, it is not possible to eliminate the bath dynamics and obtain a time-dependent equation for ψ 2 involving only a first-order time derivative.
Another remarkable situation which occurs in this SU(2) limit is the 1D dark-bright soliton, introduced by Manakov [43] and transposed by Busch and Anglin for a binary mixture of Bose gases [44].There, the majority component wavefunction exhibits a phase jump, akin to a dark soliton around which the minority component accumulates.The first observations of such solitons with cold atoms were reported in Ref. [45,46].In a 2D configuration, the equivalent situation would correspond to a vortex texture in the majority component with its core filled by the minority one.This "vortex-bright" soliton [47,48] is notably different from the stationary states explored in this article where each component exhibits a uniform phase.

III. EXCITATION SPECTRUM
We now turn to the dynamics of the localized component, restricting for simplicity to close-to-equilibrium phenomena.This problem goes by essence beyond the Townes soliton physics.Indeed it is known that the Townes soliton associated to Eq. ( 4) does not possess any localized mode with non-zero frequency [49].More dramatically, some arbitrarily small deformations, such as the multiplication by a phase factor e iαr 2 with α → 0, may lead to a collapse of the soliton.
For the two-component system of interest here, a natural approach to determine its excitation spectrum is provided by the Bogoliubov method applied to the coupled equations ( 1).This procedure is outlined in Appendix C and the results are indicated with full lines in Figs.3(a)-3(b).Note that we focus here on localized modes, i.e.Bogoliubov modal functions u i (r), v i (r) that decay exponentially to zero for r → ∞.We show in Appendix C that this constraint corresponds to a mode frequency ω smaller than the continuum set by |µ|, where µ is the effective chemical potential introduced above.We now discuss the results for the various relevant regimes.We consider simple approaches for the limiting cases of weak and strong depletion of the bath and also investigate the intermediate regime of self-evaporation, corresponding to the absence of localized excitations.For simplicity, we restrict in this section to the case of equal masses m 1 = m 2 ≡ m and equal intracoupling constants g 11 = g 22 = g.

A. Weak-depletion regime and breathing mode
In the weak-depletion regime, we expect from the results of the previous section that the dynamics of the minority component is well captured by Eq. ( 5), which takes into account the first-order corrections originating from the two-component nature of our system.Quite remarkably, the instability inherent to Eq. ( 4) does not occur for Eq. ( 5).Indeed, it was shown in Ref. [42] that the stationary solutions of Eq. ( 5) are dynamically stable.Moreover, these solutions can sustain a breathing mode, i.e. an oscillation of the system's overall size, for Breathing mode s = 0.The dotted blue line gives the perturbative limit of Eq. ( 8).The inset (same axis) shows that this limit is approached as The dotted lines show the hydrodynamic prediction of Eq. ( 9), see (c) for the color code.In (a) (resp.(b)) the dash-dotted line shows the sum-rule prediction for s = 0 (resp.s = 2), while the dashed grey line indicates the limit for localized excitations ωs = |µ|.There are no localized modes in the grey-shaded regions of (a) and (b).
any N/N T > 1.This breathing mode is the only localized mode present for Eq. ( 5).It may slowly decay (in a non-exponential way) because of nonlinear couplings with excitations in the continuum [50].Its frequency ω 0 can be obtained through perturbation theory for small |µ|'s [42] ω 0 = 0.95 ω * (N/N T − 1) with ω * = (g/g 12 ) 2 |g e |n ∞ .This prediction is shown in Fig. 3(a) (dotted line).In the limit of small depletion (typically up to N/N T < 1.05), it matches well the results of the Bogoliubov analysis.For larger depletions, Eq. ( 8) fails reproducing our results.This was expected since the stationary density profile predicted by Eq. ( 5) differs significantly from the exact one in this case.An upper bound for the breathing mode frequency can be obtained via sum rules [2].This general approach is known to provide acurate estimates of the excitation spectrum for superfluid binary mixtures [5].As detailed in Appendix D, a relevant sum rule is obtained by looking at the static response of the minority component to a loose harmonic potential (energy weighted sum rule).We also show the result obtained by this method in Fig. 3(a).

B. Strong-depletion regime and surface modes
For sufficiently large atom numbers, i.e. in the droplet regime corresponding to an almost full depletion of the bath, the two-component Bogoliubov analysis shows that there exist localized modes different from the breathing mode (solid lines in Fig. 3(b)).More precisely, a quadrupole mode (azimuthal number s = 2) detaches from the continuum for N/N T 3.5, see Fig. 3(b), and other modes with larger values of s emerge for even larger values of N/N T .We find that the localization for a mode of azimuthal number s ∈ N (see Fig. 3(c)) approximately occurs when the perimeter of the domain equals s times the spin healing length ξ s , which suggests an interpretation in terms of surface deformations, also called ripplons.Such ripplons are well known from 3D incompressible hydrodynamics [51].For a two-dimensional system, surface excitations of an incompressible circular bubble of radius R oscillate with an angular frequency ω s given by (see e.g.Ref. [52]) Eq. ( 9) features a linear tension coefficient T , which has a simple expression in the limit of nearby interaction parameters and equal masses m 1 = m 2 ≡ m (see e.g.Refs.[53,54]) In the short-wavelength limit, one retrieves the dispersion relation ∝ k 3/2 with wave number k = s/R expected for a linear (not-curved) interface subject to capillary waves [54].In Fig. 3(b), we show that the surface mode frequencies estimated using Eq. ( 9) asymptotically approach the frequencies obtained for large N/N T from the twocomponent Bogoliubov approach.

C. The intermediate regime: self-evaporation
For our choice of nearby interaction parameters, we found that the steady-states comprised in the range 1.45 N/N T 3.5 do not possess any localized excitation mode.Therefore, in this regime, any perturbation from equilibrium leads to the emission of mass to infinity, a dissipation mechanism known as self-evaporation.This situation is reminiscent of the spectrum of quantum droplets stabilized by beyond mean-field (BMF) effects [18,55,56], as well as of giant resonances observed in nuclear physics [57].
As discussed in Ref. [55], self-evaporation is not the dominant dissipation mechanism for BMF droplets, because of the prevalence of three-body losses in these large density systems.In contrast, for the two-component mixture considered here, the density of the localized component and thus the three-body loss rate can be tuned through the bath density.For a low-enough density, selfevaporation can then play a relevant role in the damping of the excitations of the system.It could be evaluated either solving explicitly the time-dependent nonlinear Schrödinger equations (1) or the extended RPA-Bogoliubov approach accounting for the coupling to the continuum (see e.g.Ref. [58]).

IV. CONCLUSIONS AND PERSPECTIVES
We have presented in this paper a mean-field study of the crossover from a solitonic to a droplet-like behavior in a 2D immiscible Bose mixture.We have determined both the steady-state of the system and its dynamics resulting from a small deviation from equilibrium.We have also proposed simple models that have allowed us to interpret the results obtained in the different limiting regimes.
Regarding the weak-depletion regime, we have shown in Eqs. ( 5) and ( 8) that the interaction mediated by the bath leads to a breaking of the scale invariance of the Townes soliton when its finite range is taken into account.The experimental observation of this emergent length scale should provide a way to discriminate between the scenario studied here and beyond-mean-field effects, i.e. quantum fluctuations, that also provide a mechanism for the stabilization of the minority component wave packet.
The study of the excitation spectrum of the system has revealed the existence of an interval for atom number (1.45 N/N T 3.5) over which no localized mode exists.This opens the possibility to study the intriguing phenomenon of self-evaporation, in a low-density regime for which other decay mechanisms may be minimized.
Other future directions of study include the setting in motion of the localized component [59,60], its link to superfluidity, and the emergence of a roton mode due to a capillary instability [61].Further clarification on the role of quantum fluctuations close to the miscibility threshold and its influence on the soliton formation may also provide additional interest, as recently discussed for immiscible mixtures in other configurations [62].

ACKNOWLEDGMENTS
This work is supported by ERC TORYD, European Union's Horizon 2020 Programme (QuantERA NAQUAS project) and the ANR-18-CE30-0010 grant.We acknowledge fruitful discussions with D. Petrov, G. Chauveau, F. Rabec and G. Brochier.We thank R. Kaiser for pointing out a related work in the context of photonics [63].

APPENDICES A. The two-body problem inside a BEC
We consider here two impurity atoms immersed in a 3D or 2D uniform BEC (see Fig. A1), and we derive at the lowest relevant order the expression of their effective interaction due to their coupling to the bath.The study of the coupling between one impurity and a bath formed by a BEC, the so-called Bose polaron, is a well-documented problem, see [64] and refs.in.The interaction between two Bose polarons was recently addressed in [36][37][38].We will thus keep our treatment quite brief, and focus on the specificity of the problem addressed in this article.
a. Yukawa potential for fixed impurities.We first consider two impurities of infinite mass located in R a and R b .They interact with the N 1 atoms of the bath by the contact interaction V = g 12 Using the second-quantized formalism for the bath variables, this interaction reads V = V a + V b with Here a k annihilates a particle of the bath with momentum k and Ω denotes the volume (resp.area) of the bath in the 3D (resp.2D) case.We assume that the bath is prepared in the T = 0, fully condensed state of density n ∞ = N 1 /Ω, denoted hereafter |Φ 0 with energy E 0 , and that its excitations can be described by the Bogoliubov approach.More precisely, we introduce for each momentum k = 0 the Bogoliubov operators b k , b † k which diagonalize the bath Hamiltonian such that where u k , v k are given by [1,2] (as in the main text, we set = 1).We can then rewrite the operators V j introduced above as V j ≈ g 12 n ∞ + V j with which is the Fröhlich Hamiltonian in a BEC [64].
We are interested here in the energy shift of the system that depends on the distance R ab between the particles, and that we will interpret as an effective potential energy between the two impurities.Treating V by perturbation theory, the shift of the ground state energy originating from the V j 's is up to second order in g 12 (see [38] for a systematic expansion starting from Eq. ( 11)): The first contribution is simply the mean-field interaction of each impurity with the BEC.In the second contribution, the sum runs over all excited states |Φ α of the bath.Here, only states with a single excitation k contribute and we find for the R ab dependent part of ∆E: where 1/2 stands for the Bogoliubov dispersion relation of the bath.Note that there are also contributions to ∆E that do not depend on R ab and that correspond to the self-energy of each polaron [64].
We now turn the discrete sum over k into a Ddimensional integral (the fact that the k = 0 contribution is missing in the sum of Eq. ( 16) does not play any role for our discussion).We find i.e. the Fourier transform of a Lorentzian function.It is equal to a Yukawa potential in 3D and to the potential given at Eq. ( 2) in the main text in 2D, with in both cases a range of the order of the bath healing length ξ = 1/ √ 2m 1 g 11 n ∞ .In practice, a quasi-2D gas is obtained by starting with a 3D gas and adding a strong confinement along the third direction, with a residual thickness z of the gas.The 2D version of Eq. ( 17) then holds when the distance R ab is large compared to z .
b. Effective interaction for finite-mass impurities.When the mass of the impurities m 2 is comparable to the mass of the bath particles m 1 , the kinetic energy of the impurities has to be taken into account in the calculation of the energy shift ∆E.Suppose for example that the two impurities are prepared in a state of well-defined momenta |p a , p b .For bosonic impurities, a calculation similar to the one given above leads to the second-order energy shift [38]  b − p 2 a )/2m 2 .Note that as in Eq. ( 16), we omitted the self-energy terms which do not play any role in the present work.We recover the Lorentzian momentum dependence of Eqs.(16)(17)) either if we take the limit m 2 → ∞ or, for an arbitrary value of m 2 /m 1 , in the case of a zero center-of-mass momentum p b = −p a [38].An equivalent, alternative approach consists in calculating at second order in g 12 the scattering amplitude for the collision of two impurities with momenta p a , p b in the presence of the bath [37].
c. Born approximation for the mediated interaction.For low-temperature gases, it is common to replace the "true" interaction by a (regularized) contact interaction g med δ(R), where g med is obtained by taking the zeroenergy limit of the scattering amplitude.In this limit, the contribution ∆ in Eq. ( 18), which is quadratic with respect to momenta p a,b , is negligible in front of ω k which varies linearly with k.We can then use the result (17) obtained for fixed impurities and get: which is the result used in Eq. ( 3).In the 3D case or for the quasi-2D situation where z a ij , the validity of the Born approximation used here requires the scattering length a med associated to the mediated interaction to be much smaller than its range (here ξ).Away from a scattering resonance and for g 12 ∼ g 11 , the scattering length for − g 2 12 g11 δ(r) is comparable to the van der Waals length, which is indeed much smaller than ξ for a weakly interacting gas.

B. Adiabatic elimination of the bath field
Here we consider the two coupled equations (1) for the classical fields ψ 1,2 and we explain how they can be simplified into Eq.( 5) involving only the minority component, when the density n 2 = |ψ 2 | 2 is everywhere small compared to the asymptotic bath density n ∞ .
We recall that the stationary solution of the equation for the bath field ψ 1 in the absence of the minority component (ψ 2 = 0) reads ψ 1 (r, t) = √ n ∞ e −iµ1t with µ 1 = g 11 n ∞ .Here, we treat the field ψ 2 in Eqs.(1) as a perturbation and we write the field ψ 1 as where δψ 1 is supposed to be a small correction, meaning that the bath is everywhere only weakly depleted (n 1 (r, t) ≈ n ∞ everywhere).We now detail how to reduce the initial system of equations to a single closed equation for ψ 2 : -First, by keeping all the terms in the first equation of system (1) up to order 2 in δψ 1 , we are left with: -Second, we note that the characteristic time scale for the evolution associated with the minority component is expected to be much longer than the intrinsic time scale µ −1 1 of the bath.Therefore we assume in the following that the state of the bath follows adiabatically the slow motion of the minority component, which amounts to fully neglecting the term ∂ t δψ 1 in Eq. ( 21).This approximation will be justified a posteriori at the end of this appendix.
-Third, we assume that we can perturbatively expand δψ 1 / √ n ∞ in terms of two small parameters.The first one is n 2 /n ∞ and is associated with the weak-depletion hypothesis mentioned above.The second small parameter is ξ 2 ∇ 2 and originates from the fact that the spatial variations of the fields (δψ 1 , ψ 2 ) occur on a scale much larger than the bath healing length ξ.For clarity, we now introduce the following combinations: We obtain the first-order contribution S (1) to S by keeping only the terms gathered in the first line of Eq. ( 21) (except the time derivative that we dropped): To determine D (1) , we consider the difference between Eq.( 21) and its complex conjugate: In this equation valid up to order 2, we can replace S by its first-order approximation (23), since it is multiplied by D which is itself at least of order 1.We therefore obtain at the second order in the small parameters which implies that D (1) = 0 since we consider only localized perturbations.Using the fact that δψ 1) +D (1) )/2 = S (1) /2, we can then extract the secondorder contribution to S from Eq. ( 21): -Finally, we inject the previous results into the equation giving the time evolution of ψ 2 .More precisely, we expand the density field n 1 up to second-order: This leads to Eq. ( 5), up to the contribution of the constant energy shift g 12 n ∞ .Note that the first term of the right-hand side of Eq. ( 26) that could give rise to a quadratic dependence in n 2 (quintic nonlinearity) eventually cancels with the contribution of |δψ 1 | 2 in Eq. ( 28).One may wonder if it is legitimate to fully neglect the time evolution operator ∂ t in Eq. ( 21), while keeping the laplacian operator in the perturbative expansion.This can be justified a posteriori using the dependence of the breathing frequency ω 0 with the small parameter = N/N T − 1.This frequency varies as 3/2 (see Eq. ( 8)), whereas the wave-packet size is ∝ 1/ √ .The Laplacian term ∇ 2 δψ 1 ∼ δψ 1 is thus large compared to ∂ t δψ 1 e iµ1t ∼ 3/2 δψ 1 e iµ1t and the procedure outlined here is legitimate close to the Townes threshold.However, one should expect significant corrections due to the non-adiabatic following of the bath variables as soon as N deviates significantly from N T (see also the discussion after Eq. ( 7)).

C. Excitations of the two-component system
We consider the two coupled equations (1) that give the evolution of the two classical fields ψ 1,2 , choosing for simplicity m 1 = m 2 ≡ m.We assume that the minority component contains N > N T atoms, so that there exists a stable localized state for this component.The steady-state of the system is thus characterized by the real radial wave functions R 1,2 (r).We look for perturbations around this steady-state by setting where the small perturbations α 1 , β 1 , α 2 , β 2 are by construction real functions.The evolution of the α j and β j is given by the linear system with the following differential operators The operators (L 0 , L 1 ) are deduced from (L 0 , L 1 ) by exchanging the indices 1 and 2 in these last equations.We also introduced the operator L 12 coupling the two components For r large compared to the extension of the localized component, L 12 vanishes and one can check that a localized excitation of component |2 , varying as e −κr / √ r for large r, will have a frequency ω = g 12 n ∞ − µ 2 − κ 2 /2m = |µ| − κ 2 /2m, with µ = µ 2 − g 12 n ∞ < 0, whereas a delocalized excitation varying as e ±ikr / √ r will correspond to ω = |µ| + k 2 /2m .This means that the condition ω < |µ| is a necessary condition for the excitation of component |2 to be localized.Numerically, the localized excitations as shown in Fig. 3 are identified by noticing that their frequency and functional form do not depend on the extension of the calculation grid.

D. Sum-rules
a. Monopole mode.The general formalism of sum rules provides sharp upper bounds for the excitation spectrum of many-body systems [2].We can estimate the frequency of the monopole breathing mode by calculating the ratio M 1 /M −1 between the energy-weighted and the inverse-energy-weighted sum rules relative to the operator F 0 ≡ x 2 + y 2 of the minority component.The energy-weighted sum rule is easily calculated using basic commutator rules: where the average should be taken by integrating the density of the minority component using the ground state wave function of the mixture.The inverse-energyweighted sum rule requires the calculation of the static response of the system δ x 2 +y 2 to a perturbation of the form m 2 λ 0 (x 2 + y 2 ) (again applied only to the minority component).One then obtains In conclusion, a rigorous upper bound to the frequency of the lowest monopole mode is given by b. Surface modes.The same approach can be employed to estimate the surface mode frequencies.For example, the quadrupole mode can be usefully described using the excitation operator F 2 ≡ x 2 − y 2 .In this case, the energy-weighted sum rule is still given by Eq. ( 34) holding for the monopole excitation, while the inverseenergy-weighted moment requires a calculation based on a perturbation of the type m 2 λ 2 (x 2 − y 2 ) applied only to the minority component.The result for the quadrupole frequency is then given by ω 2 2 ≤ −4λ 2 x 2 + y 2 δ x 2 − y 2 . ( It is easy to check that, when applied to a harmonicallytrapped 2D single-component BEC with repulsive interactions, the monopole and quadrupole frequencies estimated above coincide exactly with the hydrodynamic values ω 0 = 2ω ho and ω 2 = √ 2ω ho , the latter result holding in the Thomas-Fermi limit.In the calculation of the quadrupole static response, one should pay attention to the fact that the addition of the perturbation m 2 λ 2 (x 2 − y 2 ) may induce a collapse of the system at large distances, where the potential becomes deeply attractive along x (or y, depending on the sign of λ 2 ).The simplest way to evaluate the quadrupole response function while avoiding the risk of collapse is to add a perturbation of the form 2mλx 2 = mλ(x 2 + y 2 ) + mλ(x 2 − y 2 ), (38) with λ small and positive.In this way, there is no collapse at large distances and one can simultaneously calculate both the monopole δ x 2 + y 2 /λ and quadrupole δ x 2 − y 2 /λ responses, thereby giving access to both ω 0 and ω 2 with the same simulation.

FIG. 1 .
FIG. 1. Equilibrium state of a 2D minority superfluid component immersed in a bath.(a) Sketch of the interactions in the mixture, with the intra-component couplings g11, g22 and the inter-component coupling g12.Here, we consider the immiscible regime g12 > √ g11g22.(b) Two-dimensional density profile of the minority component for m1 = m2, g11 = g22 and N/NT = 1.5, where NT ≡ 5.85/|ge| is the atom number at which the Townes soliton exists and ge = g22 − g 2 12 /g11.(c) Density cuts through the center of both components for the same parameters as in (b).At a large distance from the center, the bath density takes the asymptotic value n∞.
the 2D atomic density in component |i (we set = 1).In this work, we look for configurations such that component |2 -the minority component -contains a finite number of atoms, N = |ψ 2 | 2 d 2 r, and is localized within the other component |1 -the bathwhich extends to infinity with the asymptotic density n ∞ .

6 FIG. 3 .
FIG.3.Frequencies of localized modes with azimuthal number s, for m1 = m2 and for nearby interaction parameters, here g12 = 1.01g.The solid lines in (a,b) show the predictions of the Bogoliubov approach for the two coupled equations (1).(a) Breathing mode s = 0.The dotted blue line gives the perturbative limit of Eq. (8).The inset (same axis) shows that this limit is approached as N → N + T .(b) Surface modes s ≥ 2. The dotted lines show the hydrodynamic prediction of Eq. (9), see (c) for the color code.In (a) (resp.(b)) the dash-dotted line shows the sum-rule prediction for s = 0 (resp.s = 2), while the dashed grey line indicates the limit for localized excitations ωs = |µ|.There are no localized modes in the grey-shaded regions of (a) and (b).
FIG. A1.Model studied in Appendix A. Two impurities located inside a bath formed by a Bose-Einstein condensate interact by a potential U (r) created by the emission and the absorption of virtual quasi-particles in the bath.
) where we have set k = |p a − p b | and ∆ = (p 2 1 g 11 n ∞ .In pratice, we perform an imaginary time evolution for given n ∞ and N 2 , from which we determine µ 2 .The resulting phase diagram and a few examples of steady-state density profiles are given in Figs.2(a)-2(b).Before commenting on them, we discuss hereafter various approaches that allow to draw simple physical pictures for this binary mixture.
B. The Townes soliton threshold NT