Impurities in a one-dimensional Bose gas: the flow equation approach

A few years ago, flow equations were introduced as a technique for calculating the ground-state energies of cold Bose gases with and without impurities. In this paper, we extend this approach to compute observables other than the energy. As an example, we calculate the densities, and phase fluctuations of one-dimensional Bose gases with one and two impurities. For a single mobile impurity, we use flow equations to validate the mean-field results obtained upon the Lee-Low-Pines transformation. We show that the mean-field approximation is accurate for all values of the boson-impurity interaction strength as long as the phase coherence length is much larger than the healing length of the condensate. For two static impurities, we calculate impurity-impurity interactions induced by the Bose gas. We find that leading order perturbation theory fails when boson-impurity interactions are stronger than boson-boson interactions. The mean-field approximation reproduces the flow equation results for all values of the boson-impurity interaction strength as long as boson-boson interactions are weak.

For a single mobile impurity, we use flow equations to validate the mean-field results obtained upon the Lee-Low-Pines transformation. We show that the meanfield approximation is accurate for all values of the boson-impurity interaction strength as long as the phase coherence length is much larger than the healing length of the condensate.
For two static impurities, we calculate impurity-impurity interactions induced by the Bose gas. We find that leading order perturbation theory fails when bosonimpurity interactions are stronger than boson-boson interactions. The mean-field approximation reproduces the flow equation results for all values of the bosonimpurity interaction strength as long as boson-boson interactions are weak.

Introduction
The flow equation approach is an ab-initio method for solving many-body problems [3]. A related method, the in-medium similarity renormalization group (IM-SRG), was recently developed and successfully used in nuclear physics (see, e.g., Refs. [4,5]) 1 . The main concept behind both methods is the same. Therefore, we use 'IM-SRG' and 'flow equations' interchangeably in this paper.
IM-SRG was recently extended to cold Bose gases [1,2]. It was tested by calculating the ground-state energies of the Lieb-Liniger model and a one-dimensional (1D) Bose gas with an impurity atom ('Bose polaron') [1,2]. Those works demonstrate that flow equations allow one to go beyond mean-field approximation without relying on many-body perturbation theory.
It has been noticed that the mean-field approximation (MFA) in a frame co-moving with the impurity can accurately describe the self-energy of the impurity in a weakly-interacting Bose gas [2,17,19,22]. This observation is somewhat counter-intuitive, since strong phase fluctuations in one spatial dimension require a beyond-mean-field treatment. A counterargument to this point can be based on the observation that (when solving a Bose-polaron problem) one is usually interested only in what happens to a Bose gas in the vicinity of the impurity, and therefore, the absence of long-range order is not necessarily relevant. Therefore, the mean-field approach can be useful (cf. [26]) as long as the phase coherence length, ξe √ π 2 /γ , is larger than the length scale associated with the polaron, which is of the order of ξ [2]. Here ξ is the healing length of the condensate, and γ is the dimensionless Lieb-Liniger parameter which characterizes the boson-boson interaction strength, see Sec. 3. This argument implies that as long as e √ π 2 /γ 1 one can use the MFA to study the Bose-polaron problem. The discussion above is based on perturbation theory, and further work is needed for its rigorous proof. Here, we use the IM-SRG method to provide numerical evidence for its validity.
In this work, we benchmark mean-field results against those obtained with IM-SRG, and find that the MFA can describe the density of a Bose gas with an impurity particle accurately. To confirm the absence of boson-boson entanglement in the frame co-moving with the impurity, we calculate the phase fluctuations. They turn out to be negligible for the considered systems. Finally, we use the Born-Oppenheimer approximation to estimate the potential between impurities supported by a Bose gas. Our IM-SRG results show that the mean-field approximation is useful to study mesoscopic and large systems with two impurities. In particular, the MFA allows one to calculate the induced impurity-impurity interaction potential beyond first-order perturbation theory in a simple manner, i.e., without including any information about beyond-mean-field boson-boson correlations.
The paper is organized as follows: In Sec. 2, we give a short review of the IM-SRG, and explain how to calculate observables using this method. In Sec. 3, we discuss a Bose gas with a single impurity assuming repulsive boson-impurity interactions. There, we benchmark the IM-SRG results for the density against the exact Bethe-ansatz solution. Furthermore, we calculate the density and phase fluctuations of the Bose gas for systems that do not allow for an exact analytic treatment. In Sec. 4, we consider a Bose gas with two impurities, and calculate induced impurity-impurity interactions. We show that two impurities attract each other for repulsive impurity-boson interactions, whereas two impurities repel each other if one impurity attracts and another repels bosons, in accord with Refs. [27,23]. We summarize our results and give an outlook in Sec. 5. Further details on the IM-SRG method in our implementation are given in Appendix A. For convenience of the reader, we present some additional results for a system with attractive boson-impurity interactions in Appendix B.

In-Medium Similarity Renormalization Group
For convenience of the reader, we shall present in this chapter the main ingredients of the IM-SRG method for bosons, see also Ref. [1], Appendix A, and Fig. 1.

Flow equations
The IM-SRG is an extension of the SRG [28,29,3] based upon the flow equation which transforms the Hamiltonian matrix into a block-diagonal form, i.e., it decouples the "ground-state" matrix element from all excitations (see Fig. 1). The flow equation is defined once the initial condition, H(s = 0), and the generator of the transformation, η, are specified. It is worth noting that Eq. (1) is equivalent to the unitary transformation H(s) = U (s)H(s = 0)U † (s), assuming that the antihermitian operator η and the unitary operator U are connected as We prefer to write the unitary transformation in the form of Eq. (1) because it allows us to choose the operator η(s) during the flow, i.e., for every parameter s, and, hence, to steer the flow in the desired direction. In our work, η(s) is chosen from the matrix elements that describe the couplings between the 'condensate' and its excitations such that these couplings become weaker as the flow progresses, see Fig. 1. A detailed construction of η(s) is presented in Appendix A.
In practice, the following steps constitute the IM-SRG method: (i) find a one-body basis to write a Hamiltonian matrix in second quantization (see Appendix A.2), (ii) find a reference state to normal order the operators (see the next subsection), (iii) solve the flow equation (1) (we use the explicit Runge-Kutta method 5(4)).

Normal ordering
In general, it is impossible to solve Eq. (1) for a many-particle system without approximations. The complexity is due to the commutator [η, H]: It leads to many-body terms, which are not present in the initial Hamiltonian H(s = 0). To solve Eq. (1), the many-body terms must be truncated at some order. To define a truncation hierarchy, we write the Hamiltonian H in second quantization using normal ordering with respect to a reference function, Ψ ref , which is a product state, see the next subsection. Upon normal ordering, we truncate threebody excitations and beyond, see Fig. 1. To estimate the introduced truncation error, we use the three-body elements and second order perturbation theory for matrices, see Fig. 1 and Appendix A.1.  (1) on the Hamiltonian. The Hamiltonian matrix is unitarily transformed to a block-diagonal form, such that the ground state (hatched green) becomes decoupled. The upper row illustrates the exact transformation without a truncation. Note that many-body excitations appear during the flow. The bottom row illustrates a truncation scheme adopted to circumvent this problem. The induced three-body terms (checkered blue) can be estimated and used to evaluate the accuracy of the truncation scheme. The excitations are defined with respect to the adopted normal ordering, see Appendix A.
To normal order a Hamiltonian with one-and two-body operators, we define the contractions 2 , following Ref. [1] : : where a † α is a bosonic creation operator, The parameter N is the number of bosons. I is the identity operator; the operator P α 1 α 2 swaps the indices α 1 and α 2 .
A generic Hamiltonian with one-and two-body operators, reads in the normal ordered prescription as with N = E is the energy of the ground state, f α 1 α 2 (Γ α 1 α 2 α 3 α 4 ) describes one-particle (twoparticle) excitations.

Reference state
The reference state, Ψ ref , should approximate an eigenstate (here the ground state) of the Hamiltonian well, otherwise the IM-SRG transformation cannot map Ψ ref onto the exact state. Since we are interested in ground-state properties of a bosonic system, it is logical to use a product state as a reference state, i.e., where x α is the coordinate of the αth boson, and f is some function whose form we discuss below. The choice of a product-state ansatz is natural for cold Bose gases with macroscopic population of a single mode, i.e., with a large condensate fraction. However, for an interacting one-dimensional Bose gas, the reference state of Eq. (10) is in general not accurate. In the thermodynamic limit, correlations in the 1D Bose gas decay algebraically, precluding condensation [30,31,32,33,34]. However, as our analysis below shows, the product state is a useful starting point for analyzing properties of Bose polarons.
In this work, we construct Ψ ref using either f 1b or f GP . The function f 1b is the ground-state wave function of the one-boson Hamiltonian as in Refs. [1,2]. The second function is obtained within a mean-field approximation, i.e., f GP is the solution of the Gross-Pitaevskii equation.
f 1b and f GP are real functions in our work. This choice does not affect the generality of our results, since the ground state of our problem can be described using a real wave function.
To distinguish the IM-SRG method with f 1b from IM-SRG with f GP , we introduce the notation IM-SRG(f 1b ) and IM-SRG(f GP ), respectively 3 . It is worth noting that one can rely on an iterative procedure to find a good reference state, starting from any reasonable initial guess f a , which signals that the results are converged. In addition, this procedure can be used to validate the convergence of our results. We have checked that the results within the zeroth-order iteration (i.e., of IM-SRG(f

Observables
In nuclear physics, the IM-SRG method was used not only to calculate the energy, but also to estimate other observables [35,5]. One of the goals of the present paper is to develop (and test the accuracy of) the IM-SRG method for calculating the density and phase fluctuations of cold Bose gases.
In order to calculate observables other than the energy, the corresponding Hermitian operator O should be transformed together with the Hamiltonian. To this end, we write O in second quantization, normal order it with respect to the reference state Ψ ref , and solve the flow Equations (1) and (11) are solved simultaneously since the generator η depends on H.
In this work, we focus on calculating one-body observables. The commutator in Eq. (11) leads to two-and higher-order terms for such observables at s > 0, which should be truncated according to our scheme. We cannot estimate the associated error using the strategy adopted for the energy (see Appendix A.1), as in general the operators O and H do not commute. Instead we define the "relative truncation error" as where e is the energy calculated using flow equations and δe is our estimation of the truncation error for the energy, (see Appendix A.1). We estimate the error due to truncation for O as By comparing to the exact density, we will show below that δO can estimate accurately the error of the IM-SRG. However, we do not expect this always to be the case. In general, one cannot infer the accuracy of an observable from the accuracy of the energy using a linear approximation 4 , which means that δO is no more than a useful phenomenological estimate.

A Bose Gas with a Single Impurity Atom
To illustrate calculations of observables based upon flow equations, we investigate a onedimensional system of N bosons and a single impurity atom. This system, which is often referred to as the 1D Bose-polaron problem, is one of the simplest models where our approach is useful. The corresponding Hamiltonian in first quantization reads as where y is the position of the impurity, x i is the position of the ith boson, m is the mass of the impurity atom, and M is the mass of a boson. To model atom-atom interactions, we use zero-range potentials [36]: where c defines the strength of the boson-impurity interactions, and g determines the bosonboson interactions. We consider periodic boundary conditions, i.e., particles are confined to a ring of length L, see Fig. 2. The average density of the Bose gas is ρ = N/L. We introduce the dimensionless set of parameters: where E is the energy of the system. The parameter γ is also known as the Lieb-Liniger parameter. In this new set of units the HamiltonianH reads as: 4 To illustrate this statement, let us assume that a numerical method produces the following approximation to the ground state: ψ0 + αf , where ψ0 is the exact ground-state wave function, and f is an element of the Hilbert space, which is orthogonal to ψ0, i.e., ψ0|f = 0. If the numerical method is accurate, then α → 0. In this case, it is easy to see that the error produced by the numerical method is proportional to α 2 for the expectation value of the energy. However, the error for a general observable can be much larger, as it scales as α. For convenience, the problem is solved using the set of coordinates {z i }, which describe the relative distances between the impurity and the bosons.
In the following we will omit the tilde for better clarity. We focus on weakly-interacting Bose gases (i.e., γ 1) that are 'large enough' (i.e., the healing length, ξ = 1/( √ γρ), is smaller than L). Furthermore, we assume that the phase coherence length is larger or comparable to L, so that the Bose gas is in a (quasi)-condensed state and our reference state is accurate.
In this work, we mainly focus on repulsive boson-impurity interactions (c > 0). In Appendix B we present preliminary results for the Bose-gas density and phase fluctuations for attractive interactions. However, further studies are needed for c < 0 and will be addressed in a future publication [37]. The self-energy of the impurity for c > 0 was already calculated with IM-SRG in Ref. [2]. We are now going to compute other observables using this numerical method.
To that end, we transform the Hamiltonian (14) to the system of coordinates co-moving with the impurity, i.e., we introduce a new set of coordinates where θ(x) is the Heaviside step function. The coordinates z i allow one to calculate meanfield properties of the Bose-polaron problem in a simple manner, see Refs. [2,17]. The transformation {y, x i } → {z i } is related to the unitary Lee-Low-Pines transformation [38] performed in coordinate space [22,19]. In the new coordinates, the Hamiltonian (14) is written as where P is a quantum number -the total (angular) momentum of the system. We consider the case P = 0, as it corresponds to the ground-state manifold. The Hamiltonian H P describes a system of N bosons, and can be easily written in the language of second quantization using the annihilation and creation operators, a i and a † i , which are defined in the frame co-moving with the impurity.

Reference state
In this section, we introduce the reference states, which are needed to solve the Bose-polaron problem with flow equations. The first reference state is built upon the ground state, f 1b , of H P =0 for a single boson, i.e., N = 1. The corresponding Schrödinger equation reads as where κ = m/(1 + m) is the reduced mass, and z ∈ [0, N ] is the coordinate of a boson. The ground-state solution that satisfies the boundary conditions The second reference state, f GP , solves the Gross-Pitaevskii equation that corresponds to H 0 : where µ is the chemical potential. The solution to this equation is given by where sn is the sn−Jacobi elliptic function, and K(p) is the complete elliptic integral of the first kind [39]. The parameters p ∈ [0, 1) and δ are fixed by the boundary conditions due to the delta-function potential , and by the normalization condition [ |f GP (z)| 2 dz = 1]. The corresponding chemical potential µ reads as: The mean-field solution f GP is discussed in more detail in Ref. [2], see also Refs. [19,22] for the discussion of the thermodynamic limit, and Refs. [40,41] for the discussion of the limit c → ∞.

Density
In this section, we calculate the density, ρ(z) = Φ gr |ρ(z)|Φ gr , of the Bose gas in the frame co-moving with the impurity: where Φ gr is the ground state of H 0 . ρ(z) should not be confused with the density of the Bose gas without the impurity, ρ = N/L. To use the density operator ρ(z) in the flow equation approach (see Sec. 2.4), we write it in second quantization as where φ i (z) is the ith element of the one-body basis employed for writing the Hamiltonian in second quantization. Note that in our implementation the basis {φ i (z)} depends on the used reference state, see Appendix A.2.
To test flow equations, we first calculate ρ(z) assuming equal masses (m = M ) and equal interactions (c = γ). For the ground state, these assumptions turn our system into the exactly solvable Lieb-Liniger model for N + 1 particles [43] in a ring, since we can no longer distinguish between the impurity and a boson. To calculate the density of the bosons in the frame co-moving with the impurity (see Eq. (23)), we relate it to the pair correlation function of the Lieb-Liniger model, g LL : To derive this relation, we have used the following definition of g LL : where Ψ LL is the ground-state wave function of the Lieb-Liniger model. For z > 0, we can write that which in combination with Eq. (23) leads to Eq. (25). For certain parameter regimes, the function g (2) LL is known for few-body systems (see, e.g., Ref. [42]), and we use those results to benchmark our findings, see Fig. 3 for N = 6. The density calculated using flow equations agrees well with the exact values for all values of z. Near the impurity, the density of the bosons is suppressed, since the boson-impurity interaction is repulsive. The presented error bars show the error due to the truncation of the Hilbert space (see, Appendix A.3), and due to the truncation of many-body forces in flow equations, see Eq. (13). For the considered parameters, the latter dominates. All in all, the comparison to the exact results allows us to conclude that our error estimate is accurate. Figure 3 shows that IM-SRG(f 1b ) and IM-SRG(f GP ) agree, which means that both f 1b and f GP are suitable reference states for the considered parameters. In our studies, we noticed that the reference state f GP is generally a better choice than f 1b . In comparison to IM-SRG(f 1b ), the scheme IM-SRG(f GP ) allows us to obtain converged results for a larger range of parameters. In particular, IM-SRG(f GP ) is more reliable for large systems, and large boson-boson interactions. Figure 3 explains this observation: The more complicated meanfield function f GP provides a better approximation of the exact density, and, hence, it is easier for the flow equation method to map this reference state onto the real ground state of the system. In what follows, we present our results only for IM-SRG(f GP ).
Finally, we calculate the density for parameters for which the system is no longer exactly solvable. Our goal here is to test the mean-field treatment of H 0 . It is already known that this treatment can produce accurate results for the energy of the impurity and the contact parameter [2,19,22]. Here, we show that it can also be used to calculate the density, see Fig. 4. As in Fig. 3, the density of the Bose gas is suppressed near the impurity. The length scale that characterizes this suppression is given by ξ/ √ κ [2,22] (for γ = 0.02, ξρ 7.1 and for γ = 0.1, ξρ 3.2). All in all, we observe that the mean-field approximation describes ρ(z) accurately for all values of the boson-impurity interaction strength c, as long as the parameter γ is small 5 .
We have checked that the mean-field approximation is accurate for up to γ 0.5 and c → ∞ by comparing to the Monte-Carlo results presented in Ref. [14]. The comparison of our IM-SRG results to the RG results 6 of Ref. [14] suggests that it is more advantageous to work with a real-space formulation of the Bose-polaron problem. Large values of c require beyond-Fröhlich-polaron treatment of the problem in momentum space, in particular, one should include phonon-phonon interactions. In contrast, in our implementation, already mean-field results are accurate. The accuracy of MFA is probably not surprising after we notice that the (phase) coherence length is larger than the length scale we are interested in. For instance, the phase coherence length for γ = 0.1 is about 20000ξ. We illustrate this statement further

Phase fluctuations
Phase fluctuations are strong in the one-dimensional world [33,34,46,47], incapacitating the mean-field treatment. However, as long as one is interested in the physics on the length scales smaller than the coherence length, the mean-field approach can give accurate results. We will now calculate phase fluctuations using flow equations, and explicitly justify the use of the MFA 7 . Another way to validate the mean-field ansatz could be based on calculating the condensate fraction in the considered mesoscopic ensembles. For example, the flow equation approach predicts that the condensate fraction 8 for systems with N = 50 and γ ≤ 0.1 is always large 7 Phase fluctuations vanish in the MFA. 8 By condensate fraction, we mean here the expectation value of the operator a † 0 a0/N .
(∼ 95%), thus endorsing the use of the mean-field ansatz for these systems. In particular, this large condensate fraction suggests that if fifty bosons can screen the impurity, then the MFA should yield accurate results for the properties of the impurity in the thermodynamic limit. However, the condensate fraction can only be used to provide a phenomenological argument in support of the mean-field approach. In one spatial dimension, the condensate fraction depends on the considered number of particles, and vanishes in the thermodynamic limit. Therefore, we do not discuss it further.
To calculate phase fluctuations, we first compute the one-body density matrix and then extract the phase fluctuations δΦ 0z using the expression [46]: The result is shown in Fig. 4. For γ = 0.02, phase fluctuations are negligibly small 9 . For γ = 0.1, phase fluctuations play a more important role, however even then they can be neglected so that ρ(0, z) ρ(0)ρ(z). This implies that for these parameters the Bose gas can be described using a mean-field ansatz. One could anticipate that phase fluctuations depend noticeably on the value c, which determines the density of bosons and, hence, the effective Lieb-Liniger parameter in the vicinity of the impurity. However, we observe that phase fluctuations depend only weakly on the boson-impurity interaction strength, c. Note that in practice it is difficult to extract phase fluctuations according to the definition (29) in systems with c/γ 1 for which ρ(0) → 0. Therefore we do not compute phase fluctuations in this limit.

Contact parameter
This section focuses on repulsive boson-impurity interactions, and leaves a thorough IM-SRG study of the attractive case for the future [37]. However, we must note that our conclusion that the mean-field approximation describes accurately an impurity in a Bose gas should no be straightforwardly extended to a Bose gas with an attractive impurity. There is an important difference between the cases with c > 0 and c < 0: The latter supports a formation of tightly bound states at large negative values of c. Many-body bound states are highly-correlated and cannot be accurately described by the mean-field ansatz, unlike the opposite case with c → ∞. This implies that the region of applicability for c > 0 must be larger (see also the discussion on the effective mass in Ref. [22]). To elaborate slightly more on this difference, we calculate the density of the Bose gas at the position of the impurity, i.e., the contact parameter, C = ρ(0)/ρ, in the mean-field approximation in the thermodynamic limit: where the positive sign of the exponent is for c > 0 (see Ref. [2]) and the negative sign should be taken for c < 0. We compare the parameter C M F to the contact parameter calculated using Quantum Monte-Carlo [12] in Fig. 5. For c > 0, the agreement between Monte-Carlo and the mean-field approximation is reasonable for all available data points. For attractive interactions, the difference between the results is more noticeable, which implies that the MFA leads to less accurate results for c < 0, see also Appendix B, where we present some additional data for the case with attractive interactions. Note in particular Fig. 15, which indicates large phase fluctuations for moderate impurity-boson interactions, in contrast to the repulsive case.

A Bose Gas with Two Impurity Atoms
In the previous section, we considered a single impurity in a Bose gas, which is the standard starting point in the analysis of systems with impurities. However, the physics of systems with many impurity atoms can be drastically different from that of a system with a single impurity atom, in particular, because the Bose gas mediates interactions between impurities. Therefore, the next step for a reliable description of a (quasispin-)polarized systems must be an assessment of the strength of the induced impurity-impurity interaction potential. One possible way to do this is to consider a Bose gas with two mobile impurities, see, e.g., Refs. [25,20,48,49]. We choose another approach. We estimate the induced impurity-impurity interaction using the Born-Oppenheimer approximation [50,51,52] (see Refs. [53,54,55,56,57] for related studies in three spatial dimensions), i.e., for m → ∞. It is known that the Born-Oppenheimer approximation captures short-range correlations, which define overall properties of impurityimpurity interactions [50,58]. Long-range impurity-impurity correlations induced by quantum fluctuations are beyond the range of applicability of the Born-Oppenheimer approximation, in particular, since they can depend on the mass of the impurity [27,59,21]. These long-range Figure 6: Illustration of a system of N bosons and two static impurities confined to a ring of length L. The ith boson has the coordinate x i . One impenetrable impurity (c 1 → ∞) is placed at y 1 = 0, and another one with the interaction strength c 2 is at y 2 = d.
correlations are not relevant for our discussion -they are weak for the considered parameters regimes, and can be neglected. In particular, these correlations are relevant only at distances 5ξ (see, e.g., [51]), i.e., they can be extracted only by considering systems that are larger than studied here.
The dimensionless Hamiltonian for two static impurities is written in first quantization as: where c 1 and c 2 describe the strength of the impurity-boson interactions, and y 1 and y 2 are the positions of the impurities. Without loss of generality, we place one impurity at y 1 = 0 and the other at y 2 = d, see Fig. 6. We assume that the impurity at y 1 is impenetrable, i.e., 1/c 1 = 0. In other words we consider an impurity in a semi-infinite Bose gas [23,24]. This assumption allows us to simplify the presentation. We will show that for c 2 > 0 the impurities attract each other, whereas if c 2 < 0 the impurities repel each other. In the former case, the energy is minimized when two impurities are on top of each other, whereas in the latter case the attractive impurity wants to be far from the repulsive one, see also Ref. [23]. We shall consider these cases separately.
The goal of this section is to calculate induced impurity-impurity interactions in the Born-Oppenheimer approximation using the flow equation approach. To this end, below, we compute with IM-SRG the ground-state energy of the Hamiltonian H 2 . We also show that the induced interactions can be accurately calculated using the mean-field approximation, at least for weak boson-boson interactions. Finally, we estimate when first-order perturbation theory, which is commonly used to estimate impurity-impurity interactions [24,51], fails.

Reference state
To use the IM-SRG(f GP ) scheme, we first find the reference state, f GP , by solving the Gross-Pitaevskii equation: The solution to this equation reads where the parameters p 1 , p 2 , δ 1 , δ 2 , a are determined by the conditions: The chemical potential is µ = 2 p 1 +1 δ 2 1 N 2 K(p 1 ). The mean-field solution presented here is valid for any positive values of c 1 . For the special case 1/c 1 = 0 that we consider in this section, one should set the parameter a to 0, and solve only the boundary conditions (34)- (37). Once the function f GP is obtained, the mean-field energy of the system is calculated as

Results
We first calculate the induced interaction potential, 2 = E 2 (c 2 , d)−E 2 (c 2 = 0, d), where E 2 is the ground-state energy of H 2 . Our results for this quantity for N = 60 particles are presented in Fig. 7. Note that E 2 (c 2 = 0, d) = E 2 (c 2 , d = 0), hence 2 (d = 0) = 0. The potential 2 is attractive, because the boson-impurity repulsion is minimized when the two impurities are on top of each other. The considered system is not yet in the thermodynamic limit (see the next subsection), but it is large enough to see the saturation of the impurity-impurity interaction at large values of d, i.e., far from the impenetrable impurity. In Fig. 7, we compare the results obtained via IM-SRG(f GP ) and perturbation theory. The latter assumes that the second impurity does not affect the density of the Bose gas, and therefore 2 = c 2 n(d), where n(d) is the density of the Bose gas for c 2 = 0. This is the common assumption for estimating the induced interaction [23,51,24]. The perturbation theory leads to the Yukawa-type potential when both impurities are weakly interacting [50,58,25,51]. Hence, our results indicate the limits of applicability of that standard potential. The figure also presents the mean-field approximation, which uses Eq. (39) to estimate 2 . We conclude that perturbation theory can be used if c 2 < γ. However, it fails already for c 2 γ, and more involved calculations are required to find the induced potential in this regime. Perturbation theory implies much stronger impurity-impurity interactions. Its use will lead to wrong predictions for a number of experimentally relevant observables, such as the limits of stability of the polaronic gas [60]. For all considered parameters, the mean-field approximation agrees with the IM-SRG(f GP ) method. This implies that the MFA can be used to calculate the induced potential beyond first order perturbation theory. Finally, we note that far from the impenetrable impurity, the energy 2 does not approach the self-energy of a single impurity, E 2 (c 1 = 0, c 2 ) − E 2 (c 1 = 0, c 2 = 0), (dotted lines in Fig. 7) 10 . This is a finite-size effect. Far from the impenetrable impurity, the density of the Bose gas in a finite system is larger than ρ, see Fig. 4, which leads to the difference at d = 20 between the dotted lines and dots in Fig. 7.
We also employ the IM-SRG(f GP ) to calculate the density and phase fluctuations of the Bose gas for N = 60, γ = 0.02, c 2 = 0.1 and d = 11, 20. The results are presented in Fig. 8. The density vanishes at x = 0 and x = 60 due to the impenetrable impurity and periodic boundary conditions. The Bose gas is strongly affected by the second impurity, which is located at x = d. This explains the discrepancy between perturbation theory and the flow equation approach in Fig. 7. All in all, the IM-SRG(f GP ) results for ρ(x) agree well with the mean-field prediction.

Figures 8 c) and d) show phase fluctuations.
Here, we choose the position of the second impurity, d, as the reference point. As for a single impurity in a ring, phase fluctuations increase far from d 11 . The maximum value of δΦ dx is small (in agreement with the results of the previous section), and does not depend strongly on d. A condensate fraction for the considered system is about ∼ 95%. Our findings presented in this subsection imply that the physics of two strongly interacting impurities in a weakly-interacting Bose gas can be conveniently studied using the mean-field approximation.

Approaching the thermodynamic limit
In this subsection, we study the transition to the thermodynamic limit, i.e., we increase N and L, while keeping the density constant, ρ = N/L. The key question here is how many bosons are needed to simulate the infinite system. For a single impurity in a weakly-interacting Bose gas, this was briefly considered in Ref. [2]. Here, we discuss what happens for two impurities. The IM-SRG results for c 2 = 0.1 and c 2 = 0.02, γ = 0.02 and different numbers of particles are shown in Fig. 9 a) and b). For the considered values of N , the induced interaction is far from the thermodynamic limit, i.e., it changes with the numbers of particles. The thermodynamic limit is reached for system sizes which are beyond the flow equation approach. Still, we can use IM-SRG to predict the induced potential in the thermodynamic limit by fitting the IM-SRG energies to the function C 1 + C 2 N C 3 , where C 1 , C 2 and C 3 are fitting parameters. The parameter C 1 defines the potential in the thermodynamic limit. The fitting parameters C 2 and C 3 have no direct physical interpretation 12 . In Figs. 9 a) and b), we present the value of C 1 . The truncation error in the IM-SRG method grows rapidly with the number of particles. This rapid growth rules out a reliable extrapolation of the errorbars to the thermodynamic limit. Therefore, we give no estimate for the accuracy of C 1 , which leads to an apparently oscillating character of the potential in the thermodynamic limit. We expect that the exact potential is a monotonically increasing function of the distance between the impurities, d.
In Figs. 9 c) and d), we compare our estimate for the potential in the thermodynamic limit with the result based upon perturbation theory, which includes quantum fluctuations [23,24]. For weak boson-impurity interactions, c 2 = γ, both curves agree well. For larger interactions, however, the curves deviate. In this case, the density of bosons is strongly influenced by the second impurity (see Fig. 8), and, therefore, perturbation theory is not a valid approximation. We draw two conclusions from Fig. 9. First, high compressibility of a weakly-interacting Bose gas leads to a large number of particles needed to reach the thermodynamic limit. This should be contrasted with systems of strongly interacting bosons or fermions, for which a handful of particles can screen the impurity [61], for more detail, see Refs. [62,63] and references therein. Here we find that for γ = 0.02, one needs more than 100 particles to capture the effective short-range interaction between two static impurities in the thermodynamic limit. This implies that any study that aims to relate measurements in current cold-atom set-ups with small number of particles to the thermodynamic limit must provide an estimation of finite-size effects. This is especially important for any prospective experimental study of induced long-range interactions. Second, calculations beyond first-order perturbation theory are required to estimate induced impurity-impurity interactions also in the thermodynamic limit. As we show here, these calculations can be based upon the mean-field approximation.  We also show our prediction for the induced interaction potential in the thermodynamic limit. No errorbars are shown for this case, since we cannot estimate them reliably, see text. In the lower panels, we compare our prediction for the thermodynamic limit (orange squares) with the corresponding result based upon perturbation theory, which includes quantum fluctuations [23] (solid blue curve), and the mean-field result (green dots). Panels c) and d) are for c 2 = 0.02 and c 2 = 0.1, respectively. The dotted lines in c) and d) show the self-energies of a single impurity, E 2 (c 1 = 0, c 2 ) − E 2 (c 1 = 0, c 2 = 0) in the thermodynamic limit while the dashed curves are added to guide the eye.

Attractive case, c 2 < 0
In this section, we consider attractive boson-impurity interactions, c 2 < 0. This case is more complicated because all bosons are bound to the impurity if γN 2|c 2 | in the limit L → ∞ (assuming that N is fixed), see [64,37]. Therefore, to obtain a meaningful estimation for the induced interactions, we must consider N 2|c 2 |/γ.

Reference state
For attractive interactions, we do not obtain the reference state f GP from the Gross-Pitaevskii equation. Instead, we choose the reference state as where f one−rep is the mean-field solution for a Bose gas with one repulsive impurity, and f one−attr is the mean-field solution for a Bose gas with a single attractive impurity. Therefore, the function f GP in Eq. (40) is the full mean-field solution for two impurities in the limit of large separation between the impurities, i.e., d 1.
Otherwise, the function f GP is an approximation to the solution of the Gross-Pitaevskii equation. We observe that f GP is an accurate approximation, see, in particular, Fig. 11, where f GP is plotted together with the density of the bosons calculated using flow equations.
The function f one−rep is given by (see Sec. 3.1) where δ 1 = 1, since we work with an impenetrable impurity, c 1 → ∞. The function f one−attr reads as [37] f one−attr (x) = 4K(p 2 ) 2 γδ 2 2 N 2 (N − 1) The parameter N is given by the normalization condition:

Results
We compare the induced interaction potential obtained with the IM-SRG(f GP ) to the one obtained using first-order perturbation theory, see Fig. 10. In contrast to the case of c 2 > 0, now the induced impurity-impurity potential is repulsive. The attractive impurity maximizes its energy by being far from the hole in the density of bosons created by the repulsive impurity.
In agreement with the case c 2 > 0, perturbation theory fails to describe impurity-impurity interactions when |c 2 | γ. The important difference is that now perturbation theory leads to a weaker induced interaction in comparison to the IM-SRG result. As in the previous case, perturbation theory fails because the density of the Bose gas is strongly modified in the vicinity of the impurity for |c 2 | γ. To illustrate a strong modification of the density, we calculate ρ(x), and phase fluctuations of the Bose gas for N = 60, γ = 0.02, c 2 = −0.1 and d = 11, 20 within the IM-SRG(f GP ), see Fig. 11. Figure 11 demonstrates that the reference state f GP gives an accurate approximation to the exact density of the Bose gas. Phase fluctuations are small, and we have checked that the condensate fraction for these parameters is above 95%. Therefore, our conclusion for the considered parameters is similar to the case with c 2 > 0: The mean-field approach can be used to describe the impurity-impurity mediated interactions as long as the boson-boson interaction is weak. Although, we do not demonstrate this here, we expect that f GP and the mean-field approach in general to be less accurate for one attractive impurity in comparison to one repulsive impurity when |c 2 | γ. This expectation is based on our analysis of the contact parameter, see

Approaching the thermodynamic limit
Finally, we calculate the induced interactions for different values of N , in order to understand the few-to many-body transition in this system. We illustrate our findings in Fig. 12 for γ = 0.02 and c 2 = −0.02, c 2 = −0.1. Like in the case with c 2 > 0, we fit the energy to estimate the induced impurity-impurity potential in the thermodynamic limit, and compare the estimate to the result of Ref. [23]. Just like in the case with c 2 > 0, we see that perturbation theory is accurate for weak boson-impurity interactions but fails to describe stronger interactions.

Summary and Outlook
The paper explores the possibility to calculate observables for degenerate Bose gases using the flow equation approach. For illustration purposes, the focus is on a one-dimensional Bose gas with one and two impurity atoms. The considered system allows us to benchmark the IM-SRG results against the existing exact data based upon the Bethe ansatz, and to study in detail the Bose-polaron problem and polaron-polaron interactions, topics of current theoretical interest. Pert. + QF., Ref. [23] Thermod. Limit c 2 > 0, Thermod. Limit Figure 12: Interaction potential between two impurities induced by the Bose gas for γ = 0.02 in the attractive case, c 2 < 0, with two different strengths a) c 2 = −0.02 and b) c 2 = −0.1. The curves show different values of N listed in the insets. We also give our prediction for the induced interaction potential in the thermodynamic limit. No errorbars are shown for this case, since we cannot estimate them reliably, see text. In the lower panels, we display our prediction for the values in the thermodynamic limit (orange dots) and the corresponding result based upon perturbation theory, which includes quantum fluctuations [23] (solid blue curve). Panels c) and d) are for c 2 = −0.02 and c 2 = −0.1, respectively. For comparison, we also show the negative IM-SRG result for c 2 > 0 (green squares). The dashed curves are plotted to guide the eye while the dotted lines in c) and d) give the self-energy of a single impurity E 2 (c 1 = 0, c 2 ) − E 2 (c 1 = 0, c 2 = 0), in the thermodynamic limit.
In the single impurity case, we consider repulsive boson-impurity interactions (c > 0). It turns out that the density of the Bose gas for the repulsive Bose-polaron problem can be calculated accurately using the mean-field approximation in the coordinate frame, which is 'co-moving' with the impurity. To explain the validity of the mean-field approach, we show that the condensate fraction is large and phase fluctuations are small for the considered mesoscopic ensembles (N < 100), and weak boson-boson interactions. For attractive interactions (c < 0), the possibility of deep impurity-boson bound states complicates the analysis. This issue will be addressed in a future publication [37].
For two impurities, we calculate induced impurity-impurity interactions in the Born-Oppenheimer approximation. For simplicity, we assume that one impurity is impenetrable. The other one either attracts or repels bosons. The former case leads to repulsive, and the latter one to attractive induced interactions. We find that the mean-field approximation describes accurately these interactions, and can be used in straightforward calculations based upon the well-studied Gross-Pitaevskii equation. We also show that first order perturbation theory is valid when boson-impurity interaction is smaller or equal to the boson-boson interaction, however, fails otherwise. Finally, we discuss the few-to many-body transition, and show the importance of finite-size effects for impurities in weakly-interacting Bose gases.
Our results show that the mean-field approach is a robust tool to study weakly-interacting Bose gases with impurities. In the other limit of strongly interacting Bose gases, one can study the self-energy and the density of an impurity using tools developed for Fermi gases, e.g., based upon variational wave functions [65] or the Bethe ansatz [66]. In the future, it might be interesting to investigate the transition between these two limits, and the evolution of the Bose polaron into the impurity in a Tonks-Girardeau gas. A modification of the flow equation approach with two reference states might be useful to study this transition. One could compare the obtained IM-SRG results to those calculated using beyond-Gross-Pitaevskii effective theories (as, e.g., introduced in Ref. [67]).
The present work paves the way for IM-SRG studies of Bose gases in higher spatial dimensions. A starting point for such an extension might be a study of dilute bosonic droplets, for which a number of exciting analytical predictions exist [68,69]. It should also be possible to investigate two-and three-dimensional systems with impurities. The relevant mean-field solutions can be found in the literature [70,71,72,73], giving reference states for flow equations. A modification of the IM-SRG, which takes into account Hilbert space associated with an impurity can allow one to study composite impurities in Bose gases, and corresponding quasiparticles, in particular, angulons [74,75].
Note added after ArXiv submission: After submission of this manuscript, we learned of recent works [76,77], where the accuracy of the mean-field approximation to polaron-polaron interactions is also discussed.
for sharing with us the data for the densities presented in Ref. [14]. This work has received funding from the DFG

A Details on the Method
Below, we explain the IM-SRG method in more detail. We start by presenting the form of the flow equation (1) after our truncation of many-body terms, and discuss our estimate of the truncation error. After that, we explain further technical details of our implementation of the IM-SRG method. First, we discuss the one-body basis, which is used to represent the Hamiltonian in second quantization. Then, we discuss the truncation of the Hilbert space.
In the Appendix, the Einstein summation rule is implied, when Arabic indices are used. There is no summation for Greek indices.

A.1 Flow equation
Our flow equation reads as dH ds where the generator η is written as The matrices ξ ij and η ijkl must be chosen such that the couplings to the ground statethe matrix elements f i0 and Γ ij00 -vanish (see Fig. 1). Note that the commutator in the right-hand-side of Eq. (1) induces many-body terms: This three-body operator induces a three-body operator in the Hamiltonian at s > 0. Since all couplings from the ground state need to vanish, we would now need a three-body operator in the generator as well, which would in turn generate a four-body operator and so on. It is therefore impossible to treat the flow equation exactly and we need to truncate manybody terms. We choose to truncate at the two-body level, keeping only the terms from the three-body operator that contain at least one a † 0 a 0 operator (see Ref. [1] for a more detailed discussion). We neglect the remaining pieces (called W ), see also Fig. 1, which illustrates the used truncation scheme.
Upon truncation, we derive a closed system of coupled differential equations (note that Ref. [1] has typos in this system of equation, which we correct here): with D α 1 = 2 − δ α 1 0 , I α 1 α 2 = 1 + δ α 1 0 δ α 2 0 − 2δ α 2 0 , and To estimate the truncation error for the ground-state energy we use second-order perturbation theory [1] where Φ p is a state that contains three-body excitations.
The explicit choice of the generator can be justified a posteriori if the couplings to the excited states vanish [3]. Our choice of the generator is For this generator, the couplings (the matrix elements f i0 and Γ ij00 ) vanish, and the energy of the ground state, N , converges, see Fig. 13 for an illustration of a typical convergence pattern in our study. We expect that other standard choices of the generator, see, e.g., Ref. [5], lead to similar results, see also Ref. [1].

A.2 One-body basis
To write the Hamiltonian in second quantization, H = A ij a † i a j + 1 2 B ijkl a † i a † j a k a l , we need to calculate the matrix elements A ij and B ijkl using some one-body basis. In this section, we discuss the one-body basis used in the present study. First, we solve the single-boson problem whose eigenstates produce the basis set {φ i }. This set is used as a basis set when we work with the (single-body) reference state f 1b = φ 0 , where φ 0 is the ground-state of the single-boson problem. The corresponding contractions enjoy the simple form C α 1 α 2 = δ α 1 0 δ α 2 0 N and C α 1 α 2 α 3 α 4 = δ α 1 0 δ α 2 0 δ α 3 0 δ α 4 0 N (N − 1). To keep this simple form of contractions also when we work with the state f GP , we construct another onebody basis set: We take f GP as the zeroth element of our basis, and use the Gram-Schmidt process to build an orthogonal basis set from f GP , φ 1 , φ 2 , ....

A.3 Truncation of the Hilbert space
We truncate the one-body basis to solve the flow equations numerically. Let us denote the size of the truncated basis by n. As argued in Ref. [1], to calculate the energy of the system in the limit n → ∞, one should compute for a few values of n, and then fit the obtained sequence using the function where (n → ∞), A, δ are fit parameters. (n → ∞) is the value of the energy in the limit n → ∞. This value is presented in the figures of the main text. In our calculations, we fit results for n = 13 − 23. We estimate the uncertainty by the standard deviation error of the fit. We show convergence of other observables in Figs. 14 b), c), and d).
Observables that depend on the density of the bosons in the vicinity of the impurity [see Figs. 14 b) and d)] approach the limit n → ∞ in a similar fashion to the energy. Such behavior is the result of a slow convergence of the wave function due to the delta-function potential. To calculate such observables in the limit n → ∞, we use the fitting procedure described above.
Observables that do not depend on the density of bosons in the vicinity of the impurity, e.g., the density of the Bose gas far away from the impurity [see Fig. 14 c)], are virtually converged for n = 23. Therefore, to estimate the value of the observable in the limit n → ∞, we take the result for n = 23. To estimate the uncertainty of this value, we use the largest difference between the results obtained with n = 13 − 23.

B An Attractive Impurity in a Bose gas
We leave a rigorous study of a Bose gas with a single attractive impurity to a future publication [37]. However, for the sake of completeness, we briefly discuss here properties of a system defined in Sec. 3 with c < 0. In Fig. 15, we present the density of the Bose gas. The figure shows that the mean-field approximation works well for c = −0.1, and c = −0.2. However, there is a difference between the MFA and IM-SRG results for two largest values of |c|. The difference is most noticeable for the density at z = 0, which defines the contact parameter  Finally, we calculate phase fluctuations, see Fig. 15 e). The figure demonstrated significant phase fluctuations for c 0.25. This increase signals that the Bose gas and the impurity form a many-body bound state. Large phase fluctuations also suggest that the mean-field ansatz is no longer a valid approximation.