Multi-impurity polarons in a dilute Bose-Einstein condensate

We describe the ground state of a large, dilute, neutral atom Bose- Einstein condensate (BEC) doped with N strongly coupled mutually indistinguishable, bosonic neutral atoms (referred to as"impurity") in the polaron regime where the BEC density response to the impurity atoms remains significantly smaller than the average density of the surrounding BEC. We find that N impurity atoms (N is not one) can self-localize at a lower value of the impurity-boson interaction strength than a single impurity atom. When the 'bare' short-range impurity-impurity repulsion does not play a significant role, the self-localization of multiple bosonic impurity atoms into the same single particle orbital (which we call co-self-localization) is the nucleation process of the phase separation transition. When the short-range impurity-impurity repulsion successfully competes with co-self-localization, the system may form a stable liquid of self-localized single impurity polarons.


Introduction
Cold atom traps offer intriguing advantages for studying strong interaction physics in quantum many-body systems [1]. We focus on the cold atom prospect for conducting strongly coupled polaron studies. The polarons we consider are neutral atoms (referred to as "impurity atoms" in this paper) distinguishable from the BEC atoms, i.e., atoms in a different spin state or atoms of a different species embedded in a large dilute gas Bose-Einstein condensate (BEC). The BEC polaron is the impurity atom accompanied by its phonon cloud created by the interaction of the impurity with the BEC atoms. By varying an external magnetic field near a specific Feshbach resonant value [2,3,4] experimentalists can tune the strength of the impurity-boson interactions and bring the polaron into the strongly coupled regime. In that regime, the BEC-polaron self-localizes [5,6]: the ground state properties become consistent with an impurity wavefunction of finite extent. The recently demonstrated species-specific potentials (e.g. the impurity atoms), such as laser generated potentials experienced by one kind of atoms [7,8], provides unparalleled prospects for the manipulation of such polaron objects.
Attractive impurity-boson interactions give a collapsed ground state in 2D and 3D [9]. The cold atom realization of self-localized polarons with repulsive impurity-boson interactions can demonstrate polaron physics in regimes that are difficult to access in condensed matter [10]. Such cold atom experiments can then address fundamental questions pertaining to the nature of self-localization, the description of heavy polarons with a mass that is comparable to the energy scale of the boson modes and the manybody structure of systems with many polarons. At sufficiently low temperatures and sufficiently high polaron densities, cold atom traps may realize polaron superfluids without the complications that arise in describing the inevitable Coulomb interactions in electronic polarons [11].
Arguably, the polaron system provides the paradigm of a strongly interacting manybody environment modifying the properties of a single particle [10]. In the original polaron problem -the description of an electron in an ionic crystal -Landau [12] and later, Landau and Pekar [13] showed that the electron self-localizes in the strong coupling regime. In an attempt to understand the role of backflow in the roton feature of the helium superfluid dispersion, Feynman and Cohen studied a different polaron: an impurity atom in a condensed 4 He superfluid [14]. The paper by Pines, Miller and Nozieres [15] explicitly referred to this system as a polaron and pointed out the similarity between Feynman's ansatz wavefunction and the wavefunction implied by a second-order perturbation calculation. In this sense, the polarons here can be defined as the excitations created through the interaction of the impurity with the BEC atoms. Cold atom physics is now poised to realize polarons of this type at strong coupling, creating structures that are smaller than the coherence length of the BEC. These may provide a new handle with sub-coherence length resolution.
The description of impurity-BEC polarons brings in aspects of related problems. Self-localization occurs at intermediate coupling, which is the challenging regime to describe. A Feynman path integral description [16] of the BEC-impurity polaron showed a surprising, non-monotonic dependence of the polaron size with coupling strength near the coupling strength at which the strong coupling description predicts self-localization [17]. In this paper, we will work within the strong coupling description in the assumption that this description correctly captures the important trends. We describe the co-selflocalization of multiple impurity atoms in which the many-body ground state properties become consistent with N bosonic impurity atoms occupying the same single particle wavefunction of finite extent. The starting point of [5] and [6], the assumption of a product state, was designed for the description of electron bubbles in helium superfluids [18]. Below, we show that the linearization of the BEC response to the impurity wavefunction reproduces the results of the Landau-Pekar elimination of the boson modes which couple as in the Frohlich Hamiltonian to the single impurity atom in the Bogoliubov approximation of the surrounding BEC. That description shows that the energy reduction produced by the BEC density response can be cast into the form of an attractive self-interaction potential. The linearization also shows that the Frohlich impurity-BEC coupling, the Landau-Pekar elimination and the resulting self-mediated interaction fail when the BEC-density response becomes comparable to the density of the surrounding BEC. When the BEC-density is not affected much by the localized impurity, the structure of the impurity-BEC system is that of a polaron. When the BEC density is expelled by the impurity, the impurity becomes a 'bubble' separated spatially from the BEC. The physics of multiple, bosonic self-localized impurities is connected with the phenomenon of zero-temperature Bose-Einstein condensate phase separation [19], [20], [21], as we point out below. The response of a dilute BEC to a localized object was studied in the context of a moving object that is significantly smaller than the healing length of the BEC [22] and in the context of a trapped ion interacting with the BEC atoms [23] via the polarization interaction of charged and polarizable particles. The polarization potential can also give rise to self-localization of the ion as shown in [24].
In this paper, we describe the ground state of a large BEC with multiple, strongly coupled, mutually indistinguishable, bosonic impurity atoms embedded in a large dilute gas Bose-Einstein condensate. In the absence of the impurity atoms we take the BEC to be homogeneous. The boson-boson, boson-impurity and impurity-impurity interactions are described by repulsive short-range effective interactions (contact interactions). As mentioned, we assume that the impurity atoms co-self-localize: the self-localized impurity atoms do not only overlap, they occupy the same single-particle orbital. We find that multiple impurity atoms can self-localize at smaller values of the impurityboson coupling constant, than a single impurity atom can. We identify the minimal number of impurity atoms needed to stabilize the self-localized polaron structure for fixed coupling constant as a property of the droplet formation process in the nucleation of the BEC phase separation transition. When the short-range impurity-impurity interaction is sufficiently large, its competition with the self-localizing BEC-mediated interaction can lead to a regime in which one impurity atom self-localizes but not two or more.
The paper is organized as follows: in section 2, we discuss the magnitudes, dimensions and coupling constants of single BEC-impurity self-localization. In section 3, we use the Landau-Pekar method for describing the single neutral atom BEC-impurity problem so that we can later proceed beyond the homogeneous phonon description and investigate multiple impurities in the strong coupling limit. In section 4, we examine the system of strongly coupled multi-impurity boson particles, followed by section 5, where we examine the multi-impurity polaron structure variationally and discuss the connection with BEC-phase separation dynamics. Finally, in section 6 we conclude.

Single impurity self-localization in BEC
First, we consider the zero-temperature self-localization of a neutral impurity atom embedded in a dilute gas BEC. The impurity atom couples to the elementary BEC excitations and forms a polaron. Strictly speaking, the resulting object resembles a piezo-polaron since the BEC excitations are acoustic phonon modes. The term 'piezopolaron' was originally used by Mahan and Hopfield [25], [26], for conduction electrons interacting with the acoustic modes in a piezoelectric crystal [27], [28]. To bring the impurity atom to the strong coupling regime of self-localization, we assume that either the impurity-boson scattering length a IB is Feshbach tuned to a large positive value or the boson-boson scattering length a BB is tuned to a small positive value. We describe the orders of magnitude and discuss the experimental challenge of avoiding large BEC depletion when impurity self-localization occurs. Since many BECs can be described as locally homogeneous (i.e., they are in the Thomas-Fermi limit), we consider the ground state of a homogeneous BEC of N B indistinguishable boson atoms of mass m B , confined to a macroscopic volume Ω, corresponding to an average density ρ 0 = N B /Ω. Two bosons at position x i and x j effectively interact via a contact potential The BEC is dilute, so that the gas parameter γ = ρ 0 a 3 BB is much smaller than unity, i.e., γ ≪ 1. The chemical potential of the BEC bosons is µ B = λ BB ρ 0 and the BEC's healing length -the distance over which the BEC-density tends to its asymptotic value if the condensate vanishes at a planar boundary -follows from equating the corresponding kinetic energy to µ B . We will work with a coherence length ξ for which 2 The excitation of a collective phonon mode of momentum q then costs energy where c = µ B /m B denotes the BEC sound velocity. Throughout the paper, we refer to the BEC atoms as "B" and to the impurity neutral atoms as "I" . When impurity atoms of mass m I are added to the BEC, the impurity atom at r i interacts with a BEC atom at r j as well as with another impurity atom at r l via short-range interaction potentials and λ II = (4π 2 a II /m I ) and a IB , a II are the scattering lengths describing the low energy impurityboson and boson-boson scattering processes. We take all inter-particle interactions to be repulsive, i.e., a IB , a BB , a II > 0, where a BB is the scattering length of the potential λ BB = (4π 2 a BB /m B ). If the impurities are bosonic and coexist with the BEC-bosons at sufficient density and at sufficiently low temperature, the system is a mixture of a "B" and an "I" condensate. The "B"-BEC mediates interactions between the impurity atoms that are described by an attractive Yukawa potential of range ξ where we characterize the strength by an effective charge Q to emphasize the useful analogy with Coulomb interactions. The appearance of a Yukawa interaction is not surprising. It is known in quantum field theory [29] that the interactions mediated by a scalar boson field (provided by the BEC) takes the form of an attractive Yukawa potential in a non-relativistic case. In accordance, the impurity atoms in the "I"-BEC distributed at density ρ I (r) experience a local chemical potential µ I (r) with a ρ I dependent contribution if the impurity BEC is distributed homogeneously. The above chemical potential term implies a diverging "I" compressibility (∼ [∂µ I /∂ρ I ] −1 ) when λ II → 4πξ 2 Q 2 . It is the diverging compressibility that causes the zero-temperature BEC phase separation that we describe in section 5. This transition occurs when the Fetter-Colson relation [19] is satisfied λ II → λ IB (λ IB /λ BB ), where λ BB = (4π 2 a BB /m B ) so that is the value of the Yukawa strength parameter.
The BEC mediated interaction is a consequence of BEC density correlations: an impurity atom i at r i causes a density variation in the BEC that is experienced by impurity atom j at r j as a change in its BEC mean field energy. Describing the impurity atoms by a quantum wavefunction, we find that even a single impurity is affected by BEC density correlations. The impurity interaction energy density at r is influenced by the impurity density at r ′ as the BEC responds to the density of the entire impurity wavefunction. As we will see below, a single impurity then experiences an effective self-interaction V s (r − r ′ ) = − Q 2 e −|r−r ′ |/ξ / |r − r ′ |. The integration of V s (r − r ′ ) over r and r ′ weighted by 1 2 ρ I (r)ρ I (r ′ ) gives the gain in impurity-boson interaction energy by the adjustment of the BEC density to the impurity wavefunction.
In the strong coupling limit, the effective self-interaction exceeds the kinetic energy cost of localizing the impurity. For the Yukawa self-interaction to overcome the kinetic energy, the impurity has to localize to a size comparable to or less than ξ. In that case, the self-interaction is Coulomb-like in most of the impurity-region and, hence, is efficient at binding. As in the Hydrogen atom description, the length scale R 0 is determined by comparing Coulomb and kinetic energies Using equation (4), µ B = 2 /[4m B ξ 2 ], and ξ −2 = 16πρ 0 a BB , we find that the effective Rydberg length is equal to Aside from the mass factor 1 + m I , R 0 is the mean free path of an impurity atom that encounters hard-sphere scatterers of radius a IB distributed at average density ρ 0 . The energy scale of impurity self-localization is given by the corresponding Rydberg energy The mediated self-interaction can only bind with sufficient efficiency if the natural impurity size R 0 is comparable to or shorter than ξ. The ratio of the BEC-coherence and Rydberg lengths, ξ/R 0 , is the only parameter left in the description of the impurity wavefunction after scaling the energy and length by E 0 and R 0 , respectively. We call this dimensionless ratio, the impurity-boson coupling parameter, Calculations [5] indicate that a single impurity atom can self-localize if β > β c (1) = 4.7, where β c (1) stands for the critical β for a single impurity atom. As β further increases, the self-localized impurity size shrinks to become a point-like object (smaller than the BEC coherence length). The calculation of a point-like potential moving with velocity v [22] gave an energy a BB , so that we might expect the effective polaron mass m * to be equal to m * = m I + m D or m * m I = 1 + 2β However, the calculation in [22] was based on perturbation theory so that its validity does not extend to the large coupling limit. For the classic polaron treatment that describes an electron coupled to optical phonons, Feynman estimated an effective mass that varies as the fourth power of the coupling constant [10]. It is useful to express the observables in terms of the impurity-boson coupling constant β. For instance, with τ 0 = /E 0 representing the time scale relevant to impurity self-localization. The polaron description of the BEC-impurity assumes that the phonon modes are excitations of a homogeneous BEC. However, as the impurity-boson interaction is increased to reach the critical coupling β c (1), the impurity can deplete the "B" condensate in its vicinity. The above description breaks down when the local BEC-density variation δρ B becomes comparable to ρ 0 . For repulsive impurity-boson interactions we expect the single impurity wavefunction to take on a phase-separated bubble profile when δρ B ≃ ρ 0 . While this is interesting by itself, the impurity loses its Frohlich-coupled polaron character in the bubble limit. As the self-interaction energy is the change in boson-impurity energy caused by the impurity induced variation in where γ represents the BEC gas parameter, γ = ρ 0 a 3 BB , defined earlier in this section (above equation (1)).
Can the self-localization of a single BEC-impurity be realized in cold atom experiments? At the end of the next section, we discuss whether cold atom technology can access the polaron regime -the parameter region where the density variation, δρ B , in equation (12) remains significantly smaller than ρ 0 . Here we discuss the experimental accessibility of the strong impurity-boson coupling regime, β > 5. Three experimental challenges must be met. Two of them relate to the requisite increase in impurityboson scattering length. First of all, if that increase is effected by a magnetically controlled Feshbach resonance, the magnetic field has to be sufficiently homogeneous and constant to avoid significant variations of the impurity-boson interactions. Secondly, three body recombination into the large, two atom impurity-boson bound (dimer) state that becomes degenerate with the incident channel impurity-boson continuum as a IB diverges, should not occur before the BEC impurity has self-localized. As a third condition, we mention that the BEC-temperature should be sufficiently low to prevent thermal fluctuations from de-self-localizing the strongly coupled BEC impurity.
How large does a IB have to be? Estimating the necessary a IB from with γ ∼ 10 −4 and a realistic mass ratio of ∼ 10, we find that a IB a BB ∼ 10 2 if β ∼ 10. The Feshbach scattering length with magnetic field dependence is given by where a IB,0 denotes the background scattering length, ∆ is the resonance width, and B 0 is the resonant field strength. Hence, a slow magnetic field variation δB induces an a IB variation Requiring |δa IB /a IB | < ε leads to |δB/∆| < ε (a IB,0 /a IB ). Assuming, for instance, a IB,0 /a IB = 1/30, ε ∼ 0.03, we obtain |δB/∆| < 10 −3 . For ∆ = 1G, the magnetic field variation is δB < 1mG, which can be achieved in today's laboratories. In contrast, a broader resonance may require δB < 10mG (∆ = 10G) or δB < 0.1G (∆ = 100G). We now estimate the recombination time of the BEC impurity atom in the worst case scenario. Applying the same reasoning as in [30] where the three-body recombination rate for a BEC of large scattering length was first calculated, we expect a recombination rate of the form For a BEC of large and positive a BB , substituting m I → m B and a IB → a BB in equation (17) gives an overestimate of the recombination rate (see [31]). The authors of reference [31] found that the a 4 power has to be multiplied by an oscillating function of magnitude less than unity. Near a node of the oscillating function, the recombination rate is much reduced. Nevertheless, using equation (17), we find where τ B = /µ B represents the time scale on which the BEC can respond to the impurity and form the phonon cloud that accompanies the impurity in the polaron state. Note that for m I ∼ 10m B , β ∼ 10, τ R can exceed τ B by more than two orders of magnitude. Finally, we comment on the temperature requirement for observing BEC impurity self-localization. Determining the polaron properties in a finite temperature Feynman path integral calculation, reference [16] identified polaron self-localization with a sudden and unexpected increase in polaron size as β increases near the critical value β c obtained from the strong coupling description. The polaron size increase was observed for a 6 Li-impurity in a 23 Na-BEC in the case where the temperature (k B T ) was smaller than half the boson chemical potential. Therefore, we expect that self-localization takes place when k B T < µ B f , where f can be as large as 0.5. For a BEC of µ B = × 300Hz, the temperature requirement of T < 7.5nK can be easily achieved in today's cold atom experiments.

Landau-Pekar elimination of phonon modes
where ψ and χ denote the normalized single-particle impurity and BEC boson wavefunctions, respectively. In the absence of an impurity, the BEC-field ϕ (r) = √ N B χ (r) is ϕ = √ ρ 0 , position independent for a homogeneous BEC. Assuming that the BEC response to the impurity remains smaller than √ ρ 0 , ϕ = √ ρ 0 + δϕ, the equations can be linearized in δϕ. The resulting description is equivalent to the Landau-Pekar treatment we describe below. The product wavefunction reveals the connection with the analogous "fluid with embedded droplet" system. The linearization procedure indicates how to proceed beyond the homogeneous phonon description. On the other hand, the Landau-Pekar method is most easily generalized to describe multiple impurities in the strong coupling limit.
In the absence of an impurity, the excitations of the dilute BEC are well-described in the Bogoliubov approximation. This description transforms the operators that create and annihilate boson particles of momentum q,b † q ,b q to quasi-particle (phonon) creation and annihilation operatorsâ † q ,â q , that diagonalize the linearized Hamiltonian. Adding a constant term, the transformed BEC-Hamiltonian reads where ω q is stated in equation (1). The first step of the Bogoliubov procedure replaces theb † q=0 ,b q=0 -operators by √ N B and expands the energy in powers of √ N B . Accordingly, the boson density operator becomeŝ where the last step followed from implementing the Bogoliubov transformation. Introducing the impurity density operatorρ I,q , the impurity-boson interaction takes the form wheren I represents the number operator of the I-atoms,ρ I,k=0 =n I . The second term on the right hand side resembles the electron-phonon interaction term in the Frohlich where The phonon part of the impurity-BEC Hamiltonian,Ĥ p , readŝ The strong coupling approximation assumes that the impurity is described by a welldefined wavefunction that is not entangled with the boson B state. As a consequence, the coupling to other impurity states can be neglected and the impurity density operator ρ I,q can be replaced by its expectation value,ρ I,q → ρ I,q = d 3 r exp(−iq·r) |ψ (r)| 2 . To bring out the oscillator nature of the boson modes, we introduce the oscillator coordinate and momentum operators that are equivalent to the creation and annihilation operators, In this notation, the phonon HamiltonianĤ p becomes a sum over second order polynomials in φ, Completing the squares in equation (24), we write the phonon Hamiltonian as a sum over displaced oscillator Hamiltonianŝ where theΦ-operator denote the displaced boson oscillator coordinateŝ As the displaced ground state energy is the same as that of the original oscillator, the energy reduction caused by the phonon response to the impurity density is given by the remainder Using ρ I,q = d 3 r exp [iq · r] |ψ (r)| 2 this expression can be cast in the form of an interaction energy where V s represents the inverse Fourier transform of −2M 2 0 ν 2 q / ω q . The expression holds regardless of the momentum dependence of ω q and ν q which are different for ionic crystal (optical) and piezoelectric (acoustic) polarons. In the case of the BEC-impurity or "Bogoliubov polaron", V s is the attractive Yukawa interaction, where we used equation (4) to obtain the above expression. The total change in ground state energy E caused by doping the BEC with a strongly coupled impurity atom, includes the kinetic energy cost of localizing the impurity is The optimal wavefunction ψ follows from minimizing E while ensuring that ψ (r) is normalized. We include the normalization constraint by introducing a Lagrange multiplier µ. The extremum of the associated functional F = E − µ d 3 r |ψ (r)| 2 − 1 is found by requiring δF/δψ * (r) = 0, which leads to where µ I = µ − λ IB ρ 0 . The impurity wavefunction that solves equation (31) is localized if and only if µ I is negative. Multiplying equation (31) by ψ * (r), integrating over r, and making use of the normalization condition, we obtain which differs from equation (30) by the factor of 2 in the self-interaction term. The numerical solution of equation (31), as calculated iteratively in [5], shows that the profile of ψ generally resembles a Gaussian function. Assuming a normalized Gaussian wavefunction, ψ (r) = e −|r| 2 /2σ 2 / (πσ 2 ) 3/4 , we solve the problem variationally by minimizing E with respect to σ. To know if the variational description predicts selflocalization, we test if at the minimal widthσ, µ I (σ) < 0. Working out the integrals we find where f , accounts for the reduced efficiency of the Yukawa self-interaction when the Gaussian width exceeds the BEC-coherence length (as f (y) → 0 when y ≫ 1). If we scale the energies by the Rydberg energy of equation (7), E → (E/E 0 ), µ I → (µ I /E 0 ), and the length by the Rydberg length of equation (6), x = σ/R 0 , β = ξ/R 0 , then E and µ I take the form   Figure 2. The variationally determined chemical impurity potential µ I . In the strong coupling description, the BEC-impurity polaron self-localizes when µ I < 0.
The BEC density response to the self-localized impurity is a consequence of the oscillator (φ) displacement in the Landau-Pekar description. We calculate the density variation in the Bogoliubov approximation, Using Φ k = 0 we find Recognizing the expression in the square brackets as the Fourier transform of V s and replacing √ ρ 0 /M 0 = 1/λ IB , we obtain which is identical to the expression found in [5] by linearizing the BEC-response. Note that the Landau-Pekar description and the linearization of the product state give the same results for all observables. We now determine the peak value of the BEC density variation. If we choose the origin at the center-of-mass of the impurity in |ψ (r)| 2 , the peak density is reached at x = 0. In the true ground state, the impurity is not localized in any specific localization but the strong-coupling product wavefunction breaks translational symmetry. While the wavefunction does not have the correct symmetry, we know that some of the relevant polaron properties such as energy, mass, and polaron size give the correct values in the strong-coupling limit as they do in the traditional optical polaron problem. Assuming the width to take on the strong coupling value, σ = 3 π/2R 0 , where the factor E 0 µ B λ BB λ IB is given by equation (12). Assuming 3 √ π/2β ≪ 1, we find In the polaron limit, the limit in which the Bogoliubov-approximated impurity-boson coupling of the Frohlich type is valid, the BEC-density response should be smaller than the average BEC density, δρ B /ρ 0 < ǫ, where ǫ is a fraction of unity. We find that the average boson BEC-density should not exceed a value ρ max with Depending on the atoms trapped in the experiment, the above condition can pose a challenge. Taking the impurity atom to be a spin-flipped atom in a 23 Na-BEC, choosing ǫ = 1/5 and assuming that a IB is Feshbach tuned while a BB = a T riplet = 3nm, we find ρ max = 3.6 × 10 10 cm −3 , which is a rather low density to make observations with. On the other hand, embedding a 87 Rb-impurity atom into 7 Li BEC of a BB = 0.5nm, with ǫ = 1/5 leads to ρ max = 1.2 × 10 15 cm 3 [β/β c (1)] −6 , yielding a maximal density of 2 × 10 13 cm −3 even for β = 10.

Co-self-localization of bosonic BEC-impurities
Most schemes for realizing strong coupling BEC-impurity polarons would create multiple impurities that are indistinguishable in the same cold atom trap. In this section, we describe the ground state of the N impurity system in a homogeneous BEC. We assume that the N impurity atoms are indistinguishable bosons and that the impurity-impurity interactions λ II δ (r i − r j ) are repulsive, λ II > 0. We will characterize the strength of the short-range impurity-impurity interactions by a third dimensionless coupling constant The system's behavior is determined by three dimensionless parameters: the BEC gas parameter γ = ρ 0 a 3 BB , which affects the BEC-density response to the self-localized impurity but not the impurity wavefunction in the polaron limit, the impurity-boson coupling constant, β = √ π ρ 0 a 4 IB /a BB (1 + m B /m I ) (1 + m I /m B ), and the above defined impurity-impurity coupling constant, which can also be written as We extend the above descriptions to treat N impurities of position r i , i = 1, 2, . . . N, and N B bosons of position x i , i = 1, 2, . . . N B . The strong coupling ground state wavefunction can be written as a product state Ψ (r 1 , r 2 , . . . r N ) χ B (x 1 ) χ B (x 2 ) . . . χ B (r N B ). Alternatively, the Landau-Pekar treatment is easily generalized as the impurity-boson coupling is described by the same Frohlich-like term of equation (25). In the strongcoupling limit, the density operator is replaced by its expectation valueρ q → ρ q where The strong coupling Landau-Pekar phonon elimination leads to a phonon-induced ground state energy reduction We can argue that the i = j terms describe self-interaction while the i = j terms describe the BEC-mediated interactions. However, this distinction becomes less meaningful when the impurities are indistinguishable and Ψ satisfies the bosonic or fermionic permutation symmetry. In this case, all i, j contributions are equal and ∆E p becomes The localization of the impurities comes at a kinetic energy cost and the impurityimpurity overlap comes at an interaction energy cost. We describe the direct impurityimpurity interactions by an effective potential v II (r i − r j ) = λ II δ (r i − r j ). Adding N strongly coupled impurities to a homogeneous dilute BEC changes the ground state energy by an amount E where The mediated inter-particle interaction potential, V s , is attractive (V s < 0) and favors impurity overlap. Maximal overlap is achieved when the N bosons occupy the same single particle orbital, ψ (r 1 , r 2 , . . . r N ) = ψ (r 1 ) ψ (r 2 ) . . . ψ (r N ). To emphasize that the boson impurity particles not only self-localize, but occupy the same orbital (localizing in Hilbert space), we will refer to the co-self-localization of the bosonic impurity atoms.
In that case, ρ I (r) = |ψ (r)| 2 and using equation (47) in the E expression we find where ψ (r) is assumed to be normalized. In determining the optimal ψ-orbital, we introduce a Lagrangian multiplier Nµ to ensure d 3 r |ψ (r)| 2 = 1. Minimizing E − Nµ d 3 r |ψ (r)| 2 − 1 with respect to ψ * yields where µ I = µ − λ IB ρ 0 . The initial N-factor in the Lagrange multiplier ensures that µ I becomes the impurity chemical potential in the limit of large impurity numbers N.
Solving the above equation for µ I by multiplying both sides by ψ * (r) and integrating over r, we obtain The physical µ I value minimizes E in equation (49). For negative values of µ I , ψ becomes a localized wavefunction. At a distance exceeding ξ from where the impurity density becomes small, ψ takes on the Yukawa profile of a bound state.

Variational description of BEC-impurity co-self-localization
Alternative to solving the non-linear Schrodinger equation of equation (50), we can approximate the impurity orbital ψ(r) variationally. As in section 3, we insert a normalized Gaussian wavefunction, ψ(r) = e −r 2 /2σ 2 (πσ 2 ) −3/4 , into the energy E of equation (49) and we minimize E with respect to the Gaussian width σ. The treatment predicts self-localization of the N co-locating impurities if µ I of equation (51), evaluated at σ that minimizes E, gives a negative value.

The variational functionals
In the scaling units, E 0 , R 0 , introduced above to treat the single impurity polaron problem, the E and µ I expressions take the form and where x = σ/R 0 and where the impurity-impurity coupling was defined above (see equation (44) Two interesting and relevant limits can be understood analytically: the 'pinprick' limit in which the size of the impurity wavefunction shrinks well below the coherence length of the BEC and the large impurity number limit in which the impurity wavefunction extent is determined by the competition of the BEC-mediated interaction and the impurity-impurity repulsion.

The pinprick polaron limit
The coherence length ξ is the natural length scale of the BEC. If an object has a subcoherence length size L, i.e., L ≪ ξ, the BEC response to its presence is indistinguishable from the BEC-response to a point object. In the presence of N coself-localizing boson impurities, the BEC-mediated impurity-impurity attractions that contribute a term ∼ N 2 can further decrease the impurity size to be sub-coherence length -the impurity becomes a 'pinprick'. In this limit, the f argument, x/β = σ/ξ, is much less than 1 and f can be Taylor-expanded in the E and µ I expressions, and where E c and µ c denote the Gaussian expectation values for Coulomb (instead of Yukawa) self-interactions, In this limit, the equilibrium Gaussian width follows from the requirement ∂E c /∂x = 0 at fixed N and η, yielding x (N, η). The critical coupling β c (N, η) follows from µ I = 0; although this approach becomes somewhat questionable for determining the impurity properties near the critical coupling since the Taylor expansion is not generally valid near the self-localization point. Furthermore, if the short-range impurity-impurity interactions can be neglected (η → 0) in the pinprick multi-impurity polaron, the condition, ∂E c /∂x = 0 leads to the equilibrium width x = 3 π/2 /N ≈ x (N = 1, η = 0) /N. The expectation value of E c scales as N 3 and the binding energy takes the form On the other hand, the expectation value of µ c scales as N 2 , µ c /E 0 = −N 2 /π and from equation (57) we estimate the critical coupling constant as Hence, N bosonic impurities can co-self-localize at much smaller values of the bosonimpurity coupling strength than the critical coupling for single-impurity self-localization. The self-localized multi-impurity polaron can be significantly smaller than the BECcoherence length. Specifically, casting the impurity wavefunction extent, r 2 = 3/2σ = 3/2xR 0 and the binding energy in terms of the BEC-length and energy scales, we find For small, but finite impurity-impurity coupling η, the impurity-impurity interaction becomes significant as N increases. Solving for the equilibrium width x from ∂E c /∂x = 0, we obtain x = 3 π/2 (1/2) 1 + 1 + (2/3) ηN(N − 1) 2/π , so that the impurity-impurity repulsion significantly increases the extent of the impurity wavefunction when N(N − 1) ≥ 3 π/2/2η . In figure 3 we show the Gaussian width of the three-impurity polaron wavefunction at coupling β = 4 as a function of η, as determined numerically from the full functionals, equations (52) and (53), and as determined from the Taylor expanded approximation equations (54) and (55). Notice that the Taylor expansion only works well for small values of η as the impurityimpurity repulsion sufficiently increases the impurity radius to invalidate the Taylor expansion assumption x/β ≪ 1. As the expansion replaces the Yukawa self-interaction by a Coulomb attraction (more efficient at binding the polaron), the E c approximation underestimates the size of the impurity wavefunction. The effect of the impurity-impurity repulsion is particularly striking in the dependence of the critical β c on the impurity numbers N, as shown in figure 4. This figure shows β c (N, η) for N = 1 (independent of η), 2, 3, ..., 10, at η = 1, 10, 30 and 40. At η = 0, β c (N, η = 0) decreases with increasing impurity numbers and β c (N, η = 0) vanishes at large N. At finite, fixed η with η < 24 the critical coupling β c (N, η) (shown by dots in the same color) still decrease with increasing N, but β c (N, η) tends to a finite value β C (N → ∞, η) that depends on η. At η > 24, β c (2, η) > β c (1), but then β c (N, η), N ≥ 2 still falls off monotonically with increasing impurity numbers N.
This dependence has interesting implications. If the BEC-impurity system has an impurity-boson coupling constant β in the range [β c (N β − 1), β c (N β )], the above description suggests that N β or more bosonic impurities would co-self-localize in the ground state, whereas N β − 1 or less impurities would not. The minimal number of particles required to make a stable, self-localized structure is reminiscent of the familiar nucleation processes that take place in first-order transitions. In that case, the formation of a stable droplet of the new phase also requires a minimal amount of matter in the new phase to overcome the relative cost of the surface energy (here, kinetic energy plays the role of surface energy). We tentatively identify the co-self-localization of impurities as the onset of the nucleation process of B-I phase separation, a connection we will further explore below.  Figure 4. Critical values of the impurity-boson coupling constant β, β c , above which N bosonic impurity atoms self-localize in the BEC-impurity ground state described in the strong coupling approximation. The graph shows β c for different impurity atom numbers N , N = 1, 2, . . . , 10 , set out on the horizontal axis, and gives the variationally determined β c -values for different impurity-impurity coupling values, η: red diamondsη = 0, black plusη = 1, black squareη = 10, black xη = 30 and blue circleη = 40.

BEC-phase separation
In a phase separation transition, initially miscible fluids become immiscible and separate from each other. At zero temperature, the ground state of a gas of I and B weakly interacting boson particles changes from a mixed to a separated configuration as the impurity-boson scattering length is increased. The natural coupling constant to characterize this zero temperature phase transition is given by As g varies from g < 1 (miscible mixture of I-and B-BECs) to g > 1 (ground state corresponding to separated BECs), the system undergoes a phase transition that appears to be second order -the compressibility and a coherence length diverge. If the I BEC is the smaller condensate, the ground state gathers the I-bosons into a spherical region as this geometry minimizes the surface energy. The surface tension itself vanishes as g → 1 + , but is finite when g > 1. Hence, we expect separation dynamics at g > 1 to be similar to that of a first-order transition.
For large I-and B-BECs, the onset of the separation dynamics is described by the growing instability of the collective excitations of the homogeneous BEC-mixture. But how does the separation dynamics set in if only a few I-atoms are present and the I-bosons have not acquired macroscopic phase coherence, whereas the temperature is sufficiently low to treat the B-BEC as a zero temperature system? We suggest that when fluctuations gather a minimal number N β of I-bosons in a B-coherence length size region, the impurities can co-self-localize into a self-localized multi-impurity polaron (at least in the parameter region in which the BEC density response to the impurities remains small). Further collisions with more impurity atoms grow the self-localized polaron impurity numbers. Finally, the relative B BEC density response becomes large, the B-BEC density is expelled from the impurity region and the impurity extent grows to exceed the B-coherence length. Below, we back up this picture by showing that the condition for first order impurity nucleation coincides (or nearly coincides) with the condition for zero temperature BEC phase separation.

The limit of large co-self-localizing impurity numbers
As the impurity number N increases at fixed β coupling, the BEC density response to the co-self-localizing impurity polaron grows. For an impurity number N D , the BEC density becomes comparable to the average BEC-density and the self-localized impurity structure crosses over to a 'bubble' instead of a sub-coherence length polaron. In that regime, the Frohlich-coupling of equation (22), which was derived in the Bogoliubov approximation and assumes low depletion, does not accurately describe the impurityboson coupling and the Landau-Pekar treatment breaks down. Hence, the large impurity limit, N → ∞, of β c further strengthens the connection between multi-impurity polaron formation and BEC phase separation. Keeping only the N 2 -terms of the E of equation (52), the kinetic energy term drops out and the polaron size is determined from the competition of the BEC B mediated interactions and the short-range impurity-impurity interactions, where Introducing y = x/β as variable, y = σ/ξ, the energy functional takes the form where A = N 2 / [2β 2 η] > 0 and κ = 8/πβ 2 η −1 . This function has a minimum in y and the local minimum has a negative value provided κ exceeds a minimal value κ c . Numerically, we found κ c = 1.0004. When κ > κ c , the value of µ I at the minimum E is negative as µ I = 2E/N in the absence of the kinetic energy term. Therefore, a large number N of bosonic impurities co-self-localize if κ > κ c or β > (π/8) 1/4 √ η √ κ c = β c (N → ∞). We assume that the nucleation occurs when a number, any number, of impurities can co-self-localize. As β c (N) decreases with N, β c (N → ∞) represents the minimal coupling value at which the transition takes place. Therefore, a transition sets in if where we used η = 8 √ 2βγ(a II /a BB ). As a consequence, the nucleation condition takes the form replacing β = √ π ρ 0 a 4 IB /a BB (1 + m I /m B ) (1 + m B /m I ) and γ = ρ 0 a 3 BB , we obtain a 2 IB a II a BB which is very close to the phase separation condition of equation (61) as κ c = 1.0004. What is N D at which the self-localized polaron structure crosses over to a phase separated bubble? The peak BEC density variation for a Gaussian impurity wavefunction takes the form In the pinprick regime, x/β ≪ 1 with weak impurity-impurity interactions, N(N −1) ≪ 3 π/2/ [2η] so that x = 3 π/2/N, we obtain Estimating N D from |δρ B (x = 0) /ρ 0 | = 1, we find The N 2 -scaling of δρ B /ρ 0 may convey the message that multi-impurity polarons with a finite number of impurities is unlikely. However, at larger impurity number, the selflocalization also sets in at a smaller value of N. The relative BEC density variation near self-localization with β ≈ β c ≈ (2π) /N is which scales as √ N and can be satisfied experimentally for N = 1.

Self-localized impurity polaron fluid
It is interesting to note that if a II is replaced by ξ in the phase separation coupling constant g, the expression goes over into the expression for the impurity-boson coupling constant β, At phase separation, g = 1 and β = β c (N → ∞) so that, up to a factor κ c /1 a II ξ = β c (N → ∞) .
From figure 4, it is clear that when η is sufficiently large to make β c (N → ∞) comparable to or greater than 5, (implying a II ≫ ξ), then β c (1) < β c (N > 1). What is the ground state of an impurity-BEC system with multiple impurities and an impurity-boson coupling constant β with β c (1) < β < β c (N → ∞)? The conditions suggest that one impurity atom self-localizes but that the co-self-localization of two or more impurities is prevented by the short-range impurity-impurity repulsion. Each impurity could selflocalize, but how do the self-localized single impurity polarons arrange themselves? As η is large, two self-localized polarons repel each other at small distances but attract each other at larger distances ranging up to ξ. In classical particle systems, these are the conditions at which a crystal can form or a liquid of self-determined density. Would the self-localized polaron liquid form a self-localized polaron crystal, liquid or BEC? This regime is, however, not easily realized experimentally. In the absence of any resonances, the cold atom BEC coherence length generally exceeds the scattering length by one or two orders of magnitude. The realization of a self-localized single impurity polaron liquid that does not nucleate into a phase separated BEC, would require the enhancement of both the a IB and the a II . If one Feshbach resonance enhances a IB to increase β, the magnetic field is 'used up' and a II cannot be varied magnetically. We suggest that in two-dimensional systems confinement-induced resonances can be used to enhance the other scattering length.
The confinement of atoms to one or two dimensions by squeezing the trap potential in the other dimensions induces resonances in the atom-atom interactions. These occur when the extent of the wavefunction in the transverse (squeezed) direction becomes comparable to the scattering length in three dimensions [32], [33]. As we expect one and two-dimensional BEC or quasi-BEC systems to undergo phase separation as well as exhibit self-localization of impurity atoms, we suggest that a combination of a magnetically controlled Feshbach resonance and a confinement-induced resonance may access self-localized single impurity polaron fluids.

Conclusions
The experimental observation of multi-impurity self-localized polarons would provide an interesting class of structures that have a subcoherence length size. BEC experiments can address fundamental, strong coupling questions: Does polaron self-localization give transition like behavior? What is the effective mass of self-localized BEC impurity polarons? Can the impurity extent, as well as the polaron size be measured directly by species specific dipole potentials? In addition to fundamental many-body research the realization of these structures may bring new methods to bear for manipulating and probing BEC systems. Self-localized multi-impurity polarons can be transported by a laser beam of moving but broad focus. Can the movable self-localized polaron carry qubit information? The creation of shock wave analogues in BEC's required density perturbations that were smaller than the BEC-coherence length [34], can other phenomena be triggered such as superradiance in phonon radiation? The realization of a Josephson junction in BEC's described by the usual linear Josephson coupling Hamiltonian requires a constriction that is smaller than the BEC-coherence length. Can a trapped self-localized impurity provide such constriction?
If many impurity atoms are present in the trap, the self-localized impurity polarons may represent transient structures. The self-localization dynamics, we suggest, describes the onset of the phase separation nucleation process. On the other hand, existing cold atom techniques can isolate the self-localized impurity droplet with a controlled number of impurities. For instance, a species specific optical lattice can contain a well determined number N of impurity atoms per site. An overlapping BEC does not feel the optical lattice but can self-localize the N-impurities. As self-localization depends on the number of co-located impurity atoms, this technique could provide another sensitive method for counting atoms in each well.
In summary, we have shown that N bosonic neutral impurity atoms embedded in a dilute gas BEC can self-localize into a collective strongly coupled polaron state. The critical impurity-boson interaction strength for self-localization can be significantly lower than the corresponding value for a single impurity polaron. We identify the self-localized impurity polarons as the initial droplets of the BEC phase separation nucleation process. In the limit of large impurity-impurity interactions, there may be a strongly interacting regime in which self-localized single impurity polarons form a stable self-localized polaron liquid. We have discussed the experimental challenge of realizing these systems in the polaron regime and not as phase separated bubbles. Finally, we have speculated on possible applications of this new class of sub-coherence length BEC structure.

Acknowledgments
ET's work was funded by the LDRD-program of Los Alamos National Laboratory. Some of this work was performed at the Aspen Center of Physics.