Ionic polaron in a Bose-Einstein condensate

The presence of strong interactions in a many-body quantum system can lead to a variety of exotic effects. Here we show that even in a comparatively simple setup consisting of a charged impurity in a weakly interacting bosonic medium the competition of length scales gives rise to a highly correlated mesoscopic state. Using quantum Monte Carlo simulations, we unravel its vastly different polaronic properties compared to neutral quantum impurities. Moreover, we identify a transition between the regime amenable to conventional perturbative treatment in the limit of weak atom-ion interactions and a many-body bound state with vanishing quasi-particle residue composed of hundreds of atoms. In order to analyze the structure of the corresponding states, we examine the atom-ion and atom-atom correlation functions which both show nontrivial properties. Our findings are directly relevant to experiments using hybrid atom-ion setups that have recently attained the ultracold regime. An impurity introduced to a many-body quantum environment gets dressed by excitations and it is of a particular interest to understand the limits of the quasi-particle description. The authors theoretically and numerically study an ionic impurity immersed in a weakly interacting gas of bosonic atoms and demonstrate the existence of two main phases of a polaronic regime for weak interactions, and a strongly correlated state with many bosons bound to the ion.

A n impurity immersed in a many-body quantum system constitutes a fundamental building block in condensedmatter physics, particularly with regards to transport properties of materials 1,2 . In order to investigate this paradigmatic problem, ultracold atoms are especially suited as they allow for experimental control of multiple parameters such as the shape of the confinement and the form of interparticle interactions 3 . Over the last decade, atomic systems have proven to be capable to observe the formation of quasi-particles such as the celebrated polaron in bosonic 4-6 and fermionic [7][8][9][10][11][12] quantum gases as well as in a Rydberg system 13 , with the possibility of exploring impurity physics also in the presence of dipolar interactions 14 .
Within the plethora of compound atomic quantum systems available, atom-ion systems provide a unique arena for investigating many-body quantum physics in the strongly interacting regime. Indeed, the interaction between the charge and the induced dipole of the neutral particle results in the asymptotic form Importantly, this polarization potential has a characteristic length scale that is about an order of magnitude larger than in the case of van der Waals interactions typical for neutral atoms and it can become comparable to the mean interparticle distance. Moreover, the characteristic interaction energy is typically in the microkelvin range and thus comparable to experimentally achievable collision energies 15 . Hence, although the polarization potential has in principle short-ranged nature one cannot replace it with a pseudopotential due to the lack of separation of length and energy scales between the two-body and the manybody systems. While this lack of scale separation gives rise to a striking competition of few-body and many-body processes, it poses a theoretical challenge due to the necessity to account for details of the potential that severely inhibits the possibility of using analytical methods to describe the properties of an ionic impurity such as its effective mass.
Recently, it has been demonstrated that for certain atom-ion combinations the ultracold regime is within experimental reach [16][17][18] . Exploiting the charge of the impurity, one appealing possibility is to study transport properties by dragging the ion by means of electric fields and detect it with high spatial and temporal resolution 19 . The long range of the atom-ion potential, on the other hand, allows to the investigation of the formation of mesoscopic bound states [20][21][22] that would have a dramatic impact on the ion transport dynamics.
Studies of the mobility of an ion moving in a bosonic medium date back to the early 1960s with the aim of explaining the small ion mobility in liquid helium [23][24][25][26][27] . Later, mean-field approaches were used to predict the formation of mesoscopic molecular ions 20 , to estimate the number of excess atoms around an ion in a homogeneous Bose-Einstein condensate (BEC) 28 , and path integral methods to determine the ion polaron properties in the strong-coupling regime 29 . While such approaches allow one to obtain a qualitative picture of the underlying phenomenology, the interplay between the long-ranged potential and strong interactions not only leads to substantial shifts to the relevant observables such as the energy of the system but also has drastic consequences for the structural properties of the ground state. In line with the recent experiments with ion-atom mixtures 16,[30][31][32][33][34][35][36][37][38] , in this letter, we study the many-body ground state properties of an ion immersed in a three-dimensional (3D) bosonic gas. To this end, we employ quantum Monte Carlo techniques that have been successful in the context of the 3D Bose polaron [39][40][41] , bipolarons 42 , as well as in two-dimensional (2D) 43 , and one-dimensional (1D) polarons [44][45][46] . The method allows us to fully take into account the quantum many-body correlations that turn out to be important for predicting how the competition between the few-body and many-body length scales gives rise to a striking impurity physics that is governed by a transition from a polaron to a many-body polaron bound state. The resulting states cannot be captured by conventional tools such as the Frölich model or Bogolyubov theory. The system properties depend not only on the scattering length and effective range of the two-body potential but rely on the presence of the long-range tail of the interaction, indicating the failure of the pseudopotential approximation.

Results
System. We consider an ion of mass M immersed into a gas consisting of N bosonic atoms of mass m at average density n = N/L 3 . For simplicity, we focus here on the mass-balanced case (i.e. M = m). We consider periodic boundary conditions in a box of size L chosen large compared to the healing length ξ = (8πna aa ) 1/2 , where a aa is the boson-boson s-wave scattering length.
The microscopic many-body Hamiltonian is given bŷ Hereafter, we denote the ion's characteristics such as position and mass with capital Latin letters, while for atom ones we use small Latin letters. Furthermore, bold symbols refer to three-dimensional vectors and cursive ones the respective norms. The first two terms in Eq. (2) represent the kinetic energy of the ion and of the atoms, respectively, and V aa (r n −r j ) is a repulsive short-range potential with coupling constant g = 4πℏ 2 a aa /m. The potential V ai describes the atom-ion interaction, for which it is essential to retain the longrange tail (1). It is further characterized by the length R ? ¼ ð2m r C 4 =_ 2 Þ 1=2 and energy scales E ? ¼ _ 2 =½2m r ðR ? Þ 2 , where m r = mM/(m + M) is the reduced mass. For the 87 Rb/ 87 Rb + system one has R ⋆ ≃ 265.81 nm and E ⋆ ≃ k B × 79 nK (k B is the Boltzmann constant). Importantly, the separation of length scales is lacking as for typical atom density n = 10 14 cm −3 the mean interparticle distance n −1/3 ≃ 0.8R ⋆ is of the same order as the interaction range as well as the healing length (ξ ≃ R ⋆ ). At short range, the real interaction potential will deviate from the asymptotic formula (1). Here, we model the short-range details by the regularization 47 : Here, the b and c parameters have units of length and control the properties of the potential such as the number of bound states and their energies as well as the scattering length, while the longrange effects of the tail (1) remain accounted for. Crucially, the properties of the system depend not only on the energies, but also on the number of the available two-body bound states, as it will be demonstrated below. In most of the calculations we tune the potential in such a way that it has only one two-body bound state. Under this assumption, there is a unique connection between the b, c parameters and the s-wave scattering length a ai of the resulting potential. Here, f B and f I account for two-particle intra-species and interspecies correlations, respectively. These functions are constructed by matching the solution of the two-body scattering problem at short distances to an appropriate tail (see the "Methods" section), i.e. phononic decay in f B and mean-field prediction for a heavy ion in f I . Both functions contain variational parameters that are optimized by minimizing the expectation value of the Hamiltonian (2). VMC calculations provide the upper bound to the ground-state energy. In contrast the DMC approach aims to obtain the exact ground state energy of the system by solving the Schrödinger equation in imaginary time. We are interested in the regime of weak atom-atom interactions and fix the gas parameter to na 3 aa ¼ 10 À6 with a aa = 0.02R ⋆ .
Ground state energy. In Fig. 1 we illustrate the 'phase diagram' of the system that is characterized by two distinct sets of ground states: many-body bound state (MBBS) and a polaronic one. A spiral is used to illustrate that the energy E depends monotonously on the s-wave scattering length a ai , although different energies correspond to the same value of a ai depending on the number of bound states. While for a ai > 0 the potential (3) always has a bound state, for negative scattering lengths the atom-ion interaction can be tuned such that either a bound state is supported (left-bottom part of the helix, MBBS), or no bound state is present (left-upper part, polaronic). The 'attractive' polaron is typically encountered in ultracold quantum gases with neutral atomic impurities 4,5,39,40,50 . On the other hand, in the MBBS regime many bosons can be bound to the ion with a large binding energy. Importantly, the spatial range of the atom-ion interaction plays a crucial role in the formation of the MBBS, while for neutral impurities physics can typically be well described by assuming an effective short-range interaction. Figure 2 provides characteristic examples of the dependence of the system's total energy on the number of bosons in the MBBS (a) and polaronic (b) regime. In the MBBS case, we find that the absolute value of the energy grows almost linearly for a sufficiently small number of bosons. The dependence can be roughly approximated by the energy of N non-interacting particles bound to the ion, as shown with a solid black line in Fig. 2a. We also have verified that the effective impurity mass approaches the total mass of the MBBS, M ⋆ ≈ Nm. As the number of bosons is increased further, the energy starts to significantly deviate from the behavior (5) and it reaches a minimum at a critical number N c , which can be estimated from the extremum condition ∂E/∂N = 0. The value N c can be interpreted as the maximal number of bosons bound to the impurity, similarly to the analysis of the 1D case 22 . We note that N c could be defined in different ways, e.g. from the form of the atom-ion correlation function, g ai 2 ðrÞ, as N c ¼ n R ½g ai 2 ðrÞÀ 1 þ N c =N4πr 2 dr, or from the atom density far from the ion, i.e. n(1−N c /N). For the chosen parameters we obtain N c ≃ 140 almost irrespective of the exact value of the atom-ion scattering length, contrarily to the mean-field prediction of ref. 20 . This indicates that while the scattering physics of quasi-free bosons is determined by the scattering length, it is the large range of the Fig. 1 Schematic phase diagram. A schematic phase diagram showing various regimes for a single ion immersed in a dilute Bose gas upon changing the atom-ion scattering length and the number of the two-body bound states. The vertical axis indicates the change of the ion energy E as the value of atom-ion s-wave scattering length a ai is varied. The spiral refers to a continuous change of the energy leading to realization of different physical regimes as the value of a ai is cyclically changed from minus to plus infinity. Departing from the non-interacting case a ai = 0, with zero ion energy, E = 0, and by making the ion attraction stronger the following regimes are observed. A weak attraction leads to the polaronic regime where the ion can be described in terms of a quasiparticle. Such a description is no longer possible in the unitary regime marked by a diverging s-wave scattering length. The unitary point, |a ai | = ∞, is the threshold value for the formation of an ion-atom two-body bound state. For stronger atom-ion attraction many-body bound state is formed in which a large number of atoms is effectively trapped by the ion. potential that determines the number of bound particles. Note that our result is significantly larger as compared to the MBBSs for the case of a neutral impurity, for which Monte Carlo calculations predict only few atoms to be bound 39 . At the same time N c is much smaller than the number of bound atoms predicted by a mean-field-based estimate, suggesting that the effective gas parameter is significantly increased in the vicinity of the ion. For N > N c the energy increases, meaning that no more bosons are able to bind to the ion and the excess atoms start to form an almost uniform gas. This view is further corroborated by the snapshots of the system taken in Fig. 3 for different system sizes, where the green sphere of radius R ⋆ depicts the position of the ion (red symbols for bosons). Figure 3a shows the snapshot of the system in the case when the number of bosons is smaller than the critical one, N < N c , and therefore all bosons are close to the ion forming a spatially localized MBBS. Contrarily, Fig. 3b depicts the case with N > N c . In this scenario, the boson density around the ion is still higher than the average one and the excess bosons form a background gas. Figure 2b shows the energy dependence on N following the "polaronic" branch where no two-body bound states are present, for three characteristic values of the atom-ion scattering lengths: a ai = −0.1R ⋆ , a ai = −R ⋆ , and the unitary case a ai → ±∞. For large system sizes it is expected that the total energy can be decomposed into two contributions, the chemical potential μ pol of the ion and the energy of a homogenous gas, E = μ pol + Ngn/2. For sufficiently small |a ai |, the polaron energy can be calculated variationally 51 where a 0 ¼ 32 is the shift of the scattering resonance due to the bosonic ensemble. Therefore, in the regime of sufficiently weak interactions the energy is universal as it depends only on a ai and no finite range corrections are required.
It is important to notice that the Fröhlich model 29 alone is not sufficient to describe the ion quantitatively, predicting μ pol = −0.096E ⋆ for a ai = R ⋆ . Instead, beyond-Fröhlich perturbation theory correctly describes the polaron energy both in the weakly interacting regime [a ai = −0.1R ⋆ data in Fig. 2a] and remarkably even for strongly interacting polarons (a ai = −R ⋆ ). In the regime of weak interactions, the ion behaves similarly to a neutral impurity with short-range interactions, for which the VMC energy is shown with short-dashed lines. The unitary regime is reached when the atom-ion scattering length significantly exceeds the mean interparticle distance. An important feature of the ion impurity is that the energy of the many-body ground state is continuous when crossing the a ai → ∞ point and connects the polaron and MBBS which are both stable branches. We note that directly at unitarity the prediction (6), which is derived within Bogolyubov approximation, is beyond its validity and thus it is not shown in Fig. 2b. In particular, the Bogolyubov approximation becomes questionable close to unitarity, where the correlation functions shown in Fig. 4 indicate a varying local gas parameter similar to the discussion of beyond Bogolyubov corrections in ref. 52 . For the same reason, the Bogolyubov-Fröhlich description of the ion polaron 29 has to be significantly revisited.
Correlation functions. In order to analyze the spatial structure of the MBBS we further turn our attention to the atom-atom and atom-ion correlation functions. Typical examples are displayed in Fig. 4. As it can be seen, the atom-ion correlation features a pronounced peak indicating a strong bunching effect at distances where the atom-ion interaction potential is strongly attractive. Moreover, for N < N c (see red lines) the atom-atom correlation function does not approach a constant at long distances, but instead it decays exponentially, supporting the interpretation that essentially all bosons are bound to the ion and are localized at distances of the order of a few R ⋆ . The width of g ai 2 can be used as the definition of the size of the MBBS that can be interpreted as a mesoscopic molecular ion. For N > N c (see blue lines) the position of the peak does not change; the atom-atom correlation function, however, converges to a constant value which is slightly below unity. This demonstrates that the excess atoms are not bound to the ion and indeed form a bosonic background for the MBBS. The atom-atom correlation functions in the presence of the ion (dashed lines) also indicate the bunching behavior close to the ion. The effect is the strongest for small systems, N < N c , where the bosons tend to stay close to each other as they are a part of the MBBS [see also Fig. 3a]. This can be interpreted as an effective interaction within the medium induced by the impurity. As the system size is increased, g aa 2 ðrÞ starts to approach a constant value at large distances, i.e. the whole volume is filled with the gas [see also Fig. 3b], and the peak at short distances is correspondingly lowered. The asymptotic value g aa 2 ðrÞ ! 1 À N c =N reflects a smaller effective density, as N c atoms are bound to the ion. Eventually, in the thermodynamic limit atom-atom correlations will coincide with those of a homogeneous Bose gas without an ion (green line in Fig. 4). We have also found that for a ai < 0 the two branches have very different behaviors in terms of coherence, which is quantified by the quasi-particle residue Z ¼ lim r!1 g 1 ðrÞ corresponding to the long-range asymptotic of the residue function g 1 ðjr À r 0 jÞ ¼ hΨ y ðrÞΨðr 0 Þi where the field operator Ψ † (r) creates an ion at position r and 〈〉 denotes the ground-state average. Indeed, the residue is finite in the polaron branch and approaches unity (full coherence) in the limit of weak attraction, a ai → 0 − . Instead, in the MBBS branch the residue vanishes exponentially fast. Figure 5 shows typical examples of the decay of the ion residue function, g 1 (r), for a fixed number of particles for several choices of the atom-ion s-wave scattering length. We observe an exponential decay, which on the semi-logarithmic scale of Fig. 5 is seen as a linear dependence. This can be understood using a simple model. As discussed in the context of Fig. 2a, the energy of small clusters of bosons in the MBBS branch can be reasonably well interpreted in terms of N non-interacting atoms bound to the ion. The twobody scattering solution (8) for each atom-ion pair scales as f ðrÞ ¼ expðÀr=a ai Þ=r with r = |R−r i | for r ≫ R ⋆ . By assuming a product over r i , i = 1, ⋯ , N we arrive at the following approximate form for the residue function at large distances For deeply bound states, a ai → 0, all N particles are bound and participate in the MBBS. In turn, for weaker interactions or larger numbers of particles, the ion is able to capture only N c atoms. We take this effect into account by substituting N by N c in Eq. (7). As it can be seen in Fig. 5, the asymptotic expression (7) captures correctly the exponential loss of coherence. Let us finally briefly discuss the dynamic properties of the system. In the polaronic case, the ion effective mass approaches its bare value M ⋆ ≈ m in the limit of weak attractions, whereas for stronger ones it gradually increases for the given boson-boson scattering length to the value at unitarity M ⋆ ≈ 6m, which is substantially larger than the neutral impurity result M ⋆ ≈ 1.65m 39 . In the MBBS regime, M ⋆ becomes exceedingly large. In particular, for large E b and small N we find that M ⋆ ≈ N c m and the total energy is given by Eq. (5). We could not verify whether the relation M ⋆ ≈ N c m holds in the thermodynamic limit due to computational limitations.

Discussion
Our calculations are focused on the ground state properties of the system. Furthermore, we have assumed that the two-body ion-atom potential only supports one or zero bound states, while a realistic interatomic potential typically features hundreds of vibrational levels, similar to the potentials with van der Waals tails describing the interactions between neutral atoms. For the latter, however, the occupation of bound states of the interatomic potential is less likely for typical quantum gas densities, unless the system is tuned close to a Feshbach resonance, since the spatial range of the potential is on the order of a few nms, and therefore it can be well described by a pseudopotential. This is not the case for the atom-ion system, whose spatial range of the polarization potential is tens of times larger. A natural question is then the experimental relevance and the prospects for observing the phenomena we have described.
The main process stemming from the existence of deeply bound states is the three-body recombination, which will inevitably lead to losses, as it also does for neutral Bose polarons. The timescale for such losses can be estimated with the classical trajectory result for the three-body recombination rate constant, which can be expressed as K 3 ' 12:52 _ m r R ? ð Þ 4 E=E ? À Á À3=4 and the decay rate given by γ = K 3 n 253, 54 . While for a thermal gas with  ) for a ai = R ⋆ ; g 2 (r) in a weakly interacting Bose gas in absence of the ion is shown with a green line. The critical number of atoms N c which can be trapped by an ion is larger than the number of particles for N = 50, which explains the pronounced maximum observed in the red lines at distances of the order of the ion potential range, r ≈ R ⋆ . The data is shown up to the half-size of the simulation box L/2 and as a result the data with N = 50 abruptly stops without yet reaching a plateau. Instead, for N = 500 particles, the number of atoms is larger than the critical number N c and the half-size of the simulation box L/2 is larger than the size of the many-body bound state. As a result the long-range plateau observed for N = 500 atoms signals presence of a homogeneous gas of atoms. density n = 10 12 cm −3 and collision energies of the order of a milikelvin this gives lifetimes of the order of a second (with γ ≈ 2.4 Hz), an ion in a high density BEC is subject to much stronger losses (for n = 10 14 cm −3 at 1 μK γ ≈ 140 kHz). This leads to submilisecond time scales, which nevertheless are sufficient to observe ion dynamics in experiments 19 . For our gas parameter γ ≈ 600 Hz while the characteristic energy E ⋆ /ℏ = 1646 Hz. We further note that the quantum three-body recombination involving an ion is still not fully understood and may deviate from the classical result, e.g. it may feature minima for certain parameters (similar to loss recombination minima found for neutral atoms 55 ). In particular, the dependence on the binding energy of the weakly bound state should be similar to the case of van der Waals interactions for which K 3 / a 4 ai . For sufficiently small loss rates, experimental detection of the signatures of the MBBS formation can be realized e.g. by injecting the ion into a cold gas and dragging it slowly using an external electric field, as has been done in ref. 19 . The response of the impurity and the measured time of arrival at the detector will then be mainly determined by the dramatically increased effective mass of the impurity. Moreover, one can use precise in situ imaging techniques with high spatio-temporal resolution such as the setup based on charged particle optics 56 to study the increase in the gas correlation functions due to the presence of the ion which would provide further information about MBBS formation dynamics. Radiofrequency and microwave spectroscopy developed for neutral gases can be used here as well, in particular to investigate the polaronic branch. Finally, quenching protocols in which one makes use of the ion hyperfine structure can be implemented. Taking advantage of the existence of Feshbach resonances, one can transfer an initially noninteracting ion to a superposition state with vastly different scattering lengths and perform Ramsey spectroscopy 57,58 to determine e.g. the quasiparticle weight. We note that most of these techniques still require some experimental progress in reaching sufficiently low temperatures to increase the interaction times and the number of partial waves involved.
In conclusion, we have investigated the ground-state properties of an ion immersed in a dilute Bose gas by means of Quantum Monte Carlo and Bogolyubov techniques. We identify three physically different regimes in the many-body system depending on the presence of the bound state in the atom-ion scattering problem: (i) polaronic branch, two-body bound state is absent; (ii) many-body bound-state (MBBS) branch, two-body bound state is present; and (iii) unitarity, at the threshold of the appearance of the bound state. In the polaronic branch, manybody dressing leads to formation of a quasiparticle (ionic polaron). In the limit of weak interactions, variational methods developed for neutral atomic polarons accurately predict the energy of the system. Close to the unitarity limit the calculations unveil strong deviations from the approximate results. Finally, the MBBS branch is characterized by the formation of a large cluster (consisting of hundreds of atoms) around the ion, which in this case possesses a large effective mass, thus providing a strong analogy between the MBBS and a localized state. These quite distinct regimes should give rise to different timescales in the impurity dynamics observed in experiment, especially when combined with Feshbach resonances that allow for tuning the position of the last bound state 59 . Our results highlight the important role of the interatomic interactions which are strongly enhanced in the proximity of the ion, driving the system away from the weakly interacting regime to a nontrivial state characterized by the interplay of long-range interaction and high local density. Apart from the atomic gases, these findings can be relevant to condensed matter systems such as electrondoped exciton gases in heterostructures of two-dimensional semiconductors, where the long-range electron-exciton interaction also has long-range character which cannot be neglected for typical experimental parameters 60 .

Methods
Values of the parameters of the regularized potential. For the sake of numerical convenience we employed in our Monte Carlo simulations the regularized atom-ion potential (3) of the main text. We only considered a few specific values of the pair (b, c) that are characteristic of the three regimes outlined in the diagram of Fig. 1 of the main text: weak-coupling Bose polaron (WCP), MBBS, and strongcoupling Bose polaron (SCP). In Table 1 we list those values in units of R ⋆ and E ⋆ .
Trial wave functions for the Monte Carlo simulations. The trial wave functions are written as a pair product of Jastrow functions for both atom-atom and atom-ion correlations, featuring appropriate short and long-range asymptotic behavior [see Eq. (4) of the main text].
The short-range part of both the atom-atom and atom-ion Jastrow function is taken from the lowest energy solution of the two-body scattering problem where V r ai ðrÞ is the corresponding interaction potential of Eq. (3) of the main text and m r is the reduced mass. For the atom-atom wave function we choose scattering states with energy E = 0, whereas for the atom-ion wave function we use the exact two-body bound state with energy E b when a bound state is present.
The long-range (large distance) part of the Jastrow term is taken from hydrodynamic theory. As shown by Reatto and Chester in ref. 61 , if phonons are the lowest-energy excitations in the system, the long-range behavior of the manybody wave function can be factorized as a pair-product of Jastrow functions.
The atom-atom potential V aa in Eq. (2) of the main text is modeled by a repulsive soft-sphere potential: V aa (r) = V 0 > 0, for r < R ss with R ss = 0.1R ⋆ and zero elsewhere. The height V 0 is chosen to reproduce the desired value of the swave atom-atom scattering length a aa = 0.02R ⋆ ≪ R ⋆ . We further choose the density nðR ? Þ 3 ¼ 0:1288, resulting in gas parameter nða aa Þ 3 ¼ 10 À6 .
Mean-field estimate of the effective mass and critical number. In order to formulate a self-consistent mean-field theory in the ion's frame of reference, the following wave function can be used 27 Ψ G ðR; k; r 1 ; ; r N Þ / e ikÁR Y N n¼1 f ðr n À RÞe isðr n ÀRÞ : Here, k is the ion momentum, f 2 the relative probability distribution of the position of the ion and the bosons, while ∇ r s(r−R) indicates the fluid velocity relative to the ion. Performing functional variation of the expectation value of the Hamiltonian (2) of the main text, one obtains the ion effective mass 27,62 : Here, R 0 is a hard-core radius physically meaning the distance at which the atomion interaction starts to deviate from its long-range C 4 r 4 asymptote. Typically R 01 0a 0 with a 0 ≃ 53 pm being the Bohr radius. Furthermore, the distance R μ is defined as |V ai (R μ )| = μ with μ = gn the chemical potential of the bosons, from Parameters b and c of the regularized atom-ion interaction with corresponding threedimensional s-wave atom-ion scattering length, aai, at zero-energy and energy of bound state. NBS means no bound state is supported. The acronyms in the last column refer to: many-body bound state (MBBS); strong-coupling polaron (SCP); weak-coupling polaron (WCP). Close to unitarity the sensitivity of the b and c parameters is higher and therefore more digits are provided. The length and energy scales R ? ffiffiffiffiffiffiffiffiffiffiffi ffi mC 4 =_ 2 p and E ⋆ ≡ _ 4 /(m 2 C 4 ), respectively, are defined with respect to m r = m/2, i.e. equal atom and ion masses. which we get For the pair 87 Rb/ 87 Rb + with an atomic density n = 10 14 cm −3 we obtain R μ ≃ 1.2R ⋆ ≃ 6061a 0 . Given this, the formula (10)  Another estimate of N c can be attained via rate and Gross-Pitaevskii equations 20 . Denoting the binding energy of the two-body bound state as E b ¼ À_ 2 =ð2m r a 2 ai Þ, it can be shown that in the presence of many weakly interacting bosons E b ðN b Þ ¼ E b ½ma ai =ð6m r a aa N b Þ 2=3 . At thermal equilibrium one would expect that For the pair 87 Rb/ 87 Rb + with a ai = R ⋆ (i.e. E b ≡ E ⋆ ), a aa = 100 a 0 , and T = 10 nK (≪E ⋆ /k B ), we obtain a critical number of N c ≃ 372. This number is much smaller than the previous estimate (10), but also does not agree with our numerical simulations. Moreover, our study predicts that N c emerges already at zero temperature. Thus, it is not only determined by charge hopping and thermal fluctuations, but also by interaction-induced correlations. Finally, the formula (12) predicts a reliance on the atom-ion scattering length as a 4 ai , while our many-body analysis [see Fig. 2a] shows that there is almost no dependence on that length parameter. This finding highlights once more how semi-classical estimates can be quantitatively erroneous.

Data availability
The data that support the findings of this study are available from the corresponding author upon reasonable request.