Ionic polaron in a Bose-Einstein condensate

G. E. Astrakharchik, L. A. Peña Ardila, R. Schmidt, K. Jachymski, A. Negretti Department de F́ısica, Universitat Politécnica de Catalunya, Campus Nord B4-B5, E-08034, Barcelona, Spain Institut für Theoretische Physik, Leibniz Universität Hannover, Appelstr. 2, 30167 Hannover, Germany Max Planck Institute for Quantum Optics, 85748 Garching, Germany Munich Center for Quantum Science and Technology (MCQST), Schellingstraße 4, 80799 München, Germany Forschungszentrum Jülich Gmbh, Institute of Quantum Control (PGI-8), D-52425 Jülich, Germany and Zentrum für Optische Quantentechnologien, Fachbereich Physik, and The Hamburg Centre for Ultrafast Imaging, Universität Hamburg, Luruper Chaussee 149, 22761 Hamburg, Germany (Dated: April 29, 2021)

Introduction. An impurity immersed in a many-body quantum system constitutes a fundamental problem in condensed-matter 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 external potential 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][5][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 new 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 many-body systems. While this lack of scale sep-aration gives rise to a striking competition of few-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. The long range of the atom-ion potential, on the other hand, allows to investigate the formation of mesoscopic bound states [19][20][21] that would have dramatic impact on the ion transport dynamics.
Studies of the mobility of an ion moving in a bosonic medium date back to the early sixties with the aim of explaining the small ion mobility of ions in liquid helium [22][23][24][25][26]. Later, mean-field approaches were used to predict the formation of mesoscopic molecular ions [19], to estimate the number of excess atoms around an ion in a homogeneous Bose-Einstein condensate (BEC) [27], and path integral methods to determine the ion polaron properties in the strong-coupling regime [28]. A many-body study of a one-dimensional system consisting of a trapped ion in a bosonic environment, however, has elucidated the role of particle correlations in the ground state that the mean-field treatment cannot capture [21]. Motivated by this as well as by recent experiments [16,[29][30][31][32][33][34][35][36][37], 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 [38][39][40], bipolarons [41], as well as in 2D [42] and 1D polarons [43][44][45]. The method allows us to fully take into account the quantum many-body correlations that turn out to be crucial to predict how the competition between the few-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 a pseudopotential approach and lead to dramatically altered impurity properties.
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 by [46] The first two terms 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 long-range tail (1). The potential is characterised by the length R = (2m r C 4 / 2 ) 1/2 and energy scales , 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). The separation of length scales is lacking as for a typical atom density n = 10 14 cm −3 the mean interparticle distance (n −1/3 0.8 R ) is of the same order 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 long-range 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.
Methods. The simulations are performed by using variational (VMC) and diffusion Monte Carlo (DMC) methods. The VMC method samples the square of the trial wavefunction that we choose in the Bijl-Jastrow [48,49] form Here, f B and f I account for two-particle intra-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 [50], 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.02 R . Results. In Fig. 1 we illustrate the 'phase diagram' of the system that is characterized by two distinct set of ground states: many-body bound state (MBBS) and a polaronic one. As the "helix" structure indicates, for negative atom-ion scattering lengths both solutions are possible, while positive scattering lengths always lead to a MBBS. Indeed, 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,38,39,51]. 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 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. 2 (a). We also have verified that the effective impurity mass approaches the total mass of the MBBS, M ≈ N m. 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 [52]. The value N c can be interpreted as the number of bosons bound to the impurity similar to the analysis of the 1D case [21]. For the chosen parameters [50] we obtain N c 140 almost irrespective of the exact value of the atom-ion scattering length. This indicates that while the scattering physics of quasi-free bosons is determined by the scattering length, it is the large range of the potential that determines the number of bound states. Note that this number is significantly larger as compared to the many-body bound states for the case of a neutral impurity, for which MC predicts only few atoms to be bound [38]. At the same time N c is much smaller than the number of bound atoms predicted by mean-field estimate [50], 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). In Fig. 2(a) N < N c and therefore all bosons are close to the ion forming a spatially localised MBBS. Contrarily, panel (b) displays the boson density for N > N c . In this scenario, the density is still higher than average around the ion, while the excess bosons form a background gas around the MBBS. Figure 2(b) 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 + N gn/2. Total energy in units of E . For the lines, the mean-field Gross-Pitaevskii energy is shifted by the polaron energy, i.e. E(N ) = µ pol +gnN/2; solid line, no ion, µ pol = 0; dashed line, variational ansatz [53], µ pol is given by Eq. (6); short-dashed line, DMC calculation for short-ranged impurity-boson interactions from Ref. [38].
For sufficiently small |a ai |, the polaron energy can be calculated where a 0 = 32 3 √ π na 3 aa is the shift of the scattering resonance due to the bosonic ensemble [53]. 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.
We find that perturbation theory correctly describes the polaron energy both in the weakly-interacting regime [a ai = −0.1 R data in Fig. 2(a)] 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 which 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 unitary limit is continuous and connects the polaron and MBBS which are both stable branches. We note that at unitarity the prediction (6), which is derived within Bogolyubov approximation, is beyond its validity and thus it is not shown in Fig. 2(b). 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. [54]. In order to analyse the spatial structure of the manybody bound state 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 atomion 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 behaviour 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. 3(a)]. As the system size is increased, the peak at short distances is lowered and 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. 3(b)], 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 find 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. 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. Also, the dynamic properties of the system are completely different. 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 [38]. 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.
Conclusions. We 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 manybody 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; (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 this regime we verify that Bogolyubov theory, applicable also to neutral polarons, describes well the regime of weak atom-ion interactions, although it fails at the unitarity. Instead, 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 [55].

MEAN-FIELD ESTIMATE OF THE EFFECTIVE MASS AND CRITICAL NUMBER
Already in the sixties Eugene P. Gross has developed a self-consistent field theory in the ion's frame of reference using a wave function of the form [26] Ψ G (R, k; r 1 , . . . , r N ) ∝ e ik·R N n=1 f (r n − R)e is(rn−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. By performing functional variation of the expectation value of the Hamiltonian (2) of the main text, one obtains the ion effective mass [26,56]: Here, R 0 is a hard-core radius physically meaning the distance at which the atom-ion interaction starts to deviate from its long-range asymptote i.e. the potential is strongly repulsive for r < R 0 . A physically reasonable value is R 0 ∼ 10 a 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 which we get Thus, R µ is the distance at which the polarization potential becomes equal to the (mean-field) interaction energy of the bosonic ensemble. 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 (8)

VALUES OF THE PARAMETERS OF THE REGULARIZED POTENTIAL
As we anticipated in the main text, 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), many-body bound state (MBBS), and strong-coupling Bose polaron (SCP). In table I 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  TABLE I. Parameters b and c of the regularized atom-ion interaction with corresponding 3D s-wave atom-ion scattering length, aai, at zero-energy and energy of the uppermost energy level E b , i.e. the most loosely bound state. For M = +∞ the ion is infinitely heavy and the binding energy corresponds to the eigenenergy of a mobile atom in an external potential provided by the regularized atom-ion interaction. The binding energies denoted with an † are the binding energies of the polarisation potential (9) below, while that of the regularized potential (3) of the main text is zero, that is, no bound state is supported by the regularization (NBS means no bound state). 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 ≡ mC4/ 2 and E ≡ 4 /(m 2 C4), respectively, are defined with respect to mr = m/2, i.e. equal atom and ion masses.
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 it was shown by Reatto and Chester in Ref. [57], if phonons are the lowest-energy excitations in the system, the long-range behavior of the many-body wave function can be factorized as a pair-product of Jastrow functions.

Atom-atom wave function
The atom-atom potential V aa in Eq. (2) of the main text is modelled by a repulsive soft-sphere potential: V aa (r) = V 0 > 0, for r < R ss and zero elsewhere. The soft sphere diameter is fix to be much smaller than the typical range of the ion potential, R ss = 0.1R and the height of the soft-sphere potential V 0 is chosen to reproduce the desired value of the s-wave atom-atom scattering length a aa . The range of the atom-ion potential in turn is chosen to be large compared to the s-wave scattering length a aa = 0.02R R , but smaller than the mean inter-particle distance, and we choose n(R ) 3 = 0.1288. The resulting gas parameter of the bosonic ensemble corresponds to a very dilute gas, n(a aa ) 3 = 10 −6 .
The Jastrow term for atom-atom correlations is chosen by matching the solution of the two-body scattering problem (as determined by the phase shift δ) at short distances (smaller than a variational parameter R par ) and phononic decay [57] at large distances, where κ 2 = mV 0 / 2 − k 2 . The scattering momentum of k of the incident boson is considered to be a variational parameter and is typically of the order of 2π/L. The remaining parameters (C 1 , C 2 , C 3 , C4) are fixed by conditions of continuity of the Jastrow term and its first derivative.

Atom-ion wave function
The atom-ion interaction is modelled by the regularized potential, Eq. (3) of, the main text. At short distances the atom-ion Jastrow term for short distances from the atom-ion scattering solution, Eq. (11). At large distances, the ion-atom Jastrow term is constructed from perturbative theory. It can be shown [58,59] that a massive impurity induces a perturbation, δψ k , to the wave function of the bath in momentum space, δψ k = −gφ 0 /( in an exponential decay of the perturbation, δψ(r) ∝ (mgφ 0 / 2 ) exp(− √ 2r/ξ) with the typical decay distance given by the healing length ξ = /( √ 2mc). Motivated by this argument, the functional form for the atom-ion Jastrow term is chosen as The variational parameter A par has the physical meaning of an inverse healing length, while B par changes the weight of the ion defect, and R * is the matching distance. The constant A is fixed from continuity of the Jastrow term at r = L/2. This defines the boundary condition for the function and its derivative at R * . The Jastrow term at shorter distances is obtained from solving the atom-ion scattering problem with the atom-ion interaction potential.

EXPONENTIAL DECAY OF COHERENCE IN MBBS
The coherence properties of the ion differ greatly in the various interacting regimes. In the polaronic branch the residue is finite and approaches unity (maximal coherence) in the weakly interacting regime, similarly to neutral impurities. On the other hand, the coherence is greatly suppressed in the MBBS branch. Figure 6 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 semilogarithmic scale of Fig. 6 is seen as a linear dependence. The exponential decay can be understood in a simple model. As discussed in the context of Fig. 2(a) in the main text, 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 two-body scattering solution (11) 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, While for deeply bound states, a ai → 0, all N particles are bound and participate in the MBBS, for weaker interactions or large number 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. (14). As it can be seen in Fig. 6, the asymptotic expression (14) captures correctly the exponential loss of coherence.  Table I).

EXPERIMENTAL CONSIDERATIONS
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 polarisation potential is tens of times larger. Here, we comment on the question of the experimental relevance of the prospects for observing the phenomena we have described.
The main process stemming from 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 K 3 2.51[α 5 /m 2 /(k B T ) 3 ] 1/4 and the decay rate given by γ = K 3 n 2 [60,61]. Here, the static polarizability of the atom, α, is connected to the dispersion coefficient as: C 4 = αe 2 2 1 4π 0 with e being the electron charge and 0 the vacuum permittivity. While for thermal gases this gives lifetimes of the order of a second, an ion in a high density BEC is subject to much stronger losses. We 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 [62]). 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 the signature of the formation of the polaron and many-body bound state can be observed in experiment in which the ion is slowly dragged through the gas due to an external electric field. The response of the impurity will then be determined by the dramatically increased effective mass of the impurity, provided that the loss rate is kept reasonably small.
We also envision the possibility for measurements of dynamical observables inspired by Ramsey spectroscopy. Specifically, given the fact that the atom-ion scattering length should depend on the electronic spin configuration of the atom and the ion, one can prepare the ion in a superposition of two internal ion states (i.e., hyperfine states), whose scattering lengths are significantly different -exploiting Feshbach resonances. The different states in the superposition then act effectively as two branches of an interferometer and following protocols originally developed for neutral impurities [63,64] one can determine the Lohschmid echo of the compound atom-ion system; e.g. to determine the ion polaron quasiparticle weight. Alternatively, schemes based on Rydberg dressing could be utilized such that the effective atom-ion interaction can be controlled by external laser fields, whose coupling to the Rydberg state depends on the atom's internal state, as proposed theoretically in Refs. [65,66]. Here, either repulsive or attractive atom-ion