Artificial atoms from cold bosons in one dimension

We investigate the ground-state properties of weakly repulsive one-dimensional bosons in the presence of an attractive zero-range impurity potential. First, we derive mean-field solutions to the problem on a finite ring for the two asymptotic cases: (i) all bosons are bound to the impurity and (ii) all bosons are in a scattering state. Moreover, we derive the critical line that separates these regimes in the parameter space. In the thermodynamic limit, this critical line determines the maximum number of bosons that can be bound by the impurity potential, forming an artificial atom. Second, we validate the mean-field results using the flow equation approach and the multi-layer multi-configuration time-dependent Hartree method for atomic mixtures. While beyond-mean-field effects destroy long-range order in the Bose gas, the critical boson number is unaffected. Our findings are important for understanding such artificial atoms in low-density Bose gases with static and mobile impurities.


Introduction
A quantum state is bound if the probability to find parts of the system infinitely far from each other vanishes. It is one of the basic problems in quantum mechanics to determine conditions for a bound state to occur. Such problems are typically encountered in fewbody settings. However, they also play an important role in many-body physics. For example, a low-energy model of dilute many-body systems may include bound states as building blocks.
Only in some special cases, there exist results that provide conditions for binding. For example, any attractive potential supports a two-body bound state in one (1D) and two (2D) spatial dimensions, whereas solely 'deep' potentials can lead to a bound state in three dimensions (3D) [1,2]. For more than two particles, general conditions for binding are not known. Moreover, there seem to be no universal theoretical approaches to find them. Typically, one has to resort to numerical calculations, and only some problems can be addressed analytically (within certain approximation schemes). The latter class of problems includes for example the Efimov effect, which provides a universal mechanism for resonant interactions in 3D [3][4][5][6]. Another example of analytically tractable model are bound states of weakly repulsive bosons attracted by a short-range potential [7][8][9][10]. That system is reminiscent of an atom where the role of electrons is played by the bosons, and the nucleus is realized by the potential. Therefore, in what follows we shall occasionally refer to the system as an 'artificial atom from bosons'.
Artificial atoms were mainly studied in one or three spatial dimensions. In 3D, different theoretical methods seem to disagree on the number of bosons that can be bound to an impurity [7][8][9][10]. In 1D, there is a similar puzzle. The outcome of the mean-field approximation [7,11] is not supported by the phenomenological argument of Ref. [9]. The latter study argues that a dilute Bose gas can always be mapped onto a system of non-interacting fermions implying that only a single boson can be bound (cf. the Pauli exclusion principle). However, the mean-field studies demonstrate that the number of bound bosons can be large if boson-boson interactions are weak. These different results certainly motivate further investigations of the 'artificial atom' problem. Besides providing insight into conditions for binding, they shed light onto the physics of the Bose polaron (see, e.g., [12][13][14]), in particular onto the polaron-to-molecule transition region. In this paper, we focus on properties of a one-dimensional artificial atom.

Main results of the paper
Our first result concerns • the mean-field solution of an artificial-atom problem in a finite ring. This solution rigorously shows that an impurity can support a many-boson bound state.
Specifically, the mean-field solution unveils the existence of three different physical scenarios depending on the number of bosons, N , see Eqs. (4) and (5). Below a critical particle number all bosons are trapped by the impurity, which confirms the previous findings of Refs. [7,11]. In this regime, the density of the Bose gas decays exponentially, see Eq. (15). Therefore, we classify the system as a many-body bound state. At the critical particle number, all bosons are also bound to the impurity. However, the corresponding density decays as 1/x 2 , see Eq. (16). Therefore, we shall say that the system is in a critical state. If the number of bosons is larger than the critical one, a certain portion of bosons occupies scattering states, i.e., there is a significant probability to find a boson far from the impurity.
The second result of this paper is • a validation of the mean-field predictions for the artificial atom problem using numerical beyond-mean-field methods.
To this end, we use a recently introduced in-medium similarity renormalization group method for bosons (IM-SRG; also called flow equations) [15] whose accuracy is confirmed here using the well-established multi-layer multi-configuration time-dependent Hartree method for atomic mixtures (ML-MCTDHX) [16]. These methods allow us to study the decay of phase correlations and demonstrate phase coherence between the bosons in the vicinity of the impurity, see Fig. 4. We conclude that the mean-field solution describes the system well as long as all bosons are bound to the impurity. When bosons populate scattering states, they occupy the whole space, which lowers the density and increases phase fluctuations. The here employed numerical methods can be used to test the argument of Ref. [9] that bosons fermionize in artificial atoms in 1D. Our results suggest that fermionization occurs only for bosons in scattering states.

Structure of the paper
The paper is structured as follows: Section 2 introduces the system under consideration. Section 3 presents the mean-field solution, which is further analyzed in the zero-density limit in Sec. 4. Further, the mean-field solution is benchmarked against the flow equations results in Sec. 5. A mobile impurity in a Bose gas is studied in Sec. 6; it is concluded that the mean-field approach describes that system also well. Section 7 contains a brief summary and outlook. For convenience, we provide five appendices that elaborate on technical details of our study. Appendix A describes the employed numerical methods. They are benchmarked against one another in Appendix B. Appendix C presents a mean-field solution for a box trap. Appendix D contains some details on the mean-field solution in the zero-density limit. In Appendix E we discuss the smallest non-trivial system -a two-boson artificial atom.

Hamiltonian
We investigate a one-dimensional system of N bosons and an impurity in a ring. The standard Hamiltonian for such a system in the context of cold-atom experiments reads (see, e.g., [17,18] and references therein) where we assume = 1; y (x i ) is the position of the impurity (ith boson), m refers to the mass of the impurity atom, and M denotes the mass of a boson. For convenience, we shall use a system of units with M = 1 in the numerical calculations. To model atom-atom interactions, we employ delta-function potentials that describe s-wave scattering [19], which is dominant in the ultracold regime. Their strengths c and g can be virtually arbitrary thanks to the possibility to tune them via external fields using the phenomenon of Feshbach resonances [20].
For simplicity, we first focus on a heavy impurity, m/M → ∞; the role of the impurity mass is briefly discussed in Sec. 6. Without loss of generality, we place the impurity at y = 0 as illustrated in Fig. 1 a). Note that from the experimental point of view, a heavy impurity can be realized using atoms with very different masses, e.g., 7 Li (bosons) and 174 Yb (impurity) [21], such that the kinetic energy of the impurity can be neglected. Alternatively, a localized external field (light blade) can be used to trap the impurity atom (see, e.g., Ref. [12]) or to even simply produce a delta-function potential, see, e.g., Ref. [22]. Note that these different experimental methods might lead to different finite range effects whose investigation we leave for future studies.
Below, we focus on attractive boson-impurity (c < 0) and repulsive boson-boson interactions (g > 0). In the main part of the discussion, we consider periodic boundary conditions, i.e., particles are confined to a ring of length L. For completeness, we also present results for a box trap in Appendix B and Appendix C. Note that such boundary conditions can also be realized in experiment, see, e.g., Refs. [23,24].

Physical Picture
Before we proceed with an analysis of the Hamiltonian H, let us provide some basic insight into the physics of the system, which is driven by the interplay between attractive impurity-boson and repulsive boson-boson interactions. First we note that a deltafunction potential supports a single bound state (see, e.g., [25]). This means that an arbitrary number of bosons can be trapped by the impurity if the bosons are noninteracting (g = 0). In contrast, if 1/g = 0, the bosons fermionize [26], and only one boson can be trapped by the impurity. [Indeed, only one fermion can be trapped by the impurity due to the Pauli exclusion principle]. This observation implies that the interplay between the attractive impurity-boson interaction and repulsive boson-boson Sketch of the density of the Bose gas for a finite value of N , and a large system size (i.e., zero-density limit) assuming that g/|c| 1. Near the impurity, at distances ∼ 1/(M c), the density of the Bose gas is high due to the impurity-boson attraction and thus the effective boson-boson interaction is weak. For larger distances from the impurity, the density is low, which implies that the Bose gas is strongly interacting there. [Note that it is specific to 1D systems that a lower density corresponds to stronger interactions. For example, in three dimensions, the situation is reversed -low densities imply weak interactions.] interaction should lead to a critical number of bosons, N cr , that can be bound to an impurity.
In this work, we estimate this critical number N cr from a mean-field approximation (see also Refs. [7,11]). We also investigate the system using two numerical approaches, namely the IM-SRG [15] and the ML-MCTDHX [16] (for a brief description of these methods see Appendix A). These methods allow us to estimate the importance of beyond-mean-field effects from the decay of the quasi-long-range order as captured by the system's reduced density matrix, see Sec. 5 and Appendix B.
To understand why the mean-field approximation is accurate, let us consider the system in the limit L → ∞ [N is fixed] and g/|c| 1, which is one of the main limits of this work. The effective strength of the boson-boson interactions can be parameterized by M g/ρ(x), where ρ(x) is the density of the Bose gas. This parameterization is natural for 1D problems, see, e.g., [27]. The value of M g/ρ(x) is the smallest in the vicinity of the impurity, and it grows towards the edge of the bound state. For example for g = 0, ρ(x) = N |c|M e −2M |cx| , see, e.g., Ref. [25]. Assuming that this density approximates also the system with g/|c| 1, we conclude that Therefore, the mean-field ansatz must describe the Bose gas well in the vicinity of the impurity as long as x is not large. The characteristic width of this 'mean-field' region is proportional to 1/(M |c|). Farther away from the impurity, the density of the Bose gas is low, hence, the boson-boson interactions are strong, and the mean-field ansatz is no longer applicable, see Fig. 1 b). We extend this line of argumentation in Sec. 4.

Relevant Length Scales
Three parameters define the length scales in our model: 1/(M g), 1/(M |c|) and L.
One can always employ one of them to define the system of units. In our work, it is convenient to use 1/(M |c|) for this purpose because it is the only relevant length scale for the impurity-boson bound state if N = 1 or g = 0 (see Eq. (2)). The corresponding two dimensionless parameters are: the relative interaction strength α = c/g and the dimensionless length LM |c|. Note that the latter will be useful sometimes to express as LM |c|/N or equivalently as M |c|/ρ, where ρ = N/L is the density of the Bose gas without the impurity. This will be especially convenient in Sec. 5 where we consider how the system approaches the zero-density limit (L → ∞ and N is finite).

Mean-Field Solution for the Heavy Impurity Problem
For a system consisting of weakly interacting bosons it is reasonable to assume that the ground-state wave function is a product state: Here, f (x) is a singleparticle function obtained by minimizing Φ|H|Φ . This minimization procedure leads to the Gross-Pitaevskii equation (GPE) [27]: where µ is the chemical potential †. By assumption, the function f is periodic i.e. f (−L/2) = f (L/2). [For a brief discussion of a system in a box trap where f (−L/2) = f (L/2) = 0, see Appendix C]. Note that some care is needed when using a meanfield approximation in 1D where quantum fluctuations destroy the condensate in the thermodynamic limit [27,[29][30][31]. We shall rely on ab-initio numerical calculations to confirm that the mean-field approximation is indeed accurate, at least for describing the Bose gas in the vicinity of the impurity. The relevant physical picture is given in Sec. 2.2. We notice that Eq. (3) with c = 0 is integrable, see, e.g., Ref. [32]. Therefore, one can follow the same strategy as when solving the Schrödinger equation with the delta-function interaction [25], i.e., use the known solutions for x > 0 and x < 0, then implement the boson-impurity interaction cδ(x) as the boundary condition at x = 0. Once the mean-field solution is obtained, it is possible to calculate any observable of interest. For example, the energy of the system is determined by Below, we present the two solutions to Eq. (3) that, by increasing L to infinity, connect adiabatically to the two different physical situations: (i) all bosons are bound †It is interesting to note that this equation was also derived and studied for a heavy atom in a strong magnetic field, see parameter regime 5 ('region 5') of Ref. [28].
to the impurity, (ii) no boson is bound, see the next section. The two solutions coincide at the threshold for binding, which we refer to as the point of transition (PoT). Note that the bound-state solution was discussed in Refs. [7,11] for L → ∞.
We focus on systems with finite values of L, since they allow us to directly benchmark the mean-field method against beyond-mean-field numerical approaches. In addition, our solution is relevant to cold-atom experiments, which typically have a finite size. Last but not least, the finite-L solution provides insight into the case with N > N cr , which is important for understanding the thermodynamic limit, as we plan to discuss in an upcoming work.

Mean-field Solutions
The first solution to Eq. (3) reads It is determined by the Jacobi elliptic function ds [32]. By construction, this solution is parity symmetric f mbb (−x) = f mbb (x). The chemical potential is where K is the complete elliptic integral of the first kind. The second solution to Eq. (3) is given by the Jacobi elliptic function ns: The corresponding chemical potential is The parameters p scatt ∈ [0, 1), p mbb ∈ [0, 1) and δ ‡ that enter in the definitions above are fixed by normalization, and the boundary condition due to the impurityboson potential It is straightforward to write these conditions in a more explicit form. For example, for f mbb , the normalization condition leads to where K = K(p), and we imply that p = p mbb [E, E(p) (not to be confused with the energy), sc and nd are standard Jacobi functions [32]]. The boundary condition can be written as Equations (6) and (7) can be satisfied only for N ≤ N cr , i.e., f mbb can describe only such systems. The solution f scatt describes systems with N ≥ N cr . The calculation of N cr will be given in the next section, see Eq. (10). The subscripts 'mbb' (many-body bound) and 'scatt' (scattering) are motivated by the observation that for a large system (L → ∞) N cr is the maximal number of bosons that can be trapped by the impurity, see Ref. [7,11] and the discussion in the next section. Note also that the chemical potential for the first (second) solution is negative (positive) for large system sizes since p mbb → 1 and p scatt → 1 for L → ∞. This is another indication that the first solution describes a many-body bound state while the second is applicable if the bosons occupy scattering states.
We illustrate mean-field solutions in Fig. 2 for different values of N and L. At the position of the impurity, any solution f reaches its maximum as a result of the attractive impurity-boson interaction. Increasing the number of particles decreases the binding energy per particle (increases the energy per particle) due to the repulsive bosonboson interaction, see Fig. 2 d), which also leads to a more flat profile of the density for the largest considered systems. The insets in panels a) and b) show the parameters p and δ as a function of the particle number. As it can be seen δ → 1 for increasing particle number. The parameter p first drops down to 0 at the critical particle number N cr = 11, see Eq. (10), and then rises again towards 1 [Note that for N ≤ 11 (N > 11), p corresponds to p mbb (p scatt ).]. For the largest ring size p is larger; for L → ∞ (not shown) we find empirically that p → 1 except in the vicinity of N cr . The chemical potential becomes negative for the largest ring size and N 11 (Fig. 2 c)), which is an indication of a bound state. We discuss this behavior in more detail in the following sections.

Point of Transition
The point of transition from one solution to another occurs at p mbb = p scatt = 0 (cf. the insets in Figs. 2 a) and b)). In this case, the functions in Eqs. (4) and (5) coincide: The corresponding chemical potential reads is for LM |c| = 5. The relative interaction strength is fixed to α = −5. Therefore the maximal possible particle number in a bound state is N cr = 11, see Eq. (10). The insets show the parameters p and δ as a function of the particle number. The lower panels depict the corresponding chemical potential c) and the energy per particle d) for different system sizes (see legends).
It vanishes for large system sizes, i.e., µ PoT = 0 for L → ∞. Normalization, and the boundary condition due to the delta-function potential determine the parameter δ and the critical number of bosons Note that N cr → ∞ when g → 0, and N cr = 1 when g → ∞, in agreement with our discussion in Sec. 2.2. Equation (9) shows that δ is determined only by the dimensionless parameter M |c|L, whereas N cr depends only on the ratio |c|/g. This decoupling of δ and N cr is unexpected for systems with finite values of L. It suggests scale invariance of the problem at the point of transition, and leads to a number of surprising consequences.
In particular, at the point of transition, the energy of the system also does not depend on L: which implies a state of zero pressure, in a sense that it costs no energy to adiabatically change the radius of the ring. This unique signature will later be used to identify the transition point in our numerical simulations. Note that Eq. (11) provides a variational upper bound on the exact value of the energy. It rigorously shows that more than a single boson can be trapped by the impurity, since this upper bound is below the ground-state energy of a single boson, −M c 2 /2, for N cr ≥ 3.

Zero-Density Limit: Mean-Field Results
In this section, we discuss the limit of vanishing density (ρ → 0) that occurs for a fixed number of bosons in a large system, i.e., L → ∞ (see also Refs. [7,11] and Appendix D). This limit provides insight into the general results of the previous section.

Many-body bound state
The limit ρ → 0 leads to p mbb → 1 (cf. Fig. 2 and its discussion), so that Eq. (4) can be written as where The quantity x mbb sets the characteristic width of the state §. If we define ζ = 0 for N > N cr , ζ can be seen as the 'order' parameter for our system. Indeed, ζ is positive for a many-body bound state, and vanishes as we approach the point of transition. The respective chemical potential reads It is negative which means that adding an additional boson lowers the total energythis is a typical characteristic of a bound state. The energy of the system is given by Additionally, for large values of |x| Eq. (12) yields §Note that x mbb is proportional to the characteristic size of a one-body bound state 1/M |c|. However, it can be much larger since x mbb grows with N . which corresponds to a typical tail of a bound state whose extension is defined by x mbb (see, e.g., Ref. [1]). Note that Eq. (15) is valid only for N < N cr . At the point of transition (ζ → 0, x mbb → ∞), another function will describe the tail of the state, see below.

Point of Transition
At ζ = 0, the mean-field solution of Eq. (12) can be further simplified We see that for our many-body problem there is a finite probability to find a boson next to the impurity even at the threshold of binding. This clearly distinguishes the manybody problem from the one-body system [see Eq. (2)] where this probability vanishes in the limit L → ∞.
The characteristic length x mbb diverges, and we need another quantity to describe the size of the state. It cannot be a root-mean-square radius, because of the 1/x tail of f PoT (x). Still, we can define a meaningful size of the state as This quantity defines the spatial region which contains half of the probability density, i.e., PoT (x)dx = 1/2. Note that x PoT is given by the size of a oneparticle bound state [Eq. (2)], which supports the physical picture provided in Sec. 2.2.
The energy of the system was already given in Eq. (11), which is independent of the system size L. The chemical potential is zero. This implies that if we add more bosons, they must occupy scattering states. Hence, Eq. (11) defines the energy for all systems with N ≥ N cr and ρ → 0.

Scattering state
The function f scatt of Eq. (5) in the limit ρ → 0 (p scatt → 1) can be approximated as follows For this solution, it is not possible to fulfil simultaneously normalization, and the delta-potential boundary conditions. If we impose the latter condition, we derive a nonnormalizable wave function for L → ∞. If we demand a normalized state, the resulting function is constant in space and does not include the impurity potential. We interpret this result as if all bosons occupy scattering states and distribute over the whole system until they are no longer affected by the impurity. This state is physical, however it is not the ground state, in which N cr bosons are bound to the impurity, and other bosons occupy scattering states. The GPE cannot describe this physics because it assumes that all bosons occupy the same orbital. Hence, as soon as there are more bosons in the system than can be supported by the many-body bound state, all occupy scattering states within the mean-field approximation.
An appropriate variational ansatz for large (but finite) L should include two parts, where the first part describes a many-body bound state, and the second one accounts for the bosons that occupy scattering states. In the low-density limit this leads to a Tonks-Girardeau gas formed outside the many-body bound state. We leave an investigation of such an ansatz for a future study.

IM-SRG Results: Approach to the Zero-Density Limit
To test our findings from the previous sections, we use the IM-SRG and ML-MCTDHX methods, which are briefly discussed in Appendix A, see the references given there for more details. Both numerical methods are able to capture corrections stemming from quantum fluctuations, and agree for the considered parameters. Therefore, in this section we illustrate our numerical results only for IM-SRG, see Appendix B for some ML-MCTDHX results. We focus on the question of approaching the limit L → ∞. This allows us not only to test the mean-field predictions but to also address beyond-meanfield corrections that must be important far from the impurity, see Fig. 1 b).

Energies
In Fig. 3, we present the energy per particle as a function of the inverse density (1/ρ) for different values of N . For the considered parameters, the critical number of bosons that can be trapped by the impurity is N cr = 11. We see that the IM-SRG data agree with the mean-field results well. Only small deviations are visible for N N cr . We attribute these deviations to residual beyond-mean-field effects naturally captured by the IM-SRG. The energies of the systems with N ≥ N cr decrease and in the limit of L → ∞, we expect them to approach the critical energy of Eq. (11). Unfortunately, we cannot follow this convergence for larger values of L; we observed that the IM-SRG method is not accurate for M |c|/ρ ≥ 1. In particular, for the largest considered particle numbers, the truncated flow equations diverge. This can be interpreted as a sign that the system becomes progressively more correlated, and IM-SRG cannot map the reference state (in our case a condensate) onto the real ground state of the system.
For N < N cr , Fig. 3 shows that the energy increases with the size of the ring. This is a typical behavior for bound states (at least for L → ∞), where the potential energy dominates the kinetic one. For N = 1, this increase can be understood using the equation for the binding energy  Figure 3: The energy per particle as a function of the inverse density M |c|/ρ for different particle numbers. The circles are calculated using IM-SRG. The dashed-dotted curves represent our mean-field results. The solid (horizontal) lines correspond to the meanfield prediction for ρ = 0, see Eq. (14). The left panel shows the data for N ≤ N cr . The right panel is for N ≥ N cr . We fix α = −5 which leads to N cr = 11. Notice that the numerical error bars on the IM-SRG data calculated according to Appendix A are smaller than the sizes of the markers.
For N = N cr the energy remains nearly constant with respect to the system size as predicted by Eq. (11). There is a very weak dependence on L pointing to beyond-mean-field effects. For N > N cr , the energy is a decreasing function of L. Our interpretation is that some bosons are now dropped out of the many-body bound state. Their kinetic energy decreases approximately as 1/L 2 , allowing us to conjecture the following behavior of the energy in the limit L → ∞ where N = (N − N cr − 1)/2 for odd values of N − N cr and with N = (N − N cr )/2 for even values of N − N cr . These expressions are the sums of the energy of the many-body bound state and the energy of the Tonks-Girardeau gas made of N − N cr particles, assuming that there is no interaction between the bound state and bosons in scattering states ¶. All in all, the IM-SRG data support the existence of different physical scenarios that correspond to bound and scattering states. However, note that our numerical analysis ¶Equations (17) and (18) assume that the size of the ring is sufficiently large in the sense that 1/k F 1/|c|M , where k F is the Fermi wavelength corresponding to the Tonks-Girardeau gas and 1/|c|M is the characteristic width of the many-body bound state. Assuming that k F = πρ, we derive the condition M |c|/ρ π. This condition implies that the ring sizes used in Fig. 3 are too small to numerically confirm Eqs. (17) and (18). cannot rule out the possibility that N cr becomes larger when L → ∞. In particular, we cannot rule out bound states with an infinite number of particles that are exponentially weakly bound in the limit L → ∞. However, one does not expect this to happen because far from the impurity the bosons interact strongly (fermionize).

Densities and Phase Fluctuations
Here, we calculate the density of the Bose gas in the ground state, Φ gr . We also investigate beyond-mean-field effects. To this end, we estimate phase fluctuations (also known as phase correlations), δΦ xx , from the one-body density matrix according to the prescription (see, e.g., [27,33,34]) The quantity δΦ xx is a measure of the off-diagonal quasi-long-range order, which vanishes for a condensate (mean-field) state. Note that δΦ xx is not only a convenient theoretical object for studying the importance of the beyond-mean-field effects. It also leads to experimental indicators of phase coherence that are observable through Bragg spectroscopy and interferometry, see, e.g., Ref. [35]. In Fig. 4, we show ρ(x) and δΦ xx for L = 0.1N/M |c| and L = 0.5N/M |c|. For all considered parameters, the IM-SRG and mean-field results agree on the density profile of the Bose gas. The density is the highest in the vicinity of the impurity, as expected. For the largest values of N , it features a weak dependence on N irrespective of the (considered) ring size. In spite of this, there is a noticeable increase of beyond-meanfield correlations as identified by the non-vanishing phase fluctuations. Their effect is more pronounced for the largest ring and N ≥ N cr , especially in the region with low densities. This observation is in agreement with the physical picture outlined in Fig. 1: low densities lead to strong boson-boson interactions, which can be quantified by M g/ρ(x). Surprisingly, IM-SRG and mean-field results are in a reasonable agreement even when M g/ρ(x) is of the order of unity, where the mean-field treatment is not expected to be valid. It is also worthwhile noting that the mean-field approximation is valid, in particular Eq. (10), even for the smallest non-trivial bound system -a twoboson artificial atom, see Appendix E.
Note that phase fluctuations are the strongest for the largest considered N . This can be rationalized in the following way. For N = 15, a few bosons are not trapped by the impurity. Therefore, the probability of strong boson-boson interactions far away from the impurity is high leading to large phase fluctuations. In contrast, for small particle numbers (e.g., N = 5), phase fluctuations may be small even if the density is low. The probability of finding two bosons outside the many-body bound state in this case is exponentially suppressed.

Mobile Impurity
A single mobile impurity atom in a weakly interacting Bose gas is an experimentally relevant system [12,36], which motivated various theoretical studies of a 'Bose-polaron', see, for example, Refs. [37][38][39][40][41][42][43]. Here, we complement those studies by considering the many-body bound state that follows from our results in the previous sections. The dashed-dotted curves represent the mean-field results. The solid (horizontal) lines correspond to the mean-field prediction for ρ = 0, see Eq. (14). The left (right) panel shows the energy for N ≤ N cr (N ≥ N cr ). We fix α = −5 which leads to N cr = 11. Notice that the numerical error bars calculated within IM-SRG according to Appendix A are smaller than the sizes of the markers.

Mean-field analysis
To investigate an impurity with a finite mass, we use the mean-field ansatz in the frame 'co-moving' with the impurity [39][40][41][44][45][46]. This frame is introduced via the set of new coordinates where θ(x) is the Heaviside step function. These coordinates allow us to exclude the position of the impurity from the Hamiltonian (similarly to the Lee-Low-Pines transformation [47]) . In the new coordinates, the Hamiltonian reads as follows where P is a quantum number -the total (angular) momentum of the system. For simplicity, we consider the case P = 0, which corresponds to the ground-state manifold. The Gross-Pitaevskii equation that follows from Eq. (22) reads (see, e.g., [40]): where κ = mM/(m + M ) is the reduced mass. This equation is equivalent to Eq. (3) up to a change of the mass of the boson M to κ. In this sense, all of our mean-field results from Secs. 3 and 4 also apply to a mobile impurity.
Transformation to the 'co-moving' frame allows us to use the analytical results of the previous sections. The mean-field approximation in the laboratory frame will lead to a system of coupled Gross-Pitaevskii equations, which one should solve numerically, see Appendix B.

IM-SRG results
We use IM-SRG to validate the mean-field predictions of Eq. (23). We focus on the case with m = M . In Fig. 5, we show the ground-state energy as a function of 1/ρ -similar to Fig. 3. First of all, we see that the mean-field and IM-SRG results are in agreement. Furthermore, just like before, the energy of the system with N < N cr increases as a function of L. For the critical number of bosons, the energy remains nearly constant. Note that according to Eq. (10) the critical number of bosons does not depend on the mass of the impurity. Our numerical simulations confirm this result. Finally, we used IM-SRG to calculate the density and phase fluctuations of the Bose gas in the presence of a mobile impurity. The comparison of the mean-field predictions to the IM-SRG is similar to the one presented in Fig. 4. Therefore, we refrain from discussing it further.

Summary & Outlook
We studied a one-dimensional artificial atom made of bosons. First, we analyzed this system within the mean-field approximation, and presented two possible solutions. In the limit L → ∞, the solutions correspond to two different physical scenarios with the bosons bound (or not) to the impurity. The critical state in between these scenarios is a zero-pressure state, meaning that its energy does not depend on the radius of the ring. We presented analytical expressions that describe this state.
Second, we investigated the system numerically using beyond-mean-field methods (IM-SRG and ML-MCTDHX). Our numerical simulations justified the use of the meanfield approximation for studying artificial atoms from bosons in one dimension. They confirmed the existence of bound, critical and scattering states in the system. Despite the validity of the obtained mean-field solutions, we argued that quantum fluctuations are present in the tail of the wavefunctions. Therefore, only the bosons near the impurity are described with a mean-field ansatz well. Bosons far away from the impurity are strongly interacting, supporting the phenomenological argument of Ref. [9]. However, their influence on the system can be neglected for particle numbers smaller than the critical one, because the attraction from the impurity assures a sufficiently large region with high density where particles are weakly interacting * * . Although, we mainly focused on a heavy impurity, we also showed that our results are applicable for a mobile one.
Further studies are needed to understand Bose systems with N cr + 1 particles in the limit L → ∞. Our results indicate that the mean-field approach is not suitable for such studies. In particular, it cannot be used to calculate the effective boson-artificial-atom interactions. The knowledge of this interaction will simplify the analysis of low-density Bose gases with attractive impurity-boson interactions.
Our results pave the way for investigations of many-atom physics using artificial atoms as elementary building blocks. For example, a lattice of heavy impurities * * Note that in Appendix E we validate the mean-field solution even for a two-boson system where the strength of the boson-boson repulsion can be of the same order of magnitude as the one of the boson-impurity attraction.
immersed in a Bose gas may feature different phases (e.g., Mott insulator and superfluid) depending on the strength of the boson-impurity and boson-boson interactions. Dilute systems of artificial atoms based upon mobile impurities can enjoy the physics of cold gases. To explore that context, one needs to understand the effective interaction between two artificial atoms.

Appendix A. Numerical Methods
In this Appendix, we briefly discuss the two numerical methods used in this work. The first method is called the flow equation approach or IM-SRG ("in-medium similarity renormalization group"). Our numerical implementation of this approach is based upon previous works [40,48] (see also Ref. [15] for a study of the Lieb-Liniger gas), which are inspired by the methods known in condensed matter and nuclear physics (see e.g. [49][50][51]). The second method is called the multi-layer multi-configuration time-dependent Hartree method for atomic mixtures (ML-MCTDHX) [16] (see also a relevant review on the topic [52]). It is a variational approach that has been extensively used, among others, for studying systems with impurities [42,43,[53][54][55][56].

Appendix A.1. Flow Equation Approach (IM-SRG)
The flow equation approach (block)-diagonalizes the Hamiltonian in second quantization, Here, s is the flow parameter, which formally plays a role of (imaginary) time. The generator of the flow η has to be chosen such that the off-diagonal matrix elements vanish in the limit s → ∞ [49].
In this work, we aim to decouple the ground state from the rest of the Hilbert space. Therefore, we normal order the Hamiltonian using a reference state following the prescription in Ref. [15]. This leads to the normal-ordered Hamiltonian where we denote normal ordered operators with : O :. The matrix elements f ij and Γ ijkl describe one-and two-particle excitations from the reference state. For the generator, we use these are the matrix elements which need to vanish in order to decouple the ground state from the excitations. Therefore, once the flow equation converges, our ground state is decoupled. The transformation governed by the flow equations can also be understood as a mapping of the reference state onto the real ground state of the system. Since we are interested in a system of bosons, it is reasonable to use condensate as reference state. Our reference state is constructed iteratively: Starting from the ground state solution of the non-interacting Hamiltonian, the density is calculated and used as the new reference state. This procedure is repeated until the density converges. Note that also other choices for the reference state are possible such as the mean-field solution, see Ref. [48]. For our system of interest such a reference state leads to the same result.
Induced higher order terms make it impossible to solve Eq. (A.2) exactly, and should be truncated. In our truncation scheme, we truncate at the two-body level, while keeping three-body operators which contain at least one a † 0 a 0 operator. This leaves us only with zero-, one-and two-body operators in Eq. (A.2) which leads to a system of coupled, closed, non-linear differential equations, which we solve numerically [15,48]. We estimate the error due to the neglected pieces (called W ) using second order perturbation theory where Φ p is a state that contains three-body excitations and Φ ref is our reference state.
To construct the Hamiltonian in second quantization we use the solution of the one-body Hamiltonian of our system. Since we can only work with a finite Hilbert space, we solve the flow equations for different numbers of basis states (in our case n ∈ [11, 13,15,17,19,21]). For the energy, we fit these values with to obtain the result in infinite Hilbert space. For other observables, such a fit is not always possible. In such cases, we take the result for the largest Hilbert space as our result and estimate the error by the largest deviation between the results for the different numbers of basis states. So there are in total two contributions to our error bars: The truncation error from neglecting higher order terms in the flow equation and the truncation error due to a finite Hilbert space. For a more detailed description of the method we refer to Ref.
[15], where the flow equations and our estimate of the truncation error are introduced, see also Ref. [48] for information about calculation of observables and a detailed explanation of our estimate of error bars.

Appendix A.2. ML-MCTDHX approach
In the ML-MCTDHX approach, the Hilbert space is truncated in a variationally optimal manner. To this end, one employs a time-dependent moving basis in which the system is instantaneously optimally represented through time-dependent permanents † †. In this sense, the many-body wave function is expressed with respect to bosonic number states | n ≡ |n 1 , n 2 , . . . , n D ; t and time-dependent expansion coefficients C n (t) as follows |Ψ(t) = n C n (t) |n 1 , n 2 , . . . , n D ; t . (A.7) Here, | n; t built upon time-dependent single-particle functions ϕ i (t) with i = 1, 2, . . . , D. The summation in Eq. (A.7) is performed over all possible combinations n i such that the total number of bosons N is conserved. In our numerical implementation, the single-particle functions ϕ i (t) are expanded in a time-independent primitive basis of dimension M ‡ ‡ that is based upon a sine discrete variable representation for the box potential with hard-wall boundary conditions at ±L/2. To calculate the ground-state wave function of the many-body setting, we determine the underlying equations-ofmotion for the coefficients C n (t) and the single-particle functions ϕ i (t) following the Dirac-Frenkel [57] variational principle. An imaginary time propagation method is used to obtain the system's ground state configuration. More details on the ingredients of this variational method can be found in Refs. [16,58].

Appendix B. ML-MCTDHX results
In the main text, the analytical solution for the bound state has been benchmarked against IM-SRG data. We have checked that these results are in agreement with the predictions of the well-established ML-MCDTHX approach. This is illustrated in Fig. B1 where the densities and phase fluctuations for the largest value of N considered in the main text are shown (cf. Fig. 4). Below, we study a system in a box potential, thus, exploring the formation of the artificial atom from bosons in the presence of hard-wall boundary conditions. This allows us to further understand the validity of the relatively novel IM-SRG method. Afterwards, we discuss the mean-field approximation to a mobile impurity in a Bose gas without the transformation to relative coordinates, Eq. (21). † †For a multicomponent setting, the variational ansatz has a multilayer structure allowing one to include both intra-and interspecies correlations, see Ref. [16]. Here, we describe a reduction of ML-MCTDHX to a single-component system that is investigated.
‡ ‡In the limit where D = M the wave function expansion of Eq. (A.7) is equivalent to a full configuration interaction approach, while for D = 1 it reduces to a single product state, which automatically satisfies symmetrisation conditions for bosons, and thus corresponds to the mean-field approximation. We show in Fig. B2 the energy, density and phase fluctuations of the Bose gas for different particle numbers N . The used values of N correspond to bound, critical, and scattering states discussed in the main text. Note, however, that the box trap modifies all properties of the system if L is of the order of 1/M |c|. For example, we noticed that we need to use stronger impurity-boson interactions (and therefore larger numbers of particles) than in the main text to be able to observe significant beyond-mean-field effects.
The ML-MCTDHX and IM-SRG results for the energy (panel a)) and the density (panel b)) are in agreement. However, phase fluctuations (panel c)) show some deviations for larger particle numbers. We notice, that while the density and the energy are accurate already for a small number of orbitals in ML-MCDTHX, phase fluctuations require more involved simulations. This is expected for several reasons. In particular, phase fluctuations require to determine the off-diagonal of the reduced density matrix which is a higher order observable. Note that ML-MCTDHX contains in general more information about the Hilbert space of the system in comparison to IM-SRG. Furthermore, ML-MCTDHX provides a direct access to spatially resolved observables and multicomponent settings. In that light, ML-MCTDHX calculations of certain observables are computationally more demanding than those with IM-SRG.
Nevertheless, increasing the number of orbitals leads to an agreement between the IM-SRG and the ML-MCDTHX results also for the phase fluctuations. We showcase this statement in panel d), presenting the phase fluctuations within ML-MCDTHX for an increasing orbital number in the case of N = 26 (a similar pattern is expected for N = 40). We observe a systematic convergence behavior. The main disagreement is near the boundaries of the box trap where the calculation of phase fluctuations becomes hard due to almost zero densities, see Eq. (20) especially so for ML-MCTDHX which operates in first quantization. We conclude that the decrease of phase fluctuations near the boundary is a numerical artifact caused in part by the presence of hard walls. Thus, we only show results for x < 0.35L. Overall, both numerical methods predict the same behavior for the observables of interest. For the sake of completeness, here, we apply the mean-field approach to the problem of a mobile impurity without the transformation to relative coordinates, see Eq. (21). We numerically solve the following set of coupled Gross-Pitaevskii equations for increasing particle numbers while fixing the ratio N/α = −1. Here, Ψ I (Ψ B ) is the mean-field wave function of the impurity (bosons). Below, we assume m/M = 1.
To justify this mean-field approximation, we benchmark it against ML-MCDTHX. [Note that we cannot use the current implementation of IM-SRG for such a benchmark, as it cannot be used to study multicomponent systems.] Our findings are illustrated in Fig. B3 where the one-body densities of the impurity and the bosons are shown. It becomes evident that a larger particle number results in a higher bosonic density at the position of the impurity as the effective boson-boson interaction decreases if the ratio N/α is kept fixed.
The considered weak boson-boson interactions lead to a good agreement between the mean-field and ML-MCTDHX methods, at least for the density of the Bose gas. This is expected since for these weak interactions boson-boson correlations are suppressed. For the impurity, the deviations between the ML-MCDTHX and the meanfield predictions become more noticeable for the largest numbers of bosons. Particularly, the impurity appears to be more spatially localized in the mean-field approach.
Our data allow us to conclude that the mean-field approximation is able to provide adequate results also without transformation to a co-moving frame. However, such a transformation is needed to obtain some analytic insight into the system, as we discuss in the main text. If one is simply interested in estimating lower order observables such as densities in the mean-field approximation, then it seems that it is sufficient to work in the laboratory frame.

Appendix C. Mean-field Solution for Hard-Wall Boundary Conditions
To complement the mean-field studies in the main text, here, we present a solution of Eq. (3) for a box trap, i.e., for f (−L/2) = f (L/2) = 0. For these boundary conditions, the solution that becomes the many-body bound state in the limit L → ∞ is given by the Jacobi-cs function [32]: To find the parameters p and δ, one should use the normalization condition and the boundary condition due to the delta-function potential: Note that there are other Jacobi-elliptic functions that can solve the GPE, for example, the Jacobi-sc function: This function, however, does not lead to a physical solution in the limit L → ∞, and therefore we do not consider it here. We refrain from discussing any further solutions, which may, for example, correspond to the scattering solution, Eq. (5), from the main text. It turns out that hard-wall boundary conditions make it harder to find correct solutions for systems with finite L.

Appendix D. Zero-Density Limit within Mean-field Approximation
In this Appendix, we provide some technical details for the results presented in Sec. 4.

Appendix D.1. Many-body bound state
We first of all notice that the solution from Eq. (4) presented in the main text for periodic boundary conditions and Eq. (C.1) from the previous appendix for closed boundary conditions are identical in the limit of L → ∞, p → 1: Note that by definition x > 1, therefore, we derive the condition for the existence of the solution: which is in agreement with the PoT condition Eq. (10) from the main text.
For the chemical potential, we derive: The solution with the critical particle number still supporting a many-body bound state reads f (x) = π 2 M gL 2 δ 2 (N − 1) 1 cos π(x−L/2) δL .

Appendix E. Few-body Limit of the Artificial Atom
In this appendix, we discuss the smallest non-trivial system with N = 2 assuming that only two or three bosons can be bound to the impurity. Note that a priori it is not clear that the mean-field solution is applicable to such a few-body system. For a two-particle system, the IM-SRG becomes essentially exact (there is only an error due to the finite Hilbert space, see Appendix A, which can be easily controlled). This allows us to benchmark our mean-field results for large ring sizes. The results are presented in Fig. E1 where panels a) and b) show the density of the Bose gas for the cases where N cr = 2 and N cr = 3. It can be readily seen that independently of the ring size the behavior of the density is in an excellent agreement between the mean-field and IM-SRG approach. The same holds for the corresponding phase fluctuations depicted in panels c) and d) for distinct sizes of the ring. Note that the density for the largest ring size, i.e. LM |c| = 2.5N , in both panels is too low to render meaningful values of phase fluctuations (cf. Eq. (20)). However, even for the largest ring size, phase fluctuations are still low and the Bose gas can be adequately approximated with the mean-field ansatz. This numerical observation shows that the physical picture given in Sec. 2.2 is accurate even for the smallest set-ups.  8)). Panels c) and d) demonstrate phase fluctuations whose non-zero values reveal the presence of beyond-mean-field correlations. The dashed curves are provided to guide the eye. The data show results for two bosons (N = 2) for different ring sizes. Panels a) and c) are for systems with α = −0.5 (N cr = 2), and panels b) and d) for α = −1 (N cr = 3). The numerical error bars are calculated according to the prescription given in Appendix A.