Induced correlations between impurities in a one-dimensional quenched Bose gas

S. I. Mistakidis ,1 A. G. Volosniev,2,3 and P. Schmelcher1,4 1Center for Optical Quantum Technologies, Department of Physics, University of Hamburg, Luruper Chaussee 149, 22761 Hamburg Germany 2Institut für Kernphysik, Technische Universität Darmstadt, 64289 Darmstadt, Germany 3IST Austria, am Campus 1, 3400 Klosterneuburg, Austria 4The Hamburg Centre for Ultrafast Imaging, Universität Hamburg, Luruper Chaussee 149, 22761 Hamburg, Germany

Let us briefly address what is known about impurityimpurity interactions induced by a weakly interacting Bose gas in one spatial dimension see Refs. [18][19][20][21] for recent studies on two impurities in cold three-dimensional Bose gases, and Refs. [22][23][24][25][26] for a discussion of impurities in strongly interacting Bose gases. Homogeneous systems with weakly interacting and heavy impurities are well understood by now. The interaction at short distances follows an attractive exponential function [27][28][29][30]. At long distances, quantum fluctuations lead to a power-law decay of interactions [26,30,31]. One can argue that weakly interacting mobile impurities interact similarly [26]. To observe induced interactions in cold-atom systems, one should understand finite-size as well as the trap effects. A number of works studied these effects using time-independent models at various parameter regimes [25,29,[32][33][34].
Our work is in line with the previous studies on induced interactions in harmonically trapped static systems. Moreover, it extends them by investigating a corresponding dynamical problem. There are two main motivations for this work. First, we aim at understanding for which parameters and observables an effective interaction derived for a homogeneous medium can be used to adequately describe the dynamics in a harmonic trap. Note that the induced interaction in a harmonic trap is not Galilean invariant [29], which complicates the use of a Galilean invariant effective potential characteristic for a homogeneous environment. As we show, this complication can be avoided by considering weakly interacting systems. Second, we want to estimate the effect of the induced interactions on the dynamics.
In this paper, we introduce an effective model for impurityimpurity interactions and compare it with the multilayer multiconfiguration time-dependent Hartree method for bosons (ML-MCTDHB) [35][36][37][38][39][40][41]. It is shown that a two-body effective model with a zero-range potential is able to explain qualitatively the time evolution of two impurities that interact weakly with a trapped Bose gas. For moderate interactions, we observe that the impurity-impurity interaction cannot be modeled using a delta function. Our findings demonstrate that even weak interactions noticeably affect the dynamics. We conclude that the quench dynamics may be used to observe induced interactions experimentally.
The paper is organized as follows. Section II formulates the time-dependent problem we consider: a trapped Bose gas with two impurity atoms being initially in the ground state for vanishing boson-impurity interactions. The time dynamics is then initiated by a sudden change of the boson-impurity interaction strength. Section III introduces an effective twobody Hamiltonian that we use to study the time evolution of the impurity-impurity subsystem. Section IV presents the time evolution of the size of the impurity cloud, the entropy, and the two-body coherence function. In that section, we compare the numerical data with the effective model from Sec. III. Section V summarizes our results and provides an outlook. Three Appendices contain the technical details of our work. Appendix A reviews the ML-MCTDHB method. Appendix B discusses the accuracy of the numerical data, and Appendix C elaborates on the effective two-body Hamiltonian.

II. HAMILTONIAN
We consider two bosonic impurities in a system of N bosons. The corresponding Hamiltonian reads as δ(x i − y j ) + g II δ(y 1 − y 2 ), (1) where y j is the position of the jth impurity, and x i is the position of the ith boson. The parameters g IB , g BB , and g II determine the strengths of the zero-range interactions. Without interactions, all particles are described by the identical one-body Hamiltonians h(z) = −¯h 2 2m ∂ 2 ∂z 2 + kz 2 2 . Here, m is the mass of a particle, k = m 2 is the spring constant, and is the frequency of the external harmonic trap. The groundstate properties of the Hamiltonian (1) were studied using a variational wave function in Ref. [29] (see also [42]). In this work, we focus on the dynamics following a sudden change of g IB (g IB = 0 at t = 0, g IB = 0 at t > 0) assuming that the system is in the ground state at t = 0. The time evolution of the system obeys the time-dependent Schrödinger equation (2) where the many-body wave function ψ satisfies the initial condition: ψ (t = 0) is the ground state of H (g IB = 0). We tackle Eq. (2) using the ML-MCTDHB method [35][36][37] (see Appendices A and B for a brief description of the method and the accuracy of the calculations). ML-MCTDHB is a variational method for the many-body dynamics that allows us to numerically study beyond mean-field correlations between particles. In the ML-MCTDHB approach, the many-body wave function is expanded in a time-dependent and variationally optimized basis that spans the most important part of the Hilbert space and disregards the remainder. To construct this basis, we first write the many-body wave function as a truncated Schmidt decomposition using D species functions for each component [see also Eq. (A1) in Appendix A]. As a next step, we expand the aforementioned species functions in a basis of d B (d I ) single-particle functions for the bosonic medium (the impurities) [see also Eq. (A2) in Appendix A]. Then, a variational principle leads to a set of nonlinear integrodifferential equations. For convenience, the single-particle functions are expanded with respect to a time-independent primitive basis. In our case, this primitive basis is given by a sine discrete variable representation of a one-body space with hard-wall boundary conditions at both edges of the numerical grid. The grid consists of 500 points.
Our system constitutes of a Bose gas with many weakly interacting particles and two impurity atoms. A large number of bosons can be adequately described using a small number of single-particle functions. In our calculations, we established that d B = 3 captures the dynamics well. The number of impurities is small and we need to use more single-particle functions, d I 8, to obtain accurate results. It is worthwhile noting that the variational nature of ML-MCTDHB allows us to estimate its accuracy by varying D, d I , and d B . For details on the precision of our numerical results, see Appendix B.
The dynamics of a single impurity has already been explored using the ML-MCTDHB [43][44][45]. In this paper, we seek insight into induced impurity-impurity correlations. We compare the dynamics of two impurities in a Bose gas to that of a single impurity presented in Ref. [43] (see also Fig. 1). Therefore, for our numerical simulations, we use the parameters of Ref. [43], namely, N = 100, = 2π × 20Hz, g BB = 10 −37 Jm, and m = m( 87 Rb). To give a sense of the size of the Bose gas for these parameters, we note that the corresponding Thomas-Fermi radius is approximately 19 μm.
To illustrate the effect of induced interactions in our model, we determine the square of the size of the cloud of two impurities placed in a Bose gas: where MB is obtained using the ML-MCTDHB method (see Appendix A). To extract information about induced impurityimpurity correlations, the quantity y 2 1 + y 2 2 MB is compared to the square of the size of the cloud of two impurities in two separate Bose gases [see Fig. 1(a)]: where ψ 1 is the wave function of a Bose gas with only one impurity atom (see Ref. [43]). We present our results in Fig. 1(b) for g II = 0 and g IB (t > 0)/g BB = 0.3. Figure 1(b) demonstrates that following the sudden change of g IB the square of the size of the impurity cloud determined by Eq. (3) (correlated impurities) features an oscillation motion similar to that of Eq. (4) (noncorrelated impurities). However, y 2 1 + y 2 2 1,MB (t ) is noticeably different from y 2 1 + y 2 2 MB (t ) [note that g IB (t = 0) = 0, hence, y 2 1 + y 2 2 MB (t = 0) = y 2 1 + y 2 2 1,MB (t = 0)]. To interpret Fig. 1, we introduce in the next Illustrates the two systems of interest: (left) a Bose gas with two impurity atoms, and (right) two Bose gases, each with a single-impurity atom. (b) Shows the time evolution of y 2 1 + y 2 2 (t )/ y 2 1 + y 2 2 (0) after a rapid change of the boson-impurity interaction strengths. The (red) dashed curve describes the impurities in different traps [43]. The (black) solid curve presents the dynamics of two impurities in the same trap. There is no free-space impurity-impurity interaction g II = 0, and we interpret the difference between the two curves as a manifestation of attractive induced impurity-impurity interactions (see the text for details). The final interaction strength is g IB (t > 0)/g BB = 0.3, all other parameters (e.g., the number of bosons N = 100) are common for all numerical data and are given in the text.
section an effective description based on dressed impurities. According to that model, the difference between the curves in Fig. 1(b) is a manifestation of induced impurity-impurity interactions. Below, we discuss this effective description, which allows us to visualize and quantify the induced impurityimpurity correlations. This description is a useful tool that gives insight into systems with multiple impurities and does not rely on elaborate many-body simulations.

III. EFFECTIVE TWO-BODY HAMILTONIAN
The simplest phenomenological approach to a single impurity in a Bose gas is the use of an effective one-body Hamiltonian. In a homogeneous case, such a Hamiltonian is written as follows: where (g II , g IB , ρ) and m eff (g II , g IB , ρ) are, respectively, the self-energy and the effective mass of the dressed impurity atom [8,43,[46][47][48][49][50], ρ is the density of the homogeneous Bose gas. In this paper, we are interested in impurity-impurity correlations in a trap where the momenta of the impurities are nonvanishing. Therefore, the emergent correlations between impurities can be considered as a perturbation, which can be accounted for by some potential V that mimics the relevant low-energy transitions between the states of h 1 . We choose the two-body effective Hamiltonian for a homogeneous case k = 0 as Our focus is on weak correlations for which the exact shape of V is not important. For simplicity, we assume that where the parameter g i II defines the strength of the induced interactions.
Trapped case. The Hamiltonian (6) is not applicable to the experiments with trapped cold atoms. To extend Eq. (6) to a harmonically trapped case, we assume that the density of the Bose gas varies slowly on the length scale given by the healing length. In this case, we can rely on the local density approximation to derive the effective two-body Hamiltonian for a trapped system where we arbitrarily choose to write g i II (ρ(y 1 ), g IB , g BB ) instead of g i II (ρ(y 2 ), g IB , g BB ); both expressions are identical because the interaction term in Eq. (8) is nonvanishing only if y 1 = y 2 . Note that the effective interaction in Eq. (8) depends, in general, on the coordinates y 1 and y 2 (cf. Refs. [29,33]), and not on y 1 − y 2 as in the homogeneous case. To simplify the description, we assume that g i II (ρ(y 1 ), g IB , g BB ) = g i II (ρ(0), g IB , g BB ). This assumption is justified only for small values of g IB , for which the impurities move close to the center of the trap, and, therefore, experience only a weak inhomogeneity during the time evolution.
It has been argued that (y i ) and m eff (y i ) in Eq. (8) are very sensitive to the density of bosons ρ(y) [43]. To circumvent this problem, another one-body effective Hamiltonian has been proposed [43], to analyze a single-impurity atom in a trapped Bose gas. The parameters k and m eff incorporate the effects of m eff (y) and (y). They can be obtained by fitting to numerical and/or experimental data. We determine them as in Ref. [43], i.e., we first use the ML-MCTDHB to calculate y 2 for a single impurity in a Bose gas, and then find the parameters in Eq. (9) that most accurately reproduce the ML-MCTDHB results. For g IB > 0, k and m eff are presented in Ref. [43]. Note that h trap 1 is accurate only for weak interactions. For strong interactions, the impurity minimizes the interaction with the bosons by residing at the edge of the Bose gas, which implies that a harmonic term ky 2 /2 cannot be used in Eq. (9) [29,43,44,51,52]. In this paper, we focus exclusively on the regime where Eq. (9) holds.
We employ Eq. (9) to mimic the effect of single-body terms in Eq. (8), which leads to the two-body Hamiltonian that we use to analyze the dynamics of two weakly interacting impurities. We denote by eff (y 1 , y 2 ; t ) the wave function that describes the time evolution governed by H eff . By assumption, eff (y 1 , y 2 ; t = 0) is the ground state of h trap 1 (y 1 ) + h trap 2 (y 2 ). To determine eff , we first separate the relative and the centerof-mass motions. To this end, we employ the coordinates y = (y 1 − y 2 )/ √ 2 and y c.m. = (y 1 + y 2 )/ √ 2. The center-of-mass dynamics is independent of the induced interactions, hence, it follows from Ref. [43]. The relative motion is described by the Hamiltonian It is worthwhile noting that the wave function that describes the quench dynamics of the Hamiltonian (11) following the rapid change of parameters at t = 0: m → m eff , k → k, and g i II = 0 → g i II = 0 for g II = 0 can be written in a closed form (see Appendix C). This allows us to calculate directly all observables of interest, e.g., We derive an estimate for the parameter g i II , employing the interaction between two weakly interacting impurities, V HH [27][28][29][30]. The potential V HH has a short-range (mean-field) part given by the Yukawa-type potential (∼e −r/r 0 ) and a long-range tail (1/r 3 ) due to quantum fluctuations. The long-range tail for our parameters is important at and beyond r LR 2 μm. We estimate this distance by comparing the value of the interaction due to the long-range part to that given by the short-range part [30] where ξ 0.4 μm is the healing length estimated using the density of the Bose gas at the origin ρ(0). Note that r LR 5ξ , which implies that the effect of the long-range tail can be neglected in our analysis. In particular, we estimate the potential volume of V HH using only its short-range part because a longrange part has a negligible potential volume. The short-range part of the potential can be derived from the density profile of a homogeneous Bose gas with a single weakly interacting impurity [48] ρ where ρ 0 is the density of the Bose gas without the impurity, and κ is the boson-impurity reduced mass. The mean-field energy of the second impurity is given by g IB ρ(y 2 ). The induced impurity-impurity interaction is determined by the second term of Eq. (13). The corresponding potential net volume is V HH (y)dy −g 2 IB /g BB . Note that it is independent of the mass of the impurity, and thus can be derived using heavy impurities. Since the potential volume is the only relevant parameter for low-energy scattering, we arrive at the expression for the strength of the induced impurity-impurity interaction Limits of applicability. The description based on Eqs. (11) and (14) fully determines the dynamics of two impurities in a trapped Bose gas. Let us summarize the crucial assumptions we use to derive this description: (i) the energy scale associated with the induced impurity-impurity interaction is much smaller than all other energy scales in the problem (e.g.,h ) so that this interaction can be parametrized by the delta function potential; (ii) the impurities move close to the center of the trap, such that the interaction depends only on the relative distance y; (iii) the interaction potential is independent of the effective mass and the spring constant, such that Eq. (14) can be used. One can argue that these assumptions are valid as long as g IB is small enough. A demonstration of this is presented in the next section, where the results of Eqs. (11) and (14) are compared to the results obtained within the ML-MCTDHB simulations.

A. Size of the impurity-impurity subsystem
For the sake of discussion, we consider only systems with g II = 0, which enjoy the most direct and clear manifestation of induced interactions. Indeed, in this case all observed correlations are induced. We determine y 2 1 + y 2 2 MB using the ML-MCTDHB method and compare it with y 2 + y 2 c.m. eff from the effective Hamiltonian [Eq. (11)]. To isolate the effect of the induced interactions, we exclude the contribution of the center-of-mass motion and compare the time evolution of y 2 1 + y 2 2 MB − y c.m.
2 eff and y 2 eff . The former quantity describes the dynamics of the relative distance between the two impurities only approximately because the center-of-mass and the relative motions of the impurities are coupled in the many-body calculations via interactions with the bath. The parameter g i II for the effective potential is obtained by fitting to the corresponding ML-MCTDHB data, and then compared to Eq. (14). To minimize the role of the effects arising due to the in-homogeneity of the Bose gas, we only fit to the data that describe the first oscillation. For consistency, the effective Hamiltonian is used throughout this work to describe the time evolution only during the first oscillation period of y 2 eff [e.g., 28 ms for Fig. 2(a)]. The discussion does not change qualitatively if we use for the fitting the first two or even three oscillation periods.
Our results are presented in Fig. 2. First of all, we observe that the values of g i II used in Figs. 2(a) and 2(b) are consistent with the estimate of Eq. (14). For the case g IB = 0.1g BB , Eq. (14) suggests that g i II −0.01g BB ; for g IB = 0.3g BB it implies g i II −0.09g BB . These estimates are in a good agreement with the values obtained by fitting to the ML-MCTDHB data, namely, g i II = −0.011g BB and g i II = −0.10g BB . We observe that even though the parameter g i II is small it significantly suppresses the amplitude of the oscillations (compare the dashed and solid curves in Fig. 2). Therefore, the induced The square of the relative distance between the impurities as a function of time for repulsive impurity-boson interactions. We choose to normalize the distance to one at t = 0. The (red) dashed curve in every panel shows [ y 2 2 eff ](t = 0) for two uncorrelated impurities [43]. The (black) solid curve demonstrates the ML-MCTDHB results for two correlated impurities [ y 2 correlations could be studied in an experiment by monitoring the quench dynamics. The challenge here is to create the initial state ψ (t = 0) at sufficiently low temperatures (see also Sec. V). We also explore the change of the bosonimpurity interaction strength to negative values of g IB , i.e., 2 eff ](t = 0) for two uncorrelated impurities [43]. The (black) solid curve demonstrates the ML-MCTDHB results for two correlated impurities [ y 2 2 eff ](t = 0). The dots in the panels show the results of the effective Hamiltonian (11) g IB (t > 0) < 0. According to Eq. (14), the strength of induced correlations is independent of the sign of g IB , if g BB and |g IB | are small. The ML-MCTDHB simulations agree with this conclusion (see Fig. 3). The dots in Fig. 3 show not a fit, but a prediction based on g i II obtained for g IB > 0, e.g, we use in the effective model g i II = −0.1g BB fitted above for g IB (t > 0) = 0.3g BB to obtain the dynamics following the change to g IB (t > 0) = −0.3g BB . The parameters m eff and k for Eqs. (9) and (11) are obtained as in Ref. [43], by fitting to the ML-MCTDHB results for a single impurity; these parameters are m eff = 1.02m, k = 1.14k for g IB (t > 0) = −0.1g BB , and m eff = 1.075m, k = 1.46k for g IB (t > 0) = −0.3g BB . All in all, the results shown in Figs. 3(a) and 3(b) confirm our expectations based on Fig. 2: the relative distance for two weakly correlated impurities is smaller than the distance between two uncorrelated impurities, which we attribute to the attractive interaction mediated by the Bose gas.
For larger values of |g IB | the quench dynamics for g IB > 0 cannot be compared with that for g IB < 0. For example, for g IB = g BB , the impurities are pushed to the edge of the trap, whereas for g IB = −g BB , the impurities move in the vicinity of the center of the trap [see Fig. 3(c)]. For comparison, we present in Fig. 3(c) also the results of the effective model with m eff = 0.83m, k = 2.56k, and g i II = −0.74g BB . As before, the parameters m eff and k are obtained from the dynamics of a single impurity in a Bose gas. The parameter g i II is now obtained by fitting to the ML-MCTDHB results for two impurities.
We remark that already for g IB = 0.5g BB the data cannot be accurately fitted with the proposed effective model. For the sake of discussion, we show the result of the effective model with g i II = −0.25g BB , which is expected from Eq. (14), in Fig. 2(c). Our failure to fit the data is not surprising: (i) The impurity starts to probe a large part of the Bose cloud [see the amplitude of oscillation of the noncorrelated impurities in Fig. 2(c)], which means that our assumption that V depends only on y 1 − y 2 must be modified. (ii) The zero-range approximation to the effective potential is also not correct, impurities start to resolve the shape of the induced potential. Some of these effects might be included by using a more elaborate form of the effective potential, as, e.g., in Refs. [29,33]. This discussion is beyond the scope of this paper, therefore, we refrain from discussing cases with g IB > 0.3g BB further.

B. Entropy
Induced impurity-impurity correlations can be studied using the one-body reduced density matrix Let us denote the eigenvalues of ρ (I ) MB as n i (t ). We arrange them in ascending order, such that n 1 n 2 n 3 . . .; n i (t ) are often referred to as natural occupation numbers. They can be seen as probabilities for finding an impurity in a given state since i n i = 1. It is clear that the impurities are correlated if two or more n i are nonzero, otherwise, the wave function MB is just a product state. To quantify the strength of impurityimpurity correlations, we introduce the entropy If n 1 (t ) = 1 and n j =1 (t ) = 0, then S(t ) = 0, signaling the absence of impurity-impurity correlations. Other possibilities lead to nonzero values of S(t ). Since the two bosonic impurities are noninteracting in free space, any positive value of entropy can be used as a witness of induced impurity-impurity correlations mediated by the bosonic gas. We first calculate the time evolution of S(t ) using the ML-MCTDHB method for different values of g IB for t > 0 (see Fig. 4). Initially, the impurities do not interact, therefore, S(t = 0) = 0. For t > 0 we observe that S = 0, which indicates the presence of impurity-impurity correlations. For g IB (t > 0) = 0.1, the entropy is close to zero, meaning that the mean-field treatment can be used to decently describe the dynamics for such weak interactions. For larger values of g IB (t > 0), beyond-mean-field corrections are important and must be taken into account (cf. Appendix A, where mean-field calculations of the size of the impurity cloud are presented). At the early stages of the dynamics, S(t ) exhibits a sudden linear increase, while, for later time instants, it shows an oscillatory behavior possessing a multitude of frequencies. The initial increase of S(t ) is more pronounced for larger postquench interspecies interaction strengths, e.g., S I (t = 3.3) = 0.032 for g IB = 0.3g BB , and S(t = 3.3) = 0.079 for g IB = 0.5g BB , suggesting stronger impurity-impurity correlations for larger values of g IB . Also, the average degree of impurity-impurity correlations quantified byS(T ) = (1/T ) T 0 dt S(t ) is larger for an increasing postquench g IB , for instance,S(120 ms) ≈ 0.027 at g IB = 0.3g BB , andS(120 ms) ≈ 0.082 at g IB = 0.5g BB . We conclude that a larger postquench g IB implies a larger degree of impurity-impurity correlations (cf. [41,44,53]).
Next, we calculate the entropy using the effective model [Eq. (11)]. To this end, we diagonalize the one-body density matrix ρ (I ) eff (y 1 , y 1 ; t ) = dy 2 eff (y 1 , y 2 ; t ) * eff (y 1 , y 2 ; t ), (17) whose eigenvalues are then used in Eq. (16) to determine the entropy (see Appendix C for details). We show this entropy for g i II = −0.10g BB [g IB (t > 0) = 0.3g BB ] in Fig. 4 for up to t 30 ms, which corresponds to one oscillation period [cf. Fig. 2(b)]. The spectrum of the effective model is close to that of a harmonic oscillator. It leads to an almost complete restoration of the initial wave packet, and, hence, S(t ) approaches a minimal value after one oscillation period. In contrast, the ML-MCTDHB results suggest a fast approach to a thermal-like state, where the entropy oscillates around its average value. However, the prediction of the ML-MCTDHB and the effective model agree on the average values of the entropies. The ML-MCTDHB method yields S(30 ms) 0.025 for g IB = 0.3g BB , which should be compared with S(30 ms) 0.029 of the effective model. We conclude that the effective model is useful to calculate average values of the entropy. We checked that this conclusion holds also for the dynamics initiated by a change to attractive boson-impurity interactions, i.e., if g IB (t > 0) < 0.
To give insight into the difference between the entropy calculated using the many-body ML-MCTDHB approach and the effective model, we note that the dynamics in the effective model is driven by the harmonic oscillator whose equidistant spectrum gives a single oscillation frequency. In reality, there are a multitude of processes associated with the impuritybath interactions, which affect the equidistant spectrum of the harmonic oscillator. The ML-MCTDHB calculations (unlike the effective model) capture these processes, which lead to the observed thermal-like state. In the future, it will be interesting to study possible extensions of the effective model that can lead to the observed approach to a thermal-like state.

C. Two-body correlation function
To reveal impurity-impurity correlations in a spatially resolved manner, we study the normalized two-body correlation function where ρ (II ) MB (y 1 , y 2 ; t ) is the two-body reduced density matrix The function ρ (II ) MB (y 1 , y 2 ; t ) represents the probability of measuring at time t one impurity at y 1 and another at y 2 . The two impurities correlate if g (II ) MB (y 1 , y 2 ; t ) = 1; g (II ) MB (y 1 , y 2 ) > 1 [g (II ) MB (y 1 , y 2 ) < 1] implies an increased (decreased) probability to observe two particles at y 1 and y 2 in comparison to two uncorrelated impurities. If g (II ) MB (y 1 , y 2 ; t ) = 1 the impurities are uncorrelated. Figure 5 shows g (II ) MB (y 1 , y 2 ; t ) for the change of the interaction strength from g IB = 0 to g IB = 0.3g BB . Initially, at t = 0, the two noninteracting bosonic impurities are uncorrelated which leads to g (II ) MB (y 1 , y 2 ; t = 0) = 1 [ Fig. 5(a)]. Two-body correlations between the impurities begin to develop at t > 0  Fig. 2(b)]. This breathing dynamics leads to the expansion and contraction of ρ (I ) (y 1 , y 1 ; t ), and to two distinct correlation patterns present in g (II ) MB . We exemplify them in Figs. 5(c) and 5(f). In Fig. 5(c), the two impurities attract each other and move together to the edge of the MB of the two bosonic impurities at different time instants in ms (see legends) during the quench dynamics that follows the change of the boson-impurity interaction strength: MB (y 1 , y 2 = y 1 ; t = 20 ms) 1. In Fig. 5(f), the impurities also move toward the edge of the trap. However, they move not only as a pair but also separately because both g (II ) MB (y 1 , y 2 = y 1 ; t = 40 ms) and g (II ) MB (y 1 , y 2 = −y 1 ; t = 40 ms) are greater than (or equal to) one.
To better understand impurity-impurity correlations depicted in Fig. 5, we use the effective model to calculate the two-body correlation function The time evolution of g (II ) eff is depicted in Fig. 6 for g IB (t > 0) = 0.3g BB . We observe a qualitative agreement between g (II ) eff [Figs. 6(b) and 6(e)] features two distinct patterns that are characterized by g (II ) eff > 1 along the diagonal and g (II ) eff > 1 across the antidiagonal. According to the effective model, these patterns represent the motion of a bound pair and the backward scattering of two dressed impurity atoms, respectively.

V. SUMMARY AND OUTLOOK
We explore the time evolution of a Bose gas with two impurities initiated by a sudden change of the boson-impurity interaction strength. Our focus is on the dynamics of the impurities, which we analyze using the many-body correlated ML-MCTDHB method [38,39] and an effective model presented in Eq. (10). We observe that the Bose gas induces impurity-impurity correlations. The strength of the induced interactions can be estimated from the attractive Yukawa-type potential between two heavy impurities in a homogeneous Bose gas [27][28][29][30]. This potential captures well the dynamics of the size of the impurity cloud. It also explains qualitatively the time-averaged value of the entropy and the correlation patterns that appear in the two-body correlation function. We observe a pronounced effect of the induced correlations on the dynamics. This means that the quench dynamics may become a test ground for studies on induced interactions provided that the initial state, the ground state of the system, can be prepared.
It will be interesting in the future to explore other types of induced interactions. For example, the long-range interactions (1/r 3 ) due to quantum fluctuations are too weak in our system, but can become important if two impurities are confined by traps whose origins are well separated. Another exciting direction is systems with repulsive induced interactions. Repulsive interactions can appear if two impurities are distinguishable with one impurity repelling the Bose gas, and the second impurity attracting it. Then, according to Eq. (13), the impurity-impurity interaction is repulsive; its strength is given by g i II = −g 1 IB g 2 IB /g BB , here g l IB characterizes the interaction between a boson and the lth impurity.
Future studies are needed to understand the effect of a finite temperature T on the discussed dynamics. The energy scale associated with the induced interaction is small (see the estimate below) in comparison to typical energy scales in current cold-atom experiments. Therefore, it is natural to expect that the induced correlations are important only at very low temperatures (cf. Refs. [54,55]). To estimate the relevant temperature scales, we analyze a dimensionless parameter constructed from quantities that enter Eq. (11): s = |m(g (i) II + g II )λ typ /(h 2 )|, where λ typ is a typical length scale. The parameter s is the only dimensionless parameter for a homogeneous system; it also determines the relative strength of the induced interactions in shallow traps. If s is large, then the induced interactions are important; in the opposite limit s → 0, the impurities correlate weakly. If we use a thermal wave length as a typical length scale, i.e., λ typ h/ √ 2mk B T (k B is the Boltzmann constant), then the parameter s is close to unity for T nK, assuming that g II + g i II are in the range −(10 −38 -10 −37 )Jm, and m = m( 87 Rb). For higher temperatures the induced interactions cannot affect the dynamics since s → 0. Future studies are needed to understand the transition from small to large values of s.
In this work, we focus on the case with g II = 0. This is an ideal scenario, which might be difficult to realize experimentally. In the future, it will be interesting to identify the most promising systems in which induced interactions can be observed even if g II = 0. Note that the limit 1/g II → 0, which corresponds to fermionic impurities, has already been studied [25,29,53]. This limit features much weaker effect of the induced interactions on the properties of the system.
In this paper we focus on the dynamics of the impurities. However, the ML-MCTDHB method also gives access to the dynamics of the Bose gas, which, for weak interactions considered here, performs a small-amplitude breathing motion at t > 0. In general, we expect a number of interesting phenomena associated with the Bose gas. For example, it is known that switching off boson-impurity interactions can lead to shock waves and solitons in the Bose gas [56]. In light of the importance of beyond-mean-field effects in our study (see Appendix A), it will be interesting to study these nonlinear objects for moderate interactions using the ML-MCTDHB method [38,39]. S.I.M. and A.G.V. contributed equally to this work.

APPENDIX A: A BRIEF DESCRIPTION OF THE ML-MCTDHB METHOD
To investigate the stationary properties of the Hamiltonian in Eq. (1), and the quench dynamics we address numerically the underlying many-body Schrödinger equation. We employ the ML-MCTDHB method [35][36][37], which is a variational approach for solving the time-dependent manybody Schrödinger equation for atomic mixtures composed either of bosonic [38,39,43,44] or fermionic [40,41,53,[57][58][59] species. The method relies on the expansion of the many-body wave function in terms of a time-dependent and variationally optimized basis enabling us to take into account both the interspecies and the intraspecies correlations of two-component systems. In this Appendix, we briefly review the method for the convenience of the reader.
The many-body wave function is first expressed as a superposition of D different species functions for each of the species, i.e., B k ( x; t ) and I k ( y; t ). Here, x = (x 1 , . . . , x N ) and y = (y 1 , y 2 ) are the spatial coordinates of the N bosons and the two impurities, respectively. Consequently, the many-body wave function MB (t ) includes interspecies correlations and has the form of a truncated Schmidt decomposition [60] with rank D, namely, The Schmidt coefficients λ k (t ) are known as the natural species populations of the kth species function [37,38,44].
The system is said to be entangled [61] or interspecies correlated when at least two distinct λ k (t ) possess a nonzero value. In this case MB cannot be expressed as a direct product of two states. Note also that λ k (t ) are the eigenvalues of the species reduced density matrices, e.g.,  Additionally, the SPFs are expanded on a time-independent primitive basis {|q } being in our case an M-dimensional discrete variable representation (see also below). Note that the eigenfunctions of the one-body reduced density matrix for bosons and impurities, i.e., ρ (1) σ (z, z ; t ) = MB (t )|ˆ σ † (z)ˆ σ (z )| MB (t ) (σ = {B, I}, z = {x, y}), whereˆ σ (z) is the bosonic field operator, are the so-called natural orbitals φ σ i (z; t ). The natural orbitals are related with the SPFs via a unitary transformation that diagonalizes ρ (1) σ (z, z ; t ) when it is expressed in the basis of SPFs. For a more elaborate discussion, see Refs. [35][36][37]. The eigenvalues of ρ (1) σ (z, z ; t ) are the natural populations n σ i (t ). They provide a theoretical measure for the degree of intraspecies correlations. When there is more than one macroscopically occupied state, the σ subsystem is called intraspecies correlated. Otherwise, it is fully coherent.
Having at hand the many-body wave-function ansatz, we can establish the ML-MCTDHB equations of motion [37]. To this end, we apply the Dirac-Frenkel variational principle [62,63] to the function determined by Eqs. (A1) and (A2). As a result, we obtain a set of D 2 linear differential equations of motion for the Schmidt coefficients λ k (t ) which are coupled to nonlinear integrodifferential equations for the species functions and d B + d I integrodifferential equations for the SPFs. We note in passing that the ML-MCTDHB method allows us to operate within different approximation orders. As an example, we retrieve the mean-field equation of motion [64,65] of the bosonic mixture when choosing D = d B = d I = 1. In this case, the many-body wave-function ansatz reduces to the well-known mean-field product state Following a variational principle for this ansatz, e.g., due to Dirac and Frenkel [62,63], we arrive at the celebrated system of coupled Gross-Pitaevskii equations of motion [64,65] that govern the dynamics of the bosonic mixture where we introduce N I = 2 to elucidate the symmetry of the equations. Within the mean-field approximation, all particle correlations of the system are neglected, implying that it can be useful only to describe the system at very weak interactions. To illustrate this, we employ the coupled Gross-Pitaevskii equations to study the time evolution of the size of the impurity cloud following a sudden change of g IB . We calculate (y 2 1 + y 2 2 )| MF | 2 dy 1 . . . dx N and 2 y 2 MF is the mean-field wave function for a single impurity in a Bose gas [43]. Following the discussion in the main text, we present in Fig. 7 the quantity y 2 MF (t )/ y 2 MF (0), where y 2 (t ) MF = (y 2 1 + y 2 2 )| MF | 2 dy 1 . . . dx N − y 2 1 | (1) MF | 2 dy 2 . . . dx N . We conclude that the parameter g i II obtained from the ML-MCTDHB calculations agrees well with that for the mean-field calculations only when g IB (t > 0) is small [see Fig. 7(a)]. Already for g IB (t > 0) = 0.3, beyond-mean-field simulations differ from the mean-field predictions [see Fig. 7(b)].

APPENDIX B: CONVERGENCE AND DETAILS OF THE MANY-BODY SIMULATIONS
As stated in Appendix A, the ML-MCTDHB method is based on the expansion of the many-body wave function with respect to a time-dependent and variationally optimized basis. The underlying Hilbert space truncation is determined by the employed orbital configuration space which we denote as C = (D; d B ; d I ), with D and d B , d I being the number of species and single-particle functions, respectively, of each species [see also Eqs. (A1) and (A2)]. For our simulations we utilize a primitive basis based on a sine discrete variable representation including 500 grid points. This sine discrete variable representation introduces hard-wall boundary conditions imposed at x ± = ±40 μm. The presence of the hard-wall boundary does not affect our results since we do not observe any significant density population beyond x ± = ±23 μm.
The many-body calculations presented in the main text are based on the orbital configuration C = (10; 3; 8). We study the convergence of the many-body simulations by varying the orbital configuration space C = (D; d B ; d I ). Let us exemplify  [43]). The interaction parameter g i II = −0.011g BB is identical to the one used in Fig. 2(a). (b) Illustrates the dynamics with g IB (t > 0) = 0.3. The parameters for the effective Hamiltonian m eff /m = 1.045 and k eff /k = 0.76 are obtained using the methods of Ref. [43]. The interaction parameter g i II = −0.1g BB is identical to the one used in Fig. 2 the convergence of our results for a different number of species and single-particle functions. For this investigation we resort to the time evolution of the variance of the two bosonic impurities Y 2 (t ) = y 2 1 + y 2 2 . We study its absolute deviation between the C = (10; 3; 8) and other orbital configurations C = (D; d B ; d I ), namely, In this expression, Y 2 (t ) C is calculated for the orbital configuration space C. Y 2 (t ) C,C → 0 implies a negligible deviation between Y 2 (t ) calculated within the C and C approximations. Figure 8 presents Y 2 (t ) C,C for the system considered in the main text with g IB (t > 0) = 0.3 [ Fig. 8(a)] and g IB (t > 0) = 0.5 [ Fig. 8(b)]. A convergence of Y 2 (t ) C,C is seen in both cases. More specifically, comparing Y 2 (t ) C,C between the C = (10; 3; 8) and C = (12; 3; 8) orbital configurations we observe that it acquires values below 0.05% (0.1%) for postquench interactions g IB = 0.3 (g IB = 0.5) throughout the time evolution. Also, the values of Y 2 (t ) C,C , when C = (10; 3; 8) and C = (8; 4; 7), imply relatively small deviations which become at most of the order of 1.4% (1.8%) in the course of the time dynamics for g IB = 0.3 (g IB = 0.5). The same observations hold also true for the variances of the bosonic gas (not shown here for brevity). Finally, we remark that a similar degree of convergence occurs also for the other observables invoked in this paper, e.g., the single-particle density of the impurity (not shown for brevity).

APPENDIX C: WAVE FUNCTION OF THE EFFECTIVE TWO-BODY HAMILTONIAN
Here, we discuss the wave function rel (t ) that evolves according to the effective Hamiltonian introduced in Eq. (11) of the main text: The initial state for the dynamics is the ground state of the one-body Hamiltonian h(y) = − 1 2 ∂ 2 ∂y 2 + Note that the Hamiltonian (C1) preserves parity, and that rel (t = 0) is an even-parity function. Therefore, rel (t ) must where eff = √ k/m eff , U is Tricomi's confluent hypergeometric function (also known as the confluent hypergeometric function of the second kind) [67], ν i = E i 2 eff − 1 4 , and the normalization constant is The parameter ν i can be found from the equation [66] −ν i + 1 The expansion coefficients can be expressed in a closed form where 2c = 1 + m m eff eff , and 2F1 is the hypergeometric function. Once the wave function is derived, we can calculate all observables of interest. For example, the size of the cloud, y 2 = y 2 | (y, t )| 2 dy, can be computed by truncating the sum where A i j = φ i y 2 φ j dy has an analytic expression [68]. Another observable of interest is the entropy. To calculate it, we need to find the spectral representation of the one-body density matrix for a system of two impurity atoms. This density matrix is derived from the total wave function, which describes both the center of mass and the relative dynamics eff = c.m. (y c.m. , t ) rel (y, t ), where the function c.m. (y c.m. , t ) equals to φ(y c.m. , t ) [43] with φ(z, t ) = √ eff m eff i √ π sin( eff t ) .
The density matrix can be written as In this equation, { f i } is an orthonormal set of functions constructed such that the function f i>0 is orthogonal to φ, At every time instant t, the parameter δ(t ) defines a time-independent harmonic oscillator whose eigenstates are given by f i (t ).
The expansion coefficients of ρ (I ) eff read as β i,i (t ) = j α * i, j (t )α i , j (t ), where where y = (y 1 − y 2 )/ √ 2 and Y = (y 1 + y 2 )/ √ 2. The expression for α i, j can be simplified as follows: 1 2 δ(t ) −1/4 √ − im eff eff cot( eff t ) rel (y, t ) f * i+ j (y, t )dy Note that α i, j = 0 if i + j is an odd number. The same is true for β i, j . To obtain the spectral representation of the density matrix, which allows us to calculate the entropy, we find and diagonalize the matrix β i,i .