We thank the Referees for taking the time to review our paper. We are happy to see the overall positive evaluation of our paper. The comments in the reports helped us to significantly improve our manuscript. We hope that the revised version is ready for publication in SciPost Physics.

Below, we provide a point-by-point reply to the comments in the reports. For convenience, we quote the comments of the Referees. The main changes in the manuscript are shown in blue. Minor changes are not highlighted in the manuscript.

Reply to Anonymous Report 1:

Major comments:

1. Referee: "The manuscript contains a large number of figures and many of them (Figs 1, 2, 6, etc) are not readable in black and white version. I suggest the Authors modify the figures and check if they are easy to read in a printed version."

Our reply: We changed all figures which were not readable in black and white. We modified the captions accordingly.

2. Referee: "I think that notation f_1 and f_2 is not as obvious as it would be with f_{2b} and f_{GP} denoting two-body and Gross-Pitaevskii solutions used as input for generating the correlations."

Our reply: We thank the Referee for this useful suggestion. We have modified the text and figures accordingly. Note that we chose to use f_{1b} instead of the suggested f_{2b} to emphasize that the reference state is the solution of a one-boson Hamiltonian.

3. Referee: "The ratios “g/ρ” and “c/ρ” which are used to quantify the interaction strength are used as dimensionless parameters, while this is not the case. Here “g” has units of the coupling constant and “ρ” of the density."

Our reply: Please note that we are using the system of units in which hbar=M=1. This implies that g/rho and c/rho are dimensionless.

We added footnote 4 to clarify this point:

"In general, the dimensionless Lieb-Liniger parameter is defined as gamma=Mg/(hbar^2rho), which leads to gamma=g/rho in our units (hbar=M=1)."

4. Referee: "“One can show that the pair correlation function of the Lieb-Liniger model, g2, is identical to the density of the bosons…” I disagree, for example, for c= ∞ and γ=∞, the pair correlation function g2 = 1 - sin(πρz) ²/ (πρz) ² has oscillations which are always below the asymptotic density while the density of bosons has oscillations which exceed the bulk density."

Our reply: We thank the Referee for alerting us of this potentially confusing statement. Please note that the density of bosons is measured in the frame co-moving with the impurity. We have modified the text to stress this point and to clarify the relation between pair correlation function and density in the co-moving frame, see the discussion around Eqs. (17-19) of the revised manuscript. We have also identified a typo in our formula. The revised text presents the correct expression.

5. Referee: "“In our studies, we have noticed that the reference state f2 allows us to investigate a larger range of parameters…” Please specify what exactly is meant by “allowing to investigate”"

Our reply: We thank the Referee for identifying this imprecise statement. To clarify it, we have re-written the discussion, which now reads as:

“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.”

6. Referee: "Fig 4c, there is a non-monotonic dependence that looks rather spurious. Is it possible to reduce the errorbars?"

Our reply: Unfortunately, it is not possible to reduce the error-bars using out truncation scheme. Note that the density of bosons at the position of the impurity is small. Therefore, large error-bars in Fig. 4c is a result of working with small numbers, which cannot be avoided, unless we significantly modify our approach.

In the revised version, we have added a sentence to clarify this statement. See footnote 8: “Note that we expect that the exact curve for c/ρ= 0.5 in Fig. 4 c) is monotonous. Our calculations of thiscurve have large error-bars, which allow for an apparently non-monotonous behavior.”

7. Referee: "Discussion below Fig. 5. “For attractive interactions, the agreement is only quantitative for large values of |c|”, this statement contradicts what is shown in Fig. 5. The agreement is not perfect, but reasonable for repulsive interactions. Exactly the same can be said about the attractive case. In order to keep the claim, I guess, one has to go to a stronger attraction"

Our reply: We thank the referee for this comment. We have modified the discussion accordingly: “ 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”

8. Referee: "Please comment on the values of C3 in the fit to the energy. Is that a linear or quadratic behavior?"

Our reply: We thank the Referee for this suggestion. The revised version now states that the parameter C3 is in between 1 and 2, see footnote 11. To be more precise, we find that for c_2/rho=0.02 C3 is between 1.25 and 1.7 depending on d. For c_2/rho=0.1, C3 is between 1.1 and 2. We also checked that these windows are consistent with the value obtained by fitting to MF results. We noticed that large values of d usually imply smaller values of C3.

Minor Comments: We thank the Referee for providing us with minor comments [The revised version addresses all of them.]. They helped us to significantly improve our manuscript.

Reply to Anonymous Report 2

Referee Report 2:

1. Referee: "The flow equation results are used to benchmark the mean-field results. There is however no proof or at least some convincing arguments why the flow equation approach should be superior to the mean-field approximation. The authors compare the flow-equation results along with mean-field results with exact data obtained from Bethe ansatz in the integrable case. It would be much more convincing if a parameter case could be found, where there is a sizable difference between flow equations and mean field. Perhaps the author can comment on this."

Our reply: We thank the Referee for alerting us of this weakness in interpreting our results. Indeed, the flow equations are more accurate than the mean-field approximation. We had demonstrated this in our previous works, for example, by calculating the energy of the Lieb-Liniger model, see [1] in the revised manuscript. In that paper, for gamma=1 and Nsim 10, the flow equations yield the exact energy with 10% accuracy, the mean-field energy for the same parameters is about 50% larger.

In the revised version, we address this issue in the Introduction: “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. In the present work, we use IM-SRG to calculate the density and phase fluctuations of the Bose gas."

In addition, following the recommendation of the Referee, we illustrate a parameter regime for which there is a sizable difference between flow equations and mean field, see new Appendix B in the revised manuscript. In particular, we demonstrate that for an attractive impurity the flow equations predict significant phase fluctuations, which are beyond the mean-field approximation.

2. Referee: "The authors have calculated the impurity-impurity interaction potential. Ref. [23] predicts an exponential behavior at short distances and a power-law at large distances. The large-distance scaling is due to Casimir-like forces induced by phonon exchange and cannot be captured in the mean field approach. It would be interesting to see if the flow equation approach follows the mean field results or is actually able to reproduce the Casimir forces. Furthermore other recent work in Ref.[76] seems to suggest a linear interaction potential at short distances. Fig. 9 seems to show something different. Can the authors comment on this?"

Our reply: We absolutely agree, a numerical validation of the long-range Casimir-like force would be an extremely interesting result. Unfortunately, that force is very weak. The Casimir-like force is important only at distances of the order of 5-10 xi (xi for the healing length), which are larger than the typical sizes we consider. We have added a corresponding remark to the manuscript, see also very recent works on the topic in Refs. [76,77], where a similar conclusion is reached.

To the best of our knowledge, currently, there are no numerical techniques capable of calculating the Casimir-like force. In particular, it is also out-of-reach of the state-of-the-art Quantum Monte Carlo results [76].

Please note that 2101.11997 suggests linear interaction only for etatoinfty, i.e., not for the case presented in Fig. 9. We also expect a linear behavior when c_2 becomes very large, since our equations are identical to those of 2101.11997. The case of large c_2 is however not discussed in our work.

3. Referee: "What is the reason that the case of attractive impurity-boson interaction is not included. While one could expect that the mean-field predictions are less accurate here, in contrast to the statement just before eq.(29), the mean-field solutions can be obtained similarly to the repulsive case and it would be interesting to see if the flow equations predict different results than the mean-field approach."

Our reply: We thank the Referee for this comment. Our plan is to consider the system with attractive boson-impurity interactions in another publication. For finite systems, the attractive case is fundamentally different from the repulsive case. In particular, there is a transition from a many-body bound state to a state in which the bosons may occupy scattering states. This transition requires a separate discussion, which is beyond the scope of the present paper. Following the recommendation of the Referee, we have added Appendix B to the revised manuscript. The Appendix shows some data for a Bose gas with a single attractive impurity. In particular, it illustrates significant phase fluctuations that happen close to the aforementioned transition.

4. Referee: "In the extrapolation of the bi-polaron energies to the thermodynamic limit obtained either from flow-equations or mean-field in Fig. 9c and d it appears that the flow equations predict an oscillatory correction on top of the mean field result. This should be discussed."

Our reply: We thank the Referee for this remark. We believe that the oscillating behavior is due to the numerical accuracy of our results. In particular, the amplitude of oscillations is within the error-bars for N=70. Unfortunately, it is not easy for us to increase the accuracy for these large systems.

We have added a clarifying discussion to the revised manuscript (see page 19): “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 C1, 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.”

5. Referee: "It would be beneficial to explain in a few words what the goal of the unitary transformations governed by the flow equation (1) is. What is meant if the authors say: "to choose the operatore eta(s) .... to steer the flow in the desired direction."? What desired direction?"

Our reply: Following the recommendation of the Referee, we have added a corresponding discussion to the revised manuscript. See the discussion around Eqs. (1) and (2):

“[…]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).”

“ In our work, eta(s) is chosen from the matrix elements that describe the couplings between the ‘condensate’ and its excitations such that these couplingsbecome weaker as the flow progresses, see Fig. 1. A detailed construction of eta(s) is presented in Appendix A.”

6. Referee: "Right after eq.(2) the reference state is mentioned without any introduction. The latter is only given in the following subsection."

Our reply: We thank the Referee for this comment. We have modified the text to make the discussion more clear. In particular, we no longer mention the reference state in Chap. 2.1 “Flow equations”. Instead, we moved the discussion of our truncation scheme (in which we mentioned the reference state) to Chap. 2.2 “Reference state”.

7. Referee: "Normal ordering with respect to the reference state is mentioned but only defined in the Appendix."

Our reply: Following the recommendation of the Referee, in the revised version, we make a more clear reference to the definition of normal ordering in the Appendix.

8. Referee: "The flow equation results show error bars from "relative truncation error". Reference to Appendix A.4 should be made here to explain how the truncation error was obtained."

Our reply: We thank the Referee for this omission in the original version of our manuscript. We have added a discussion which clarifies the calculation of errorbars, and makes a more explicit reference to the Appendix.

Chapter 2.2 now starts as: “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 in second quantization using normal ordering with respect to a reference state (for a definition of normal ordering see Appendix A.1). 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. Upon normal ordering, we truncate three-body 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.2.”