Two distinguishable impurities in BEC: squeezing and entanglement of two Bose polarons

We study entanglement and squeezing of two uncoupled impurities immersed in a Bose-Einstein condensate. We treat them as two quantum Brownian particles interacting with a bath composed of the Bogoliubov modes of the condensate. The Langevin-like quantum stochastic equations derived exhibit memory effects. We study two scenarios: (i) In the absence of an external potential, we observe sudden death of entanglement; (ii) In the presence of an external harmonic potential, entanglement survives even at the asymptotic time limit. Our study considers experimentally tunable parameters.

1 Introduction Entanglement represents a necessary resource for a number of protocols in quantum information and for other quantum technologies, which are expected to be implemented in the foreseeable future for various practical applications [1][2][3][4][5].The entangled parties of a composite quantum system evolve, under realistic conditions, coupled to external degrees of freedom, which may be treated as a bath.In this open quantum system, the bath acts as a source of decoherence, leading to the destruction of quantum coherence among the states of the entangled subsystems [6].Indeed, to reach a relevant technological level, one of the main obstacles is the difficulty to ensure such a coherence despite the interaction of the system with the bath.However, the presence of the bath can produce other potentially useful phenomena.For instance, two non-interacting particles immersed in a common bath can be entangled as a consequence of an effective interaction induced by the bath [7,8].A number of situations have been considered, where indeed entanglement is observed in the aforementioned setting [9][10][11][12][13][14][15][16][17].Here we investigate the bath-induced entanglement among two distinguishable impurities embedded in a common homogeneous Bose-Einstein condensate (BEC) in 1D.Such an issue recently attracted a lot of attention, e.g, the study of impurities in double-wells in a BEC [18], the study of entanglement and the measurement of non-Markovianity between a pair of two-level localized impurities in a BEC [19,20], particle number entanglement between regions of space in a BEC [21], or the environment-induced interaction for impurities in a lattice [22,23].Moreover, in a number of experiments [24][25][26] performed this same year, entanglement between regions of a BEC was observed.In these works discrete observables were considered, while on the contrary, in our studies we will focus on entanglement of continuous variables in BEC.
The major motivation for us to undertake this work is the contemporary progress in the manipulation and control of ultracold atoms and ions, that paves the way to new possible experiments.The behaviour of an impurity in a Bose gas has recently attracted a lot of attention, both on the theoretical end experimental side [27][28][29][30][31][32][33][34][35][36][37][38][39][40][41][42][43][44][45].The main feature of such a system lies in the creation of excitation modes (Bogoliubov quasiparticles) associated to the motion of the atoms of the gas, that dress the impurity, leading to the formation of a compound system named Bose polaron.For two such impurities within a BEC, studies have focused in the past [46,47] on the possibility to form bound states (bipolarons) for sufficiently strong interactions between the impurities and the condensate atoms.Furthermore, the Bose polaron problem was recently studied within the quantum Brownian motion (QBM) model [43,48], which describes the dynamics of a quantum particle interacting with a bath made up of a huge number of harmonic oscillators obeying the Bose-Einstein statistics [6,[49][50][51].In this analogy, the impurity plays the role of the Brownian particle and the Bogoliubov excitations of the BEC are the bath-oscillators.Here, by extending this view, we study the creation of entanglement between two different impurities in a Bose gas, as a consequence of the coupling induced by the presence of the Bogoliubov modes, which play the role of the bath.
Note that the QBM model for the motion of the two kinds of impurities in a BEC is a continuous-variable description i.e. it is expressed in terms of position and momentum operators.Thus, entanglement measures based on the density matrices are not conveniently calculable because the density matrix in this case is infinite-dimensional.We therefore use the logarithmic negativity [52] as a more fitting choice in this context.This measure is expressed in terms of the covariance matrix, namely a matrix, whose elements are all of the position and momentum related correlation functions.
To compute these correlation functions, we solve the Heisenberg equations of the system, which can be reduced to a quantum stochastic Langevin equation for each particle.This set of two coupled equations are non-local in time, namely the dynamics of both impurities in a BEC carry certain amount of memory.In this context, such a feature is often related to the super-Ohmic character of the spectral density, constituting the main quantity that embodies the properties of the bath.The presence of memory effects (non-Markovianity) can also be shown to lead to the appearance of entanglement [53][54][55][56][57].The role of the memory effects in the works above is to preserve entanglement in the long-time regime [58] and in the high temperature regime [13].Nevertheless, in these works the spectral density was assumed to be Ohmic, such that the non-Markovianity is purely attributed to the influence of one particle on the other.Indeed, the disturbances caused by particles to one another is mediated through the common bath, which take a finite time to propagate through the medium, making the evolution of each particle history-dependent.This results on a decay of entanglement in several stages [13,17] or to a limiting distance for bath induced two-mode entanglement [54].In [59], the scenario of an additional source of non-Markovianity, emerging from a non-Ohmic spectral density was considered.The non-ohmic spectral density resulted in more robust entanglement among the two impurities.
In this work we study entanglement as a function of the physical quantities of the system, such as temperature, impurity-gas coupling, gas interatomic interaction and density.These parameters may be tuned in experiments allowing to control the amount of entanglement between the impurities.We distinguish the situation in which the impurity is trapped in a harmonic potential and that where it is free of any trap.In the trapped case, we also study squeezing which is a resource for quantum sensing.
The manuscript is organized as follows.In Sec. 2 we present the Hamiltonian of the system, showing that it can be reduced to that of two quantum Brownian particles interacting with a common bath of Bogoliubov modes.We also write the quantum Langevin equations, find the expression for the spectral density showing that it presents a super-ohmic form, and solve the equations in order to evaluate the position and momentum variances.Finally, we review and discuss the logarithmic negativity as an entanglement quantifier and a criterion which we use to detect two-mode squeezing.In Sec. 3 we study the out-of-equilibrium dynamics of untrapped impurities, while in Sec. 4 we study entanglement and squeezing of the two impurities as a function of the system parameters for the trapped case.In Sec. 5 we offer the conclusions and outlook.We discuss details on the derivations of the spectral density and susceptibility in appendices A and B. In App.C, we comment on the difficulty of finding an analytic solution for the trapped case even when the centers for the particles potentials coincide, and in App.D we study, for the trapped case, the effective equilibrium Hamiltonian of the system reached at long-times.

Hamiltonian
We consider two kinds of distinguishable impurities of mass m 1 and m 2 , immersed in a bath of interacting bosons of mass m B enclosed in a box of volume V .This system is described by the Hamiltonian H = H I + H IB + H IB , with with a † and a being the bath creation and annihilation operators respectively, x j and p j the position and momentum of particle j = 1, 2. Here, we consider the bosons to be in a homogeneous medium, i.e.V (x) = const, and the external potential experienced by the impurities is that is a 3D parabolic potential centered at d j = (d j,x , d j,y , d j,z ) and with frequency Ω j = (Ω j,x , Ω j,y , Ω j,z ).We will consider both the trapped Ω j > 0 and untrapped cases Ω j = 0.The densities of the impurities in the momentum domain ρ (j) The Fourier transforms of the interactions among the j th impurity and the bath and among bath particles themselves are, respectively The coupling constants in Eq. ( 5) are where m , is the reduced mass, and a (j) IB and a B are the scattering lengths between the impurities and the bath particles and between the bath particles themselves, respectively.
By performing a Bogoliubov transformation, with the terms related purely with the bath particles transform to the following non-interacting term where we neglected some constant terms and where is the coherence length for the BEC and is the speed of sound in the BEC.Finally, n 0 is the density of the bath particles in the ground state, which turns out to be constant as a result of the homogeneity of the BEC.For a macroscopically occupied condensate, i.e.N j =0 ≪ N 0 where N j is the occupation number of the j th state, we obtain the following expression for the interaction part of the Hamiltonian which after the Bogoliubov transformation, becomes Here, the couplings are After this procedure the Hamiltonian is transformed into I + H The above Hamiltonian with a non-linear interacting part is in general difficult to treat, and requires the usage of influence functional techniques [8].To simplify the situation we will resort to the so called long-wave or dipole approximation.This is expressed by the assumption such that we can approximate the exponentials by a linear function.The validity of the linear approximation was studied in [43].For the untrapped impurities, it turned into a condition of a maximum time as a function of T for which the linear approximation holds.For the trapped case, the condition was over the trapping frequency as a function of T .The results we present here fulfil these conditions.With this, the Hamiltonian reads I + H (2) with I the identity and where and

I, one obtains the following Hamiltonian
with where R jq = |d j − d q |, and g k .At this point we note that, to preserve the bare oscillator potential in Eq. ( 8) and hence ensure a positively defined Hamiltonian, it is conventional to introduce a counter term to the Hamiltonian.In this way the Hamiltonian at hand can be interpreted as a minimal coupling theory with U (1) gauge symmetry [60].This term is important because its introduction guarantees that no "runaway" solutions appear in the system, as shown in [61].Nevertheless, we are committed not to introduce any artificial terms in the Hamiltonian and maintaining the fact that we are considering a Hamiltonian that describes a physical system.We will however identify a condition under which the problem of "runaway" solutions will not appear.
For the trapped impurities case, Hamiltonian (18) shows similarities to the Hamiltonian that describes the interaction of three harmonic oscillators in a common heat bath [59].However, the two Hamiltonians differ in the following three aspects: First, our Hamiltonian describes the interaction as a coupling between the position of the particle and a modified expression of the momentum of the bath particles while [59] describes a position-position interaction.Second, the term W (d 1 , d 2 ) [x 1 + x 2 ] is absent in the Hamiltonian in [59]; Finally, our Hamiltonian lacks the counter term which is artificially introduced in [59] that results in a renormalization of the potential of the harmonic oscillator.

Heisenberg equations
From here on we treat only the one dimensional case, i.e. we assume that the BEC and the impurities are so tightly trapped in two directions as to effectively freeze the dynamics in those directions.In practice, the one dimensional coupling constant has to be treated appropriately, as discussed in [62].To study the out-of-equilibrium dynamics of the system, we first obtain the Heisenberg equations, which read as Next we solve the equations of motion for the bath particles (20c) and (20d), k e ikd j x j (s) e iω k (t−s) ds.
Replacing these in Eqs.(20a)-(20b) yields the following equations of motion for the two particles where B j (t, d j ) plays the role of the stochastic fluctuating forces, given by and the memory friction kernel λ jq (t − s) reads as with, g(jq) k .In the literature, λ jq (t), is also refered to as the dissipation kernel, or also the susceptibility.The Heaviside function guarantees causality by introducing the corresponding retardation due to the finite distance d j − d q between the centers of the two harmonic oscillators.The two are related by the Kubo formula as Hence, we understand that the reason the Heaviside step function arises is due to the fact that the forces commute for time-like separations.On the left hand side of Eq. ( 21) appears a restoring force which is originated by the fact that the impurities are similar to two harmonic oscillators.Furthermore, also non-local terms appear due to the interaction of the impurities with the environment, in particular a dissipative self-force and a history-dependent non-Markovian interaction between the two impurities.Both of these non-linear terms are a consequence of the coupling of the impurities to the bath.On the right hand side, the stochastic fluctuating force appears with Gaussian statistics, i.e. the first moment, which is assumed to be B j (t, d j ) = 0 for j = 1, 2, and the second moment of the probability distribution related to the state of the stochastic driving force B j (t, d j ) is enough to describe the state of the bath.
There is another equivalent way of writing the equations of motion for the particles, in terms of the damping kernel Γ jq (t − s), related to the susceptibility as 1 m j λ jq (t − s) := − ∂ ∂t Γ jq (t − s), which will enable us to identify a condition on the range of parameters for which our system is valid.Making use of the Leibniz integral rule, one can show that With this, one can rewrite the equations of motion for each particle position in terms of the damping kernel One can identify that Γ jq (0) x q (t) in equations of motion (26) play the role of renormalization terms of the harmonic potential.Most importantly, these terms will not be present in case one includes a counter term in the initial Hamiltonian, Eq. (18).Under the assumption of weak coupling, one expects that these terms will not affect the long-time behaviour of the system, as explained in [49].At the end of this section we obtain the necessary condition that ensures the positivity of the Hamiltonian.It is useful to write Eqs.(21) as a single vectorial equation, where Equivalently, the vectorial equation that corresponds to Eqs. (26) (i.e. in terms of the damping kernel), is where We then introduce the transformation matrix Q = O X such that Eq. ( 32) transforms into where i.e. it is diagonal.In the following sections we solve Eqs ( 27) and (32) for the cases under study.
We now identify from Equation (27) the condition under which even though the Hamiltonian lacks an ad hoc introduced renormalization term, this does not affect the long-time dynamics of the particles.This condition is that both Ω1 D , Ω2 D > 0, because this way the Hamiltonian remains positive-definite and diverging solutions are avoided.In particular, it is required that This imposes a restriction on the coupling constants range.Note that if we decouple the second particle i.e. if g = 0 for all k, then we obtain the same condition as for the one particle, Ω 2 1 > Γ 11 (0) [43].

Spectral density
Let us write the dissipation kernel in Eq. ( 23) as where we identify the spectral densities as and We note that g k is an even function of k, which implies that J antisym.

jq
(ω) = 0.Then, in App.A we show that in the continuum limit of the spectrum, Eq. ( 38) takes the following integral form for a system with 1D environment with where with and represents the relaxation time, with m = m 1 m 2 m 1 +m 2 .For the single impurity case, R 12 = 0 and η 2 = 0, we obtain a cubic dependence of the spectral density on ω, in accordance with the results obtained in [43].It is worth noting here that the validity of the Fröhlich type Hamiltonian we have here imposes a restriction on η j [38,63], namely Therefore, we restrict ourselves within this limit, which for typical values of the related parameters, such as for example g B = 2.36 × 10 −37 J • m and n 0 = 7 (µm) −1 which are values that were experimentally considered in [34], becomes η (cr) ≈ 7.This condition is satisfied for all the values of η j , g B and n 0 considered here.Finally, k ω in Eq. ( 42) is the inverse of the Bogoliubov spectrum and the susceptibility is where c is the speed of sound, ξ the coherence length, and One identifies two opposite limits for χ 1D (ω), i.e. ω ≪ Λ and ω ≫ Λ. Hence Λ appears naturally as the characteristic cutoff frequency to distinguish between low and high frequencies.
Note that the spectral density exhibits a super-ohmic behaviour given by the third power of the bath frequency in the continuous limit.Same behaviour was found in [64][65][66] for analogous problems, and it is attributed to the linear part of the Bogoliubov spectrum.With such a spectral density, it can be shown that certain quantities, e.g. the momentum dispersion would be divergent, unless the high frequencies are somehow removed from the spectrum.This can be achieved by introducing an ultraviolet cutoff given by Λ such that only the part of the spectrum where ω ≪ Λ remains.Upon doing so, χ 1D (ω) → 1, and the Bogoliubov spectrum takes the linear form ω = c |k| such that the spectral density becomes We consider two different possible analytical forms of the cut-off, the sharp one provided by an Heaviside function, and the exponential cutoff In [43], it is shown that the physics of the system in the long-time limit, i.e. associated to frequency regime ω ≪ Λ, does not depend on the existence, nor the form, of the cutoff.In this paper we will use both types of cutoffs depending on the problem at hand, namely we will use the exponential cutoff whenever we study the problem for distances between the traps of each kind of impurity different than 0 and the sharp cutoff otherwise.Finally, note that in comparison to the equations of motion obtained in [59], in our case the couplings of the two particles can be different, adding an extra parameter.

Solution of Heisenberg equations and covariance matrix
To evaluate the covariances one needs to solve the above equations of motion.In particular, we solve Eq. ( 32) by first obtaining the solution for the homogeneous equation, and then adding to that the particular solution [43,49], where and with L z [•] denoting the Laplace transform.Notice that the second function is often referred to as the susceptibility, in analogy with the harmonic oscillator.It basically carries the same information as ∂ ∂t Γ (t − s), which is what we refer to as susceptibility in this paper.We now need to evaluate the covariance matrix with P (t) = (p 1 (t) , p 2 (t)) T the momentum vector and Hence, the 4×4 matrix ( 55) is constructed from the vector Furthermore, if we assume the initial state of the system to be of the Feynman-Vernon type, i.e. ρ (0) = ρ 12 (0) ⊗ ρ B , then quantities like X (0) B T (t, d 1 , d 2 ) will vanish.We note that the appearance of the extra term W (d 1 , d 2 ) in the dynamics, absent in [59], indeed does not affect the evaluation of the covariance matrix since we are only interested in averages with respect to the state of the bath.Proceeding in a similar manner as in [43] for a thermal equilibrium bath at temperature T we conclude that the equal time correlation function of position reads as where Similarly, one can obtain the position-momentum and momentum-momentum blocks of the covariance matrix, and Notice that for g 2 k = 0 ∀k , i.e. by removing the second particle from the system, one reproduces the expressions provided in [43] for the single particle case, by considering the Laplace transform of the damping kernel L z Γ (t) .Importantly, in the presence of the second particle, the expressions for ( 57), ( 59) and ( 60) can still be obtained with the method described in [43], but only when d j − d q = 0, i.e. when the centers of the particle potentials are in the same position.Indeed, for d j − d q = 0, the Laplace transform of the damping kernel can be computed as However, if d j −d q = 0, L z Γ (t) cannot be obtained as the integral does not converge.In such case, we use the method presented in [59], where the Fourier transform of the susceptibility is evaluated instead, through the usage of the Hilbert transform to solve Eq. ( 27).In the following, we find that we circumvent the aforementioned problem by using the Fourier method.
By doing so, we are avoiding the integral on the imaginary axis, since the following relation is known between the Fourier and Laplace transforms for even (+) and odd (-) functions respectively, for the functions J(ω), i.e. the spectral density, and λ(ω), i.e. the susceptibility, where the definition of the Fourier transform is [67] The correlation functions take the following form where and Im[.] is the imaginary part.Here, From [59] we have where λ (ω) is the Fourier transform of the susceptibility λ (t).We remind that the bath is assumed to be at thermal equilibrium.Equation ( 64) was proven in [68].The other autocorrelation functions can be obtained in a similar way once we relate the momentum to the time derivative of the position as and To proceed further, one needs to evaluate F ω λ (t) that appears in Eq. ( 66).We already know the imaginary part of the susceptibility from Eq. (67).The real part of the susceptibility can be obtained by making use of the Kramers-Kronig relations, which mathematically means that one has to take the Hilbert transform H of the imaginary part of the susceptibility where Re is the real part and P denotes the Cauchy principal value.
For the particular form of the spectral density in Eq. ( 51), the susceptibility takes the following form: with where dt denotes the incomplete gamma function.Finally, we remark here that for the case R 12 = 0, one could proceed by using the Laplace transform as initially intended.With the work in [57] in mind, one could think that an analytic expression could be obtained, but we explain Appendix C, that this method is not applicable for a super-ohmic spectral density, and hence numerical integration should be applied instead, as is the case for the method using the Fourier transform.

Entanglement measure
We will address the existence and dependence of entanglement between the two impurities on the parameters of the model, both for trapped and untrapped impurities.For continuousvariable systems, an entanglement measure based on the density matrix is not conveniently calculable because the density matrix in this case is infinite-dimensional.For Gaussian bipartite states however, there are a number of ways to circumvent this problem [69].In general, a state that is not separable is considered entangled.The well known Peres-Horodecki separability criterion [70,71], which poses a necessary condition for separability, was shown to be easily formulated for a Gaussian quantum bipartite state through the usage of the symplectic eigenvalues of the covariance matrix of the bipartite system.However, this criterion alone does not allow for a quantification of the entanglement in the system, because it does not offer a way to quantify entanglement that is a monotonic function of the aforementioned symplectic eigenvalues [72].For pure states a unique quantification measure exists, which is the entropy of entanglement [73].This is defined as the von Neumann entropy of the reduced states of the bipartite system.For mixed states however this is not the case.In this case, there are a number of possible measures of entanglement, such as the entanglement of formation [74], the distillable entanglement [75] and the logarithmic negativity [52].The first two are notoriously hard to calculate in general.For this reason we resort to the usage of logarithmic negativity as the most convenient measure to quantify entanglement.The logarithmic negativity is defined as where ρ 12 is the density matrix of the two impurities and ν − is the smallest symplectic eigenvalue of the partial transpose covariance matrix C T 2 (0), where the partial transpose is taken with respect to the basis of only one of the impurities.Here, we briefly present the method to obtain the symplectic eigenvalues of C T 2 (0).To do so, it will be more convenient to reconstruct the covariance matrix in Eq. ( 55) using the rearranged vector Ŷ (t) = (x 1 (t) , p 1 (t) , x 2 (t) , p 2 (t)) such that the matrix becomes Thus, the diagonal matrices C 11 (0) , C 22 (0) represent the covariance matrices of the first and second particle respectively, and C 21 (0) = C T 12 (0) represent the correlations between the two particles.It can be shown that, for such Gaussian quantum bipartite systems, partial transposition corresponds to time reversal [76].Therefore the partial transpose of C (0) with respect to the second particle degrees of freedom is obtained by assigning a minus sign to all the entries of the matrix where p 2 (t) appears.The symplectic spectrum then is obtained as the standard eigenspectrum of iQC T 2 (0) where Q is the symplectic form and |.| represents taking the absolute value.The 4×4 symplectic matrix is defined as which in our case it reads as The final result is a diagonal matrix C T 2 diag.(0) with diagonal entries diag (ν 1 , ν 1 , ν 2 , ν 2 ) from which the smallest one is selected and introduced in Eq. ( 73) to obtain a quantification of entanglement.At the Hilbert space level, this symplectic diagonalization transforms the state into a tensor product of independent harmonic oscillators [52], each of which is in a thermal state, the temperature of which is a function of ν j .We recall that the Peres-Horodecki criterion states that if a density matrix ρ 12 of a bipartite system is separable, then ρ T 2 12 ≥ 0. The logarithmic negativity quantifies how much this condition is not satisfied [69].
In [77,78] it was shown that logarithmic negativity quantifies the greatest amount of EPR correlations which can be created in a Gaussian state by means of local operations, and in [79] it was shown that logarithmic negativity provides an upper bound to distillable entanglement.The form of the symplectic eigenvalues can be explicitly given for the bipartite Gaussian system above using symplectic invariants constructed from determinants of the covariance matrix as [53,80] where In this scenario, the uncertainty relation can also be conveniently expressed in terms of the symplectic invariants constructed from determinants of the covariance matrix.In particular, it becomes equivalent to the following three conditions This can also be proven to be equivalent to the condition that the lowest eigenvalue of the symplectic matrix of C (0), is larger than 1 2 , i.e.
We notice here that the logarithmic negativity can take arbitrarily large values as ν − can in principle go to 0, with the uncertainty principle still being satisfied.To attain maximal, finite entanglement, one has to fix the values of both local and global purities of the state of the twoimpurity system [77,78].Note also that one can construct an estimate of entanglement, the average logarithmic negativity, that is a function of these purities, which are easier to measure in an experiment [77,78].However, in this work we focus on the study of the logarithmic negativity itself, because for the amount of entanglement that we find in our numerical results, average logarithmic negativity is not a good estimate, i.e. the error as defined in [77,78] is large.

Squeezing
Once the covariance matrix is obtained, squeezing in the long time limit state of the system can also be rigorously studied.To this end one can use a set of criteria identified in [80].As is noted in this work, squeezing in phase space in a system is a consequence of the appearance of non-compact terms in the Hamiltonian, i.e. terms that do not preserve the particle number, of the form This is indeed the form of the interacting part of the Hamiltonian in Eq. ( 18) once we write it in terms of creation and annihilation operators of the harmonic oscillators and observe that in the Hamiltonian, terms of the form appear.The criterion derived in [80] states that if where ν min is the smallest normal eigenvalue of the covariance matrix (not to be mistaken with the symplectic eigenvalue), then the state is said to be squeezed.
3 Out-of-equilibrium dynamics and entanglement of the untrapped impurities In this section, we study the case of untrapped impurities, Ω 1 = Ω 2 = 0. Before presenting any results, we note that the following checks are made in order to guarantee that the assumptions presented in the previous section are valid.For all the parameters for which we present results here we check: i) that the constraints for the validity of the linear approximation described in [43] hold; ii) that the condition for the positivity of the Hamiltonian holds; iii) that the interactions are not so strong as to make invalid the initial Hamiltonian, as discussed in [38,63] (see also discussion in [43]); iv) that the Heisenberg uncertainty principle condition holds; v) that the high temperature limit obeys equipartition.The fact that the latter does not hold is an indication of problems with the numerical integration.The same conditions hold for the trapped case scenario.
In the untrapped impurities case, the time dependent expressions for the Mean Square Displacement (MSD), defined in Eq.90, the average energy and the entanglement can be obtained analytically.We emphasize that some of these quantities, such as the MSD, can be measured in the lab for ultracold gases [34].Our first aim, is to solve the equation of motion (52).To obtain analytic expressions for G 1 (t), G 2 (t), we first consider the particular form of their Laplace transforms, Eqs. ( 53) and (54).They now read as Equation (80a) is independent of the specific form of the damping kernel, and can be easily inverted to get Equation (80b) depends on the form of the damping kernel.Thus, to get G 2 (t), one needs to compute the Laplace transform of the damping kernel, Eq. (61).In the long time limit, i.e.
when Re[z] ≪ Λ, the Laplace transform of the damping kernel reads and hence where This can be inverted as which we rewrite as with Entanglement is observed at the long but not infinite time limit, for the case of untrapped impurities of potassium K in a bath of particles of Rubidium Rb, at the low temperature limit.
The initial variances of position and velocity for the two particles are assumed to be 0, i.e.
x 2 j (0) = ẋ2 j (0) = 0 for j ∈ {1, 2}, however the qualitative behavior was the same for other initial conditions as well.Entanglement, decreases to zero as time passes.It is studied for a number of different coupling constants η 2 of the second impurity.The rest of the parameters are Ω = 2π • 500Hz, η 1 = 1, g B = 3 • 10 −37 J • m and n 0 = 7(µm) −1 .It is observed that increasing η 2 decreases both the value of the entanglement and the time at which it reaches zero.
Finally, the solutions of the equations of motion for the two impurities, written in the form of Eqs.(52), read as where Now, we can evaluate, first, the MSD for each one of the particles.The MSD is defined as For the sake of simplicity, and to study the dynamical evolution of the MSD of the impurities purely due to their interaction with the bath, we assume that the initial states of the impurities and the bath are uncorrelated, ρ (0) = ρ I (0)⊗ρ B , such that averages of the form ẋj (0) B q (s) , that would otherwise appear in the expression, vanish.For simplicity, we also assume that there are no initial correlations between the two impurities such that the terms x j (0) x q (0) , ẋj (0) x q (0) and ẋj (0) ẋq (0) for j = q also vanish.Furthermore, to evaluate the MSD one needs to evaluate where is the noise kernel.To prove this, we used the fact that for a bath mode at thermal equilibrium at a temperature T , Then the expression for the MSD of one of the particles, in the long time limit, takes the form In the regime of low temperatures, where coth Therefore, we find that the particles motion is superdiffusive.We note here that the same result was found for the single particle case [43], where this effect was attributed to the memory effects present in the system.In this context the result in Eq. ( 90) represents a witness of memory effects on a measurable quantity.
In [81,82] the presence of memory effects is associated to backflow of energy.To examine whether such backflow of energy appears in our system as well, we derive an expression for the average energy of the system as a function of time.To do so, we need an expression for the time evolution of the momentum, which reads as Thus, in the low temperature limit, the average energy as a function of time, reads as The oscillatory behaviour of the energy suggests that, in addition to the traditional dissipation process where the impurity loses energy, also the environment provides energy to the impurity, i.e. we detect a backflow of energy from the environment to the impurity.We note that, in the two particles case, the diffusion coefficient in Eq. ( 90) is different for each particle [see expression for α jq , Eq. ( 85)].Thus, it depends on the interactions of each kind of particle with the BEC and the mass of each particle, together with the density and coupling constant of the BEC [see expression for the cutoff frequency, Λ, Eq. ( 48)].
As explained in Sec.2.4.1, we will use the logarithmic negativity to study entanglement and hence the covariance matrix of the two impurities is needed.To this end, we find for the low temperature case, coth ( ω/2k B T ) ≈ 1 and in the long-time limit, Re[z] ≪ Λ δ ky x k (0)x y (0) + δ ky α jk α qy ẋk (0) ẋy (0) t 2 (93a) δ ky m q α jk α qy ẋk (0) ẋy (0 δ ky m j m q α jk α qy ẋk (0) ẋy (0) (93c) where in Eq. (93a t dt is the cosine integral function.In Fig. 1 we depict the entanglement as a function of time for the low temperature scenario and for different coupling constants η 2 .We find entanglement in the long-time limit, which vanishes linearly.The maximum time at which entanglement reaches zero is increased by decreasing the interactions of the second kind of impurities.For a given time, entanglement increases with decreasing η 2 .Also, for a given time we found that increasing η 1 while keeping the ratio η 1 /η 2 constant, increases entanglement.We also studied the dependence of entanglement on density (not shown), finding that for large enough densities, the higher the density, the less entanglement was found and the faster it disappears.However, for low densities, increasing density increases entanglement.Since it is the presence of the bath that entangles the particles, small density is necessary to induce entanglement.This is in accordance with the results for the trapped impurities, presented in next section.Quantitatively, entanglement of an order of magnitude larger was found for a density of an order of magnitude smaller than the one presented in Fig. 1.We note that, for each curve, there is a minimum time at which equipartition is fulfilled according to Eqs. (77).The form of ν − can be analytically obtained from the expressions (93).We do not present it here explicitly to avoid long expressions.
Finally, the MSD for the relative motion r = x 1 − x 2 and center of mass R = (x 1 + x 2 )/2 coordinates can also be obtained.We find that they equally perform superdiffusive motions.This means that both the variance of the distance between the two particles and the center of mass increases ballisticaly on average.For the particular case where the impurities are of the same mass and interact with the same strength with the BEC, the relative distance variance is constant in time.This is showing that it decouples from the bath.The center of mass position variance instead still grows ballisticaly with time.Instead, we find that for this case still one can find entanglement.This indicates that the vanishing of the entanglement at large times is not explained solely but the increase of the relative distance variance.

Squeezing and Entanglement for Trapped impurities
We conjecture to find squeezing and entanglement between the two particles in the regime where quantum effects play an important role.Thus, we consider the low temperature regime, namely the regime where k B T ≪ Ω j with j = 1, 2. The parameters that we use are such that the condition (37) that guarantees that no runaway solutions are encountered, is satisfied.At the same time the coupling constants used are relatively strong such that the non-Markovian effects are manifested.We find that for entanglement, one has to consider coupling constants in the range that satisfies: g IB /Ω j ∈ [0.01, 1], as below this range the bath effect is not enough to create entanglement, while above this the effect of the bath destroys entanglement.In the trapped case we make use of the covariance matrix whose elements are constructed by numerical integration from Eq. ( 64), ( 68) and (69).In general, and unless stated otherwise, we make the assumption that the centers of the harmonic traps are equal, d 1 = d 2 , as entanglement is maximized in such case.

Squeezing
In this section we study squeezing as a function of the parameters of the system and the bath.To detect squeezing we make use of the condition (79).However, note that the value of ν min is not a measure of the level of squeezing in the system, but a criterion that squeezing occurs.In the numerical computations presented in Fig. 2, we take Ω 1 = 600πHz, Ω 2 = 450πHz, n 0 = 90(µm) −1 , η 1 = 0.325 and R 12 /a HO = 0, where distance was measured in units of a fixed harmonic oscillator length for both impurities equal to α HO = /(mΩ), with m = m 1 = m 2 and Ω = πkHz a typical frequency of the same order of magnitude as the frequencies considered throughout all of our studies.Without loss of generality, we assumed d 1 ≥ d 2 such that R 12 ≥ 0. In our studies, we varied the temperature T , η 2 and g B .The qualitative behaviour for a varying n 0 was found to analogous to that of g B , so the results with respect to this variable are not shown.The parameters used are within current experimental feasibility [34].In Fig. 2(a) we studied squeezing as a function of temperature The T crit at which the entanglement goes to zero increases with η 2 and seems to saturate at large η 2 .(c) the entanglement at fixed T = 4.35nK increases for a range of η 2 and then decreases to zero.In this case we considered n 0 = 350(µm) −1 and g B = 9.75 • 10 −39 J • m.In all of the graphs, we are considering impurities of potassium K in a bath of particles of Rubidium Rb.The parameters used in these plots are: for a number of different η 2 and, in panel (b), squeezing as a function of the temperature for various g B .First, we find squeezing at low temperatures, in the nK regime, that vanishes at higher temperatures.Second, we observe that the temperature at which squeezing vanishes increases with the coupling constant η 2 or g B .Furthermore, as we will show in next section, squeezing appears in the range of temperatures where entanglement appears as well, as can be seen in Fig. 3, but the range is slightly larger than that for entanglement.The existence of a relation of squeezing and entanglement is in agreement for example with the work in [83], where an ohmic spectral density was considered, such that analytic results could be obtained, and they find in the long-time limit that the logarithmic negativity was given by E LN (ρ 12 ) = 2r with r being the two-mode squeezed state squeezing parameter.
We also studied the position and momentum variances, observing that indeed in the large temperature limit they approach the equipartition theorem.This means that, in this limit, the system is formally analogous to two independent harmonic oscillators as expected.We used this as a test to verify the validity of our numerical results.In appendix D we study the equilibrium Hamiltonian in detail.This allows us to find a prediction for the large temperature limit of, e.g., x 1 x 2 or p 1 p 2 .We found that these correlation functions do not vanish at large T , not implying the presence of quantum correlations but only classical correlations in this limit.This is in agreement with the fact that one only finds entanglement for very low T .
In addition, we calculated the uncertainty δ x 1 when compared to δ p 1 , restricting only to the regime where Heisenberg uncertainty principle was fulfilled.This amounts to studying we study entanglement as a function of T for various distances between the trap of each kind of impurity.As expected the entanglement is reduced as the distance is increased.This is shown in figure (b).For each distance, at certain T the entanglement vanishes.The dependence of this critical temperature as a function of distance is shown in (c).In all of the graphs, we are considering impurities of potassium K in a bath of particles of Rubidium Rb.The parameters used in these plots are: squeezing for the partially traced state of the system, tracing out the other kind of impurities.This way we were able study how the introduction of a second kind of impurities modified the squeezing found when defined as that of only one impurity (as is done in [43]).We found that the squeezing observed for one particle reduces as η 2 increases.

Thermal entanglement induced by isotropic substrates
Here we study the appearance of thermal entanglement, i.e. assuming that the entangled resource, the bath, connecting the two impurities is in a canonical Gibbs ensemble density matrix at certain T .In general the following parameters were used: In certain cases we consider other values of n 0 and g B to study the appearance of the phenomenon of the bath causing a decrease of the entanglement.Also, some general comments about the results that will be presented below, are the following.First, it was observed that parameter regimes existed in which the uncertainty principle, translated into the condition Eq. ( 77), was not satisfied.We only present results when this condition is fulfilled.In the results presented here, entanglement is normalized based on the instance of maximum entanglement found in the system, which was obtained for the following parameters: For all figures where we show the dependence of entanglement on temperature, we emphasize that below a certain temperature, the uncertainty principle was not satisfied.This minimum temperature depends on the other parameters of the system.For example, the larger the distance considered was, the lower temperatures that one can reach under the requirement that the uncertainty principle holds.We also note that for low enough temperatures, we find numerically a saturation of entanglement in all cases.Finally, the general behaviour of entan- In (c) we illustrate that the entanglement at fixed T increases for a range of n 0 and then decreases to zero, as it was the case with η 2 in Fig. 1.We show this dependence for various values of the coupling constant among the bosons.As shown, increasing g B increases the maximum value of the entanglement reached, but also entanglement vanishes at a smaller value of the density.In all of the graphs, we are considering impurities of potassium K in a bath of particles of Rubidium Rb.The parameters are: glement with temperature is to decrease, as expected, and beyond a certain temperature it vanishes.We term this as critical temperature, T crit In Fig. 3 (a), the temperature dependence of entanglement was studied for a number of coupling constants η 2 , and it was observed that increasing η 2 implied an increase in entanglement, as well as an increase in T crit .Information about this temperature is particularly important experimentally, and for this reason we studied the dependence of T crit on η 2 in Fig. 3  (b).We see that the increase of T crit with η 2 is decreasing for larger η 2 .The dependence of entanglement was studied also as a function of η 2 .In this case we considered n 0 = 350(µm) −1 and g B = 9.75 • 10 −39 J • m which allowed us to see the diminishing effect of the bath, meaning that entanglement reached a peak value for some η 2 and was later then decreases with increasing values of η 2 .
In Fig. 4 (a), we study entanglement as a function of the temperature for various distances between the two impurities.As expected, entanglement decreases with increasing distance.Furthermore, in Fig. 4 (c) the T crit also decreases with distance and it acquires a maximum for distance equal to 0. In Fig. 4 (b) we see that entanglement drops to 0 beyond a certain distance which is R 12 ∼ 0.35a HO , which for the parameters that we have chosen results in 0.2µm.The distance at which entanglement drops to 0 depends on the other parameters of the system as well, e.g. it increases with decreasing temperature, but remains at the same order of magnitude.
In Fig. 5 (a), we show the dependence of entanglement with T for various densities.The dependence of T crit on the density in Fig. 5 (b) again shows that it increases with n 0 for small densities.However, in Fig. 5 (c) we show that while for small densities entanglement grows with the density, for larger values of the density it starts to decrease toward zero.We plot this figure for various values of g B showing that the larger the bosons interactions the smaller the value of n 0 at which the entanglement starts to decrease to zero.Similar studies with similar results were undertaken for the dependence of the system on g B and are presented in Fig. 6.There, for g B = 4.2•10 −38 J •m we see the qualitative behaviour of entanglement with a varying temperature.These results, together with the fact that we do not assume a particular form for the state of the two impurities initially, are clear indications that the induced entanglement between the two impurities is an effect of their interaction with the common bath.
In Fig. 7 (a) we present the results of the dependence of entanglement with T for various ratios Ω 2 /Ω 1 .In Fig. 7 (b) we observe is that T crit is not monotonically dependent on the ratio of trapping frequencies.This implies that there is a regime where, even though one keeps the temperature constant at a value where there is no entanglement, by increasing the trapping potential, such that the second impurity is more confined, entanglement appears for this given temperature.In Fig. 7 (c) we see that at resonance (Ω 2 /Ω 1 = 1), entanglement achieves its maximum.This was also observed in Ref. [57].The reason is that, since the values of the entries of the matrix Γ (t) we use are so small, what guarantees that the equations of motion for the two harmonic oscillators cannot be decoupled into two equations of motion for the center of mass and relative distance degrees of freedom, is the fact that the two trapping frequencies are different.To clarify this point, and draw the parallel with the results in [57], we note that at relatively high temperatures, in which however entanglement can still be observed, one could assume that the system is described by the effective Hamiltonian Eq. (105) in Appendix D, where the effect of the bath degrees of freedom is represented by an effective coupling between the two particles.In the limit K/Ω j → 0 with j ∈ 1, 2, for  In (a) we study entanglement as a function of T for various ratios between the trap frequencies of each kind of impurity.For each Ω 2 /Ω 1 , at certain T crit , the entanglement vanishes.The dependence of T crit as function of Ω 2 /Ω 1 is shown in (b).(c) illustrates that the entanglement at fixed T has a maximum value for certain value of Ω 2 /Ω 1 .The parameters n 0 = 500(µm) −1 and g B = 9.75 • 10 −39 J • m were used in this case.In all of the graphs, we are considering impurities of potassium K in a bath of particles of Rubidium Rb.The parameters used in these plots are: Ω 1 = 600πHz, η 1 = 0.325, η 2 = 0.295, n 0 = 90(µm) −1 , T = 4.35nK, g B = 2.36 • 10 −38 J • m and R 12 /a HO = 0. Ω 1 = Ω 2 , the aforementioned decoupling takes place, and the state of the system goes to a non-symmetric two-mode squeezed thermal state with infinite squeezing, i.e. the ideal EPR (maximally entangled) state.We finally note that, as discussed in Appendix D, it is possible to find a approximate prediction for the critical T , given by Eq. (107).We found that this prediction is in the same order of magnitude (nK) for all numerical results presented.

Conclusions
In this paper, we studied the emergent entanglement between two distinguishable polarons due to their common coupling to a BEC bath.To this end, we formulated the problem of two different kinds of impurities immersed in a BEC as a quantum Brownian motion model.The BEC is assumed to be confined in one dimension and homogeneous.The impurities do not interact among themselves, but only with the BEC.By means of a Bogoliubov transformation we diagonalize the part describing the BEC in the Hamiltonian.This brings the total Hamiltonian into a form in which one can identify the BEC part as a collection of oscillators with different frequencies, thus resembling a bath Hamiltonian.Also, we identify the impurities part of the Hamiltonian as the system Hamiltonian, in the usual terminology from open quantum systems.Finally, under the same physical constraints discussed in [43] for the Bose polaron problem, we linearise the interaction Hamiltonian, which brings the Hamiltonian into a conventional Caldeira-Leggett-like Hamiltonian, i.e., a quantum Brownian motion model.
We henceforth solve the associated coupled quantum Langevin equations system of motion, which encode the bath as a damping and a noise kernel.The damping kernel includes a non-diagonal term, often called hydrodynamic term in the context of Brownian particles, as it encodes the effect of the particles on one another.We find that the spectral density characterizing the bath is superohmic in 1D.We emphasize that in our work the properties of the bath, and particularly the spectral density, are not arbitrarily assumed but derived from physical considerations of the BEC.We solve these equations both for the case of untrapped and trapped impurities.We do not add artificial terms to the Hamiltonian as to make it nonnegative.Instead we find a condition on the parameters of the system, for which the energy spectrum of the Hamiltonian is positive.
For the untrapped case we were able to solve the equations of motion analytically.Hence, we studied the MSD and the diffusive properties of the impurities, which are found to perform a super diffusive motion.In addition, we studied the momentum variance of the impurities, observing an energy back-flow from the bath, attributed to its non-Markovian nature.Moreover, we obtained the covariance matrix explicitly.By using the covariance matrix, we quantified entanglement between the two types of impurities using the logarithmic negativity.This is found to decrease linearly as a function of time.What is more, the relative distance and center of mass coordinates were considered for the case of identical impurities.The former becomes decoupled from the bath, such that it no more performs a superdiffusive motion, hence the variance of the distance between the impurities stays on average constant.Yet, we can detect a linearly decreasing with time entanglement between the two types of impurities, so we conjecture that the decrease of entanglement is attributed to their interaction with the bath rather than them running away from each other.This conjecture is further enhanced from our studies of the entanglement dependence on the rest of the parameters of the system, i.e. the density of the bosons n 0 in the BEC and their interaction strength g B .In particular, we found that for any fixed finite time, and fixed n 0 (g 0 ), entanglement reaches a maximum value at some optimal g 0 (n 0 ).Increasing this value beyond the optimal reduces the entanglement until it vanishes.
For the trapped case, we obtained the covariance matrix elements numerically.We saw that the coherence correlations (off-diagonal terms of the covariance matrix) were linearly increasing with temperature, unless the parameters of the impurities were the same.Nevertheless, entanglement decreases as a function of temperature, which means that these correlations are not quantum.In the case of trapped impurities, entanglement was studied in detail as a function of the rest of the parameters of the system.It was found to decrease with increasing distance between the centers of the two trapping potentials.Furthermore, entanglement was found to increase and then decrease as a function of both n 0 and g B .Moreover, entanglement is maximized at resonance of the frequencies of the two trapping potentials as was seen in previous similar studies as well [57].
Beyond entanglement, for all of these parameters the dependence of the critical temperature was also studied.In Appendix D a rough estimate of the critical temperature is made, using the effective form of the Hamiltonian in the thermalized regime.The estimate is of the same order of magnitude as the critical temperature observed in our studies.Squeezing was also examined in this case as a function of all the parameters of the system, found to behave qualitatively the same as entanglement.In particular, it is seen that entanglement always appears if there is squeezing but the converse is not always true.
In summary, we have studied the emergence of entanglement of two types of impurities embedded in the same bath, starting from a Hamiltonian justified on physical grounds.We examined analytically the case of two untrapped impurities, and gave numerical results in the scenario of harmonically trapped impurities.The dependence of entanglement, squeezing as well as the critical temperature, i.e. the temperature beyond which entanglement vanishes, were studied as functions of the physical parameters of the system.The parameters of the system used were within current experimental settings and we believe that our results can be experimentally verified.These results on squeezing and entanglement in these setups are particularly interesting as the two phenomena represent resources for quantum information processing.
while the real part is given by where H represents the Hilbert transform, Θ is the heaviside step function and P is the principal number.We first find as we can evaluate it with the property of the Hilbert transform We apply this property three times to obtain We can then also show that: where the top case corresponds to ω ′ ∈ (0, ∞) and the bottom to the complementary interval dt denotes the incomplete gamma function.After introducing (96) in the expression for the real part of the susceptibility above, we obtain that where With this, the susceptibility takes the form in (71).

C Study of an exact expression for the covariance matrix elements
In the particular case of R 12 = 0, it can be shown that the integral in Eq. ( 57) takes the following form: where h 6 (ω) is a 6th order polynomial, and g 6 (ω) is a 7th order polynomial.An integral of this form was also obtained in [57], for an ohmic spectral density.The criterion to use this formula, is that the roots of h 6 (ω) lie in the upper half plane.Finding the roots of this polynomial requires finding expressions for 12 variables, 6 real and 6 imaginary.In our case one can show that the polynomial h 6 (ω) can be written as with the coefficients A j being all real, which implies that the roots of h 6 (ω) are symmetrically located about the imaginary axis.This reduces the problem in finding just 6 variables, 3 real and 3 imaginary.Furthermore, by making use of the Vieta relations, one can show that This implies that not all roots of the polynomial can lie in the upper half plane.The reason for this can be traced back to the fact that there is no linear term in the polynomial h 6 (ω), which is an artefact of considering a super-ohmic spectral density as can be understood by comparing to the case in [57].The issue of not all roots being in the upper half plane, can be resolved in the following way.The roots of h 6 (ω) that are in the lower half plane, have mirrored roots in the upper half plane in the polynomial h 6 (−ω).If the expressions for the roots could be derived, one could simply redefine polynomials h 6 (ω) and h 6 (−ω) to h6 (ω) and h6 (−ω) such that all the roots of h6 (ω) would lie in the upper half plane and all the roots of h6 (−ω) would lie in the lower half plane.However for a 6th order polynomial of the form we have, one can not find the roots.So this could perhaps only be applied in an algorithmic way, where one first selects a set of parameters and then makes the redefinition of the polynomials once the roots are found.

D Equilibrium Hamiltonian
Here we discuss the thermalization properties of the system.We will use the fact that the two kinds of impurities are formally analogous to two oscillators.The bath dependent part of the Hamiltonian ( 18) is For the system to thermalize, there should be no memory effects induced on the oscillator due to its coupling with the thermal bath.This means that the Hamiltonian will have to be independent of the bath variables {b † k , b k } k =0 evolution, i.e. ∀k = 0 the following should be fulfilled This results in the following conditions k e −ikd j x j (t) , k e ikd j x j (t) .
Replacing these expressions, for the bath degrees of freedom operators, in the initial Hamiltonian (18), it becomes: k cos (d j − d q ) x j x q .
This can be rewritten as where k g (2) k cos (d 1 − d 2 ) and W (d 1 , d 2 , x 1 , x 2 ) is the function of the remaining constant terms and terms linear in x 1 and x 2 .At thermal equilibrium, the terms in W (d 1 , d 2 , x 1 , x 2 ) will not affect the equilibrium state.Then, one can neglect them and consider the effective Hamiltonian This is formally analogous to two effectively coupled harmonic oscillators.To decouple them, we transform to the normal modes of the system, Q 1 , Q 2 .This is achieved by an orthogonal transformation of the form  With this, the effective Hamiltonian in terms of the normal modes is Here we made the assumption of m 1 = m 2 = m.To diagonalize the Hamiltonian, the condition on the rotation that should be performed reads as For this rotation the Hamiltonian is 2m . (105) Hence the system has decoupled into two harmonic oscillators with frequencies and with the following second order moments where the average is a statistical average over the bath variables.Also we find Q j Q q = Π j Π q = 0 for j = q.In the high temperature limit, coth (ν j /2k B T ) ∼ 2T k B ν j and therefore One can now express the coherences or off-diagonal correlation functions at thermal equilibrium for the initial set of variables as: In case that ν 1 = ν 2 , which happens when Ω 1 = Ω 2 and/or g (1) k for some k, then the coherence between x 1 and x 2 does not cancel out in the high temperature limit.Note however that this does not imply that entanglement survives at the high temperature limit, as this correlation is not necessarily a quantum correlation.This behaviour of x 1 x 2 is consistently verified in the numerical studies that we undertook using the original Hamiltonian.We note here that such behaviour, i.e. of asymptotic non-vanishing coherences at the high-temperature limit were identified in [84] but for the case of each particle attached to its own environment and both environments having an ohmic spectral density.The reason for the resemblance of the two cases is that in our case, one could argue that even though the two particles are coupled to a common bath contrary to [84], the particles effectively see the bath at different temperatures when ν 1 = ν 2 , in particular each one sees the bath at a temperature T ν 2 j .Hence our system in this case also reaches a non-equilibrium stationary state.
From these considerations, an estimate of the critical temperature can also be obtained.Notice that an important difference between the Hamiltonian in Eq. ( 105) and the original linear Hamiltonian in Eq. ( 18) is that the former has a spectrum that is bounded from below.In this case, and as the Hamiltonian is also self-adjoint, one can apply the results of Ref. [85], where it is proven that the critical temperature for such a symmetric system as the one we are considering is given by where ν max = max[{ν j }] with j ∈ 1, 2, r = ν max /ν min where ν min = min[{ν j }] and σ(r) = ts(t) with 1 ≤ t ≤ r such that s(t) = s( t r ) with s(x) := 1 For the general values of the parameters given in Sec. 4, the value of the critical temperature obtained assuming this effective Hamiltonian, was of the order of nK in agreement with our findings.

Figure 2 :
Figure 2: Temperature dependence of the squeezing between the two kinds of trapped impurities.In (a) we study the temperature dependence of squeezing for different values of η 2 with g B = 2.36 • 10 −38 J • m and in (b) for different values of g B with η 2 = 0.295.In both cases we use impurities of potassium K in a bath of particles of Rubidium Rb and we set: Ω 1 = 600πHz, Ω 2 = 450πHz, n 0 = 90(µm) −1 , η 1 = 0.325 and R 12 /a HO = 0.

Figure 3 :
Figure 3: Coupling strength dependence of the entanglement between the two kinds of trapped impurities.(a) Entanglement as a function of T for various couplings.As shown, for the values of η 2 considered here, increasing the interactions of the second type of impurities, η 2 , enhances entanglement, at fixed T .For increasing T , entanglement decreases and eventually vanishes at certain T crit .(b) The T crit at which the entanglement goes to zero increases with η 2 and seems to saturate at large η 2 .(c) the entanglement at fixed T = 4.35nK increases for a range of η 2 and then decreases to zero.In this case we considered n 0 = 350(µm) −1 and g B = 9.75 • 10 −39 J • m.In all of the graphs, we are considering impurities of potassium K in a bath of particles of Rubidium Rb.The parameters used in these plots are: Ω 1 = 600πHz, Ω 2 = 450πHz, η 1 = 0.325, n 0 = 90(µm) −1 , g B = 2.36 • 10 −38 J • m and R 12 /a HO = 0

Figure 4 :
Figure 4: Distance dependence of the entanglement between the two kinds of impurities.In (a) we study entanglement as a function of T for various distances between the trap of each kind of impurity.As expected the entanglement is reduced as the distance is increased.This is shown in figure (b).For each distance, at certain T the entanglement vanishes.The dependence of this critical temperature as a function of distance is shown in (c).In all of the graphs, we are considering impurities of potassium K in a bath of particles of Rubidium Rb.The parameters used in these plots are: Ω 1 = 600πHz, Ω 2 = 450πHz, η 1 = 0.325, η 2 = 0.295, n 0 = 90(µm) −1 , T = 4.35nK and g B = 2.36 • 10 −38 J • m.

! 1 0Figure 5 :
Figure5: Dependence of the entanglement between the two kinds of impurities on the density of the bosons in the bath.In (a) we study entanglement as a function of T for various densities of bosons.As shown, in this range of densities, increasing the density of bosons enhances entanglement.For each density, at certain T crit , the entanglement vanishes.The dependence of T crit as a function of density is shown in (b).In (c) we illustrate that the entanglement at fixed T increases for a range of n 0 and then decreases to zero, as it was the case with η 2 in Fig.1.We show this dependence for various values of the coupling constant among the bosons.As shown, increasing g B increases the maximum value of the entanglement reached, but also entanglement vanishes at a smaller value of the density.In all of the graphs, we are considering impurities of potassium K in a bath of particles of Rubidium Rb.The parameters are: Ω 1 = 600πHz, Ω 2 = 450πHz, η 1 = 0.325, η 2 = 0.295, T = 4.35nK, g B = 2.36 • 10 −38 J • m and R 12 /a HO = 0.

Figure 6 :
Figure6: Dependence of the entanglement between the two kinds of impurities on the coupling constant between the bosons g B .In (a) we study entanglement as a function of T for various g B .As can be seen for the largest considered value of g B , namely g B = 4.2 • 10 −38 J • m, the qualitative behaviour of entanglement with temperature, changes.For each density, at certain T crit , the entanglement vanishes.The dependence of T crit as a function of density is shown in (b).In (c) we illustrate that the entanglement at fixed T increases for a range of g B and then decreases to zero, as it was the case with η 2 in Fig.1.We show this dependence for various values of the density of the bosons.As shown, in this range of parameters, decreasing n 0 increases the maximum value of the entanglement reached, but also entanglement vanishes at larger values of g B .In all of the graphs, we are considering impurities of potassium K in a bath of particles of Rubidium Rb.The parameters used in these plots are: Ω 1 = 600πHz, Ω 2 = 450πHz, η 1 = 0.325, η 2 = 0.295, n 0 = 90(µm) −1 , T = 4.35nK and R 12 /a HO = 0.

2 F 2 G 1 /Figure 7 :
Figure7: Trapping frequency dependence of the entanglement between the two kinds of impurities.In (a) we study entanglement as a function of T for various ratios between the trap frequencies of each kind of impurity.For each Ω 2 /Ω 1 , at certain T crit , the entanglement vanishes.The dependence of T crit as function of Ω 2 /Ω 1 is shown in (b).(c) illustrates that the entanglement at fixed T has a maximum value for certain value of Ω 2 /Ω 1 .The parameters n 0 = 500(µm) −1 and g B = 9.75 • 10 −39 J • m were used in this case.In all of the graphs, we are considering impurities of potassium K in a bath of particles of Rubidium Rb.The parameters used in these plots are: Ω 1 = 600πHz, η 1 = 0.325, η 2 = 0.295, n 0 = 90(µm) −1 , T = 4.35nK, g B = 2.36 • 10 −38 J • m and R 12 /a HO = 0.

p 1 p 2
Figure 1: Time dependence of the entanglement between the two kinds of untrapped impurities.