Heat current control in trapped Bose–Einstein Condensates

We investigate the heat transport and the control of heat current among two spatially separated trapped Bose–Einstein Condensates (BECs), each of them at a different temperature. To allow for heat transport among the two independent BECs we consider a link made of two harmonically trapped impurities, each of them interacting with one of the BECs. Since the impurities are spatially separated, we consider long-range interactions between them, namely a dipole–dipole coupling. We study this system under theoretically suitable and experimentally feasible assumptions/parameters. The dynamics of these impurities is treated within the framework of the quantum Brownian motion model, where the excitation modes of the BECs play the role of the heat bath. We address the dependence of heat current and current–current correlations on the physical parameters of the system. Interestingly, we show that heat rectification, i.e. the unidirectional flow of heat, can occur in our system, when a periodic driving on the trapping frequencies of the impurities is considered. Therefore, our system is a possible setup for the implementation of a phononic circuit. Motivated by recent developments on the usage of BECs as platforms for quantum information processing, our work offers an alternative possibility to use this versatile setting for information transfer and processing, within the context of phononics, and more generally in quantum thermodynamics.


Introduction
Control of heat transport has enormous potential applications, beyond the traditional ones in thermal insulation and efficient heat dissipation. It has been suggested as a resource for information processing, giving rise to striking technological developments. A series of smart heat current control devices, such as thermal diodes [1][2][3][4], thermal transistors [5], thermal pumps [6][7][8], thermal logic gates [9] or thermal memories [10], have been proposed in the past decade. A key underlying idea in some of these devices is that heat transport associated to phonons can realistically be used to carry and process information. The science and engineering of heat manipulation and information processing using phonons is a brand new subject, termed phononics [11][12][13]. In the emerging field of quantum thermodynamics, the issues of heat transfer and heat rectification are basic ingredients for the understanding and designing heat engines or refrigerators at nanoscales. Besides the technological interests, highly controlled platforms for heat transport can potentially enable the study of fundamental theoretical questions, which will shed light on the study of thermodynamics of nonequilibrium systems [14]. In particular, heat conduction in low-dimensional systems has attracted a growing interest because of its multifaceted fundamental importance in statistical physics, condensed matter physics, material science, etc [15][16][17][18].
In past years, advanced experimental techniques have allowed for the miniaturization of heat transport platforms down to the mesoscopic/microscopic scale. Experimentally, a nanoscale solid state thermal rectifier using deposited carbon nanotubes has been realized recently [19], and a heat transistor-heat current of electrons controlled by a voltage gate-has also been reported [20]. Furthermore, it has been shown that thermal rectification can appear in a two-level system asymmetrically coupled to phonon baths [21,22]. In general, various models for thermal rectifiers/diodes that allow heat to flow easily in one direction have been proposed [1-4, 21, 23-25]. Among such platforms, interesting examples are those based on ultracold atomic systems, both at the theoretical level [26][27][28][29] as well as at the experimental one [30]. In these works, the paradigm of the model used is that of tailored time dependent protocols performing a heat cycle (often an Otto cycle). A cycle here consist of a series of steps performed in finite time where thermodynamic parameters are controlled such that heat is transferred from one bath into another. Frequently, the energy that one needs to spend in implementing the aforementioned protocols is too large for the amount of heat transferred. In addition, in such a heat cycle, the Hamiltonian must be modified as adiabatically as possible to avoid unwanted excitations. We remark that for the goals of the present work considering adiabatic changes in the Hamiltonian will require long timescales, and therefore it can have drawbacks, such as exceeding the life time of the Bose-Einstein condensates (BECs). Even though attempts have been made to overcome this problem [29,31,32], our proposal circumvents it altogether by considering a rather different paradigm. Instead of cycles, we consider autonomous heat platforms, where the working medium is permanently coupled to different baths and under the influence of external time-dependent driving.
The primary aim of the present work is to design a method to transfer heat between two BECs-with a comprehensive analytical description. The BECs are confined in independent one-dimensional parabolic traps, and kept at a certain distance such that the two BECs do not spatially overlap. Our second aim is to create thermal devices in our platform, specifically a heat rectifier. In our platform, the working medium is constructed with two impurities, each of them immersed in one of these BECs. The two impurities interact through longrange dipolar interactions (see figure 1 for an schematic of the system considered). Heat transfer through dipolar interactions has been also proposed recently in two parallel layers of dipolar ultracold Fermi gases [33]. We first show that the Hamiltonian of the impurities in their corresponding harmonically trapped BECs can be cast as that of two bilinearly coupled Quantum Brownian particles interacting bilinearly with the Bogoliubov excitations of each BEC-which play the role of heat baths for the corresponding impurities. We analytically derive the spectral density (SD) of the system and construct the quantum Langevin equations describing the outof-equilibrium dynamics of the coupled impurities [34]. We emphasize that, following this procedure we avoid approximations involved in an alternative conventional approach based on Lindblad master equations, such as the Born-Markov approximation and the rotating wave approximation. We solve the quantum Langevin equations and find the covariance matrix of the impurities. We henceforth focus on our two main aims. Firstly, we find exactly the steady state heat current between the BECs, and study how it can be manipulated by controlling relevant parameters, such as the trapping frequencies, dissipation strengths, and the physical distance of the BECs. Secondly, we introduce a periodic driving on the trapping frequency of one impurity and show that our setup can be used as a thermal rectifier. The setup that we present here can be used for implementing other thermal devices as well.
The paper is structured as follows: in section 2, we introduce the Hamiltonian of our system, present the main assumptions involved, and rewrite the Hamiltonian in a form analogous to that of two coupled Brownian particles. In section 3 we derive the Generalized Langevin equations of motion for the two impurities and study their solution in both static and periodically driven scenarios. In section 4 we present the relevant quantities of interest and in section 5 we present our main results, both for heat transfer and rectification. Finally we summarize the results and conclude in section 6. Figure 1. Schematic of the system. The red and blue shaded profiles (and circles) represent the hot and cold Bose-Einstein condensates, respectively kept at T L and T R . The green profiles (and circles) represent the two impurities trapped in their corresponding parabolic potentials, plotted as blue lines-we omit the representation of the trap for the BECs. The spring-like blue line represents the long-range dipolar interaction among the impurities.

The model
We consider a system composed of two interacting impurities of the same mass m. Each one of these impurities is embedded in a different BEC, which we label as L or R because they are trapped in parabolic potentials of frequencies W B,L and W B,R , respectively, but one trapping potential is centered around a minima located in d L and the other around d R , with certain distance among them. Each BEC has, respectively, N L and N R interacting atoms of the same species with mass m B . The two impurities are also trapped in a parabolic potential of frequencies Ω L and Ω R located around the same minima as their corresponding BEC. We allow for the scenario of these impurities being driven by an external periodic force. The two BECs do not overlap among themselves (we will specify later the condition over the distance between the minima of the trapping potential). With this, the Hamiltonian describing this setting is IB where x and x B are the three-dimensional position operators of the impurity and the bosons respectively. We assume contact interactions among the bosons and between the impurity and the bosons, with strength given by the coupling constants ( ) g j B and ( ) g j IB respectively. Here, H int. denotes the interaction Hamiltonian between the two impurities, which we will specify later. For the rest of the paper we assume that the trapping frequencies in two directions are much larger than in the third one. Therefore, the dynamics in those directions is effectively frozen and we can study the system in one dimension. From here on we assume that the minima of the potential are located at d L =−d/2 and d R =d/2. We comment here that from the form of the Hamiltonian, one can see that since there is no direct interaction between the two BECs, any relative phase between them should not play any role in the dynamics of the system. It would be interesting to study the effects of allowing such an interaction, e.g. by introducing a coherent coupling between the two BECs or by studying heat transport in the experimental setting of homogeneous BEC in ( [30]), but this goes beyond the scope of our work.
We next recast our initial Hamiltonian (1), in such a way as to describe the motion of two interacting Quantum Brownian particles in two separate BECs. The procedure, which is the same as that presented [35] for the case of a single impurity embedded in a harmonically trapped BEC, can be summarized as follows: (i) we make the BEC assumption, i.e. that the condensate density greatly exceeds that of the above-condensate particles, which results to only quadratic terms in the Hamiltonian; (ii) we perform the Bogoliubov transformation appropriate for the case of a harmonically trapped BEC; (iii) we solve the relevant Bogoliubovde-Gennes equations, in the limit of a one-dimensional BEC that yields a Thomas-Fermi parabolic density profile, as in [36]; (iv) we finally assume that the oscillations x j of each one of the impurities are much smaller than the corresponding Thomas-Fermi radius R j , which physically implies that we study the dynamics of the impurities in the middle of their corresponding bath traps, which allows us to obtain a bilinear interaction of the position of the impurities and the positions of their corresponding baths. These steps bring the Hamiltonian in the form The chemical potential μ j of the jth bath is and the Thomas-Fermi radius R j of the jth bath is At this point, it is important to note that Hamiltonian (7), which is indeed in the form of a Hamiltonian describing two coupled quantum Brownian particles, was derived from the physical initial Hamiltonian. Therefore, this Hamiltonian lacks the renormalization term that guarantees that no 'runaway' solutions appear in the system, which is often included in conventional open quantum system approaches [37]. We are not going to introduce such a term in order to guarantee the positivity of the Hamiltonian and cast the system in the form of a minimal coupling theory with ( ) U 1 gauge symmetry. Instead, we are going to determine in the next section a condition on the allowed parameters of the system that will guarantee this.
We assume dipole-dipole interactions among the two impurities. From [38] we know that the form of dipole-dipole interaction is the following: being the vacuum permeability, while μ is the magnetic moment of the dipole. The angle θ is that formed between the axis of the two dipoles and it can determined in experiments. Furthermore, in (14), L is the distance between the two oscillators, with d L and d R being the centers of the trapping potentials of the two oscillators. We then rewrite (14) as and we assume that the distance between the two oscillators centers d is much larger than the fluctuations x L , x R and hence the fluctuations difference . Importantly, let us emphasize that this is not an additional assumption, because we assumed right from the beginning that the two baths should not overlap. This indeed means that the sum of the Thomas-Fermi radius have to be smaller than the distance between the impurities, R L +R R <d. With the additional assumption we made before, namely that the impurities oscillations are much smaller than their corresponding Thomas-Fermi radius, equation (6), then one concludes . One could also tackle the problem of interacting baths, by making use of the work in [12] where one needs to consider the surface Green functions, and in this case then the assumption should be made explicitly.
Finally, after expanding the binomial series, equation (16) is rewritten as One can show that the first two terms in the parenthesis will not contribute to the dynamics of the impurities, since they are linear in the displacement operators x L , x R and hence they will appear only as constants in the equations of motion that we will study later on, which will be obtained from the Heisenberg equations of motion. The third term then is expanded as From here on we omit the tilde in the frequencies to avoid unnecessary complications in the nomenclature. Therefore we can rewrite the interaction as which models a spring-like interaction among the two impurities with The angle θ can be experimentally controlled as was recently shown in [39]. This is possible thanks to a rotating magnetic field B rot in the x-y plane that causes the dipoles to rotate at an angle θ with respect to a static magnetic field B z along the z axis. The rotation angle is related to the magnitude of the two components by q = B B tan z rot . In practice, to achieve this in the experiment, the angle θ is controlled by using a calibration procedure that corrects for the effect of eddy currents. This allowed the authors to determine the amplitude of the ac current required to produce a given rotation angle θ. For this reason, to simplify our system, we considered the scenario where the angle is fixed. However, it is important to note here that the results that we present in this paper are valid even if the angle between the dipoles cannot be experimentally controlled, but rather an average over the angle is considered. The constants in the interacting Hamiltonian as well as the power dependence on the relative distance will be different, but in the limit we consider, i.e. when , qualitatively the results will be the same, with the difference being that the distance at which the oscillators should be kept will change [40]. This can be understood by considering the work by Keesom in [41], where the angle-averaged interaction between two dipoles was evaluated. In this case, the initial expression for the distance dependent potential would depend on r −6 , but effectively the same procedure could be followed, that would just result in a modified expression for the spring constant k. We also note that, in most ultracold dipolar gases, the dipolar interaction is present together with the short range interactions arising from low angular momentum scattering [38,42]. Usually the latter is dominant and a Feshbach resonance is needed to probe the regimes where dipolar effects are prominent. However in the setting that we will consider in this work, as we maintain the two BEC baths spatially separated, the effect of the short range interactions is negligible, such that dipole interaction is the main process through which current is transferred.
Finally, there are a number of ways to drive our system, either by driving degrees of freedom of the central system, or degrees of freedom of the environment, or their coupling. In this work we focus on the first case. There are basically two types of driving that one could consider and would maintain the quadratic form of the Hamiltonian such that an analytic solution to the resulting equations of motion can be obtained. First we could consider applying a periodically driven ramp potential on the central particles degrees of freedom only, of the form x t x t X , 1 1 and ( ) t f is some periodic function, e.g.
This type of driving was considered in [43], and can represent the force exerted on the system by a time dependent electromagnetic field. For this kind of driving however, it was shown in [7] that with the setup we assume above, neither a heat engine nor a heat pump can be constructed. Furthermore, our goal is to construct a phononic diode with our setup which means that the driving should be able to induce a unidirectional flow of heat current, and this is not the case for this type of driving. Therefore, in our work we consider the only other possible type of driving on our system that maintains our Hamiltonian in a quadratic form, i.e.
where the driving is either on the trapping frequency of the central oscillators or on their in-between coupling. It was recently shown in [44] that in such escenario one can observe the appearance of phenomenon of heat rectification, and also there is the potential to construct a heat engine as it was shown in [45], by introducing a coherently driven coupling between the two oscillators. We also assume the driving to be periodic, t V with τ being the time period, such that it can be Fourier expanded as where Ω d =1/τ being the driving frequency. This type of coupling could also be implemented by a laser.

Quantum Langevin equations
Let us now derive the equations of motion for the two impurities. First we write the Heisenberg equations of motion for the bath and impurity particles We first solve the bath particles equations of motion and we replace these in the impurities equations of motion (26)- (27), to obtain plays the role of the stochastic force and reads as Here l ( ) t j is called the susceptibility or noise kernel. In this setting it can also be identified as the self-energy contributions coming from the bath, and it reads as , contain all the relevant information about the baths. Furthermore, let us assume the Feynman-Vermon initial state assumption, i.e. that the initial conditions of the impurities and the bath oscillators are uncorrelated where r ( ) 0 is the total density state, r ( ) 0 S is the initial density state of the system and ρ B is the density state of the bath which is assumed to be thermal and hence is a Gibbs state. Then, it can be shown that the Fourier transform of ( ) B t j obeys the fluctuation dissipation relation is the Fourier transform of l ( ) t j and obeys [46] Hence one concludes that, upon determination of the SD w ( ) J j , one determines the influence of the baths on the impurities. In [35], it was shown that in the continuous frequency limit, the SD takes the following form Note that the ultraviolent cutoff, that is usually introduced to regularize the spectrum in the conventional QBM model, is now given in terms of a physical quantity, the trapping frequency of the potential well of the jth bath. Here t L j j 3 plays the role of a relaxation time. Heat transport with a superohmic SD was considered for example in the energy transport in the phenomenon of photosynthesis [47]. Note that equation (31) can be rewritten in terms of the damping kernel g ( ) t j which is related to the susceptibility by l g i d , where the Leibniz rule was used as in [49]. Moreover, in [35], it was also shown that the form of the damping kernel for such a SD reads as In the limit  t 0 this damping kernel becomes In vector form the two coupled equations in (31) read as where  is the identity matrix and , e 0 e e e 0 , 0 0 , and .
Static case. . To this end one first diagonalizes the Hamiltonian, and then requires that the normal mode frequencies are positive. We diagonalize the Hamiltonian by making the transformation Q=O·X, that brings (43) into , and the frequency matrix is diagonalized as The positivity of the Hamiltonian condition is then guaranteed by requiring that which in terms of the original frequencies reads as g g g g g g k g g which is understood to play the role of a phonon Green function, and w ( ) where the real part of the susceptibility was obtained through the Kramers-Kronig relation . In terms of this solution of the static equations of motion, the positivity condition (46), can be interpreted in a different way: it guarantees that the phonon propagator, i.e. the Green function w ( ) G , has no poles in the lower half plane of the complex plane. This implies that there are no divergencies in the integrals that will be performed later on and will involve these Green functions [50].
Driven case. Now we consider the case where a driving is applied on the central system. In particular we assume that the driving is either on the oscillators' frequencies or on their inbetween coupling, as in [44]. The analytic treatment of this case is slightly more involved than the static one. We are now dealing with a periodic differential equation, and the analysis of the stability of the long-time steady state solution of the equations of motion is not straightforward. To be able to perform the stability analysis, what one usually does is to convert the periodic differential equation to a linear one by resorting to Floquet theory. This is done by converting first all the terms in the equation of motion into periodic ones, and then study the stability of the Floquet-Fourier components of the resulting Green function. The basic assumption that enables us to employ the Floquet formalism is that even though some function ( ) f t might not be periodic, another function defined based on this one as will indeed be periodic. Furthermore, such periodic function can always be expressed in terms of its Fourier components as By performing these two transformations on the equations of motion, i.e.
d e e t t k k k t i i d that can be selfconsistently solved. By following this procedure, in [44] the authors were able to obtain expressions for these coefficients Note that we will assume that the driving strength coefficient is sufficiently small such that we can ignore terms of the order of ( ) V j 3 or higher. Furthermore, since the Fourier coefficients w ( ) A k are related to the Green function, one can interpret them as describing the fundamental processes responsible for the phonon and hence heat transport. These coefficients tell us that the driving is responsible for a sudden change of the propagation frequency ω of the phonon by an amount of kΩ d . Finally, the solution of the equations of motion in this case would read as Finally, we comment that we check that the uncertainty relation holds for both cases that we consider, static and driven. This is simply expressed by the condition where nis the minimum standard eigenvalue of~( ) C 0 (ÿ is assumed to be equal to 1 in this case). In the above expression, =( ) · ( ) C W C 0 i 0 with W the symplectic matrix and ( ) C 0 the covariance matrix, that both are defined in the appendix. Note that the covariance matrix can be expressed in terms of the Green's function and the SD, for both static and driven cases as is shown in the appendix.

Heat current control between the BECs
Here we present the thermodynamics quantities of interest, in order to evaluate the performance of our system as a heat current control platform and as a thermal diode. These quantities cannot be expressed in terms of analytically known functions, due to the non-ohmic SD that describes the baths. As such, the results in next section are obtained by numerically evaluating the integrals.

Static case
We begin with the static scenario with = ( ) H t 0 Drive , and we study the behavior of heat current and the currentcurrent correlations.
The heat current. When studying a heat engine, a key quantity is the average current J L going from the left reservoir (assuming it to be the hot reservoir) to the left oscillator [16,51] (by conservation of current J L =−J R ), or equivalently the average rate at which the left bath does work on the left particle (power). In general, there are two ways to define the heat current. The first one is derived from considerations of energy conservation on the system  , since we assume that each system interacts only with its own reservoir. Under these criteria the two definitions of heat current are equivalent. In other words, in our model, the rate at which the bath looses/gains energy is equal to the energy that the system gains/looses. The more general scenario of strong coupling and hence non-separability of the system-bath was considered for a spin-boson model in [53], while the case of interactions among the baths was studied in [12]. Therefore the heat current considered here is where the second equality is valid in steady state with , and it is obtained after solving the bath particles equations of motion. Note that in (56) the current is a scalar, which is the average current, and it is averaged under steady state condition i.e. at the long time limit which is independent of the initial state of the system. In [16,51,54], by using a direct solution of the equations of motion for a noninteracting system of bath particles, the average current is proven to be equal to , the phonon occupation number for the respective thermal reservoir and being the transmission coefficient for phonons at frequency ω. This is called the Landauer formula. Note that  w ( ) can be related to the transmittance of plane waves of frequency ω across the system as in [55]. The transmission coefficient, interestingly, depends on all the parameters of the system. That is, it depends on the parameters of the baths (apart from temperature) as well as the parameters of the coupling between the system and the baths. Equation (57) describes the situation where the central region is small in comparison with the coherent length of the waves, which is the assumption we also abide to, so that it is treated as purely elastic scattering without energy loss. The dissipation resides solely in the heat baths. This implies ballistic thermal transport, which corresponds to direct point-to-point propagation of energy, contrary to the transport in bulk and disordered structures which is referred to as diffusive transport. Indeed, within the framework of modeling thermal baths by means of quantum Langevin equations, ballistic transport has been observed for chains of quantum harmonic oscillators [56], and hence this is also our case. Furthermore, note that the Landauer formula can also capture phonon tunneling, i.e. the case when phonons off-resonance with the systems vibrations cross the 'junction', showing features of quantum tunneling [57]. Finally, it is worth commenting on an implicit assumption we made. We assumed here that a unique steady state was reached, or equivalently that there are no bound states in our system, i.e. that no modes outside the bath spectrum are generated for the combined model of system and baths. The problem is that these modes are localized near the system and any initial excitation of the mode is unable to decay [58].
It is interesting to comment on two limits of (57), namely the linear limit D - where ω refers to the bath frequencies. In the first limit, the current reduces to 6 0 0 such that once the two baths are at the same temperature there is no current flow. In the second limit, it becomes where the current is independent of the temperature of the baths T L , T R but it only depends on their difference ΔT.
We conclude with one last comment regarding our system and the role that entanglement could play on the amount of heat current transported from one BEC to the other. From [44], it is known that the static current can also be expressed in terms of the off-diagonal elements of the covariance matrix. This, might lead one to consider that entanglement, the presence of which is understood to be related to these off-diagonal elements might play a role in the amount of heat current transported. However, it was shown in [59,60] that in a system of two harmonic oscillators coupled to distinct baths, as is our case, there is no long-time entanglement present.
Current-current correlations. Next, we focus on the current-current correlations which is an easily accessible quantity from an experimental point of view. This is because these correlations are related to noise, which can be experimentally measured. They contain valuable information on the nature of the fundamental processes responsible for the heat transport. Furthermore, from the current-current correlation many other quantities can be obtained, such as the thermal conductance [16,51] and the local effective temperature of driven systems [61].
The current-current time correlations ¢ ( ) JJ t t , LL , which is sometimes referred to as current fluctuations in time or current noise, is defined as the symmetrized correlation function of the current, that is We are interested in the steady state correlations, for which an expression for this can be obtained using the nonequilibrium Green's functions as was mentioned above. Hence the correlation function of interest is a function only of the time difference, LL 0 The first two terms of this expression correspond to the equilibrium noise, while the third corresponds to the non-equilibrium noise, also referred to as shot noise. At high energies, the latter is negligible. Note that the shot noise is negative and hence contributes to diminish the noise power in comparison with having the equilibrium noise alone. The expression above equation (63) is true only under the assumption of independence of the two baths. Finally let us address the linear and the classical limits of the correlations. In the first limit, the currentcurrent correlations read as In the classical limit it becomes Note that (64), contrary to the expression for the current (60) at the same limit, does not vanish when D  T 0, which results in nonzero fluctuations of the current even in the scenario that no average current flows in the system.
We comment here on the experimental feasibility of the measurements of the two proposed quantities above, namely the heat current and the current variance. As is proven in detail in [44] in appendix C, the heat current depends on the covariance matrix element ( ) C 0 XP defined in appendix of our paper, where the nonzero elements of this matrix are á ñ x p i j with ¹ i j. Hence one needs to measure simultaneously the position of one particle and the momentum of the other. It should be experimentally feasible to make this measurements almost instantaneously. For the former type of measurement, that of position, there are already experiments in which one is able to evaluate them [65]. The idea is that, one measures the position of the particle using a time-of-flight experiment, by implementing a resonant in situ absorption imaging technique, in a system with a two species ultracold gas, in which one of the species is much more dilute-dilute enough as to consider its atoms as impurities immersed in a much bigger BEC. A much more recent technique could be used to make measurements of both the position and the momentum. In particular, a quantum gas microscope [66,67] may be an option. This technique uses optical imaging systems to collect the fluorescence light of atoms, and it has been used in the study of atoms in optical lattices, achieving much better spatial resolution [66,67], and avoiding the aforementioned problem with time-of-flight experiments. In the past such a technique has been used to study spatial entanglement between itinerant particles, by means of quantum interference of many body twins, which enables the direct measurement of quantum purity [68]. Finally, having these measurements at hand, one can evaluate the ( ) C 0

XP
. The current-current correlation then is the variance of the current, which could be obtained by the data collected. This would be the protocol to follow in order to measure the heat current. Nevertheless, we can only give an idea of a measurement protocol, leaving open for future research the question of whether the resolution when measuring the correlations in current experimental set-ups would be enough as to infer the heat current.

The dynamic case
For the driven case, the steady-state averaged heat current is given by [44] ,r is the value of the current, once the temperature gradient is reversed, i.e. when the two baths' temperatures are interchanged. Notice that this coefficient takes values between 0 and 2, namely, R=0 when = - ,r , with the current being symmetric under reversing the temperature gradient. The upper bound is achieved when the current remains unaffected by reversing the temperature gradient. When either of the two currents is blocked, the coefficient is equal to one.

Main results
Before presenting our main results, let us summarize the major assumptions we made for our system, and the restrictions that these impose on the parameter regimes that we can consider: (i) Linearization of the impurity-bath coupling, which is achieved by assuming that the impurity is in the middle of its corresponding trap (see equation (6)). This in practice imposes a restriction on the maximum temperature we can consider [35] Note that in the scenario when the impurity is driven, W j 2 is replaced by (ii) BEC independence condition R L +R R <d.
In our analysis below we consider the BECs of Rubidium (Rb) atoms. The impurities are Dysprosium (Dy) atoms, which are the atoms with the largest magnetic moment known at present, μ=10μ B , where μ B is the Bohr magneton.

Heat current
In figure 2(a) we plot the heat current with the temperature difference ΔT=T L −T R , where we keep fixed the temperature of the left reservoir T L . As expected, the heat current increases with increasing temperature difference, while it is zero when there is no temperature gradient. We see that the current depends linearly on the temperature gradient. This is in accordance to the linear limit for the heat current in equation (60). In addition, we studied heat current in the scenario where the temperature difference between the two baths was fixed to some value, in particular ΔT=10 nK, with T L =T and T R =T+ΔT, and we considered the simultaneous variation of the temperatures of both baths, in the temperature regime T=10-100 nK. In this case we saw that the heat current remained constant as a function of T and hence we conclude that the regime in which we could study the system was that of the classical limit equation (61). A figure for this case is omitted since the current was just constant as a function of T. From figure 2(a) we also observe that increasing the distance of the two impurities results in decreasing the heat current flow (red versus blue curves), while increasing the impurity-BEC couplings results in an increase of the current (red versus green curves). Figure 2(b) depicts the heat current against the trapping frequencies of the BECs. In particular we fix one of the trapping frequencies and vary the other. Firstly, we observe that the heat current is reaches a maximum when the trapping frequencies of the two impurities match, i.e. Ω R = Ω L . This is understood as follows. The current density, which we define as den. We study the dependence on the impurity-BEC coupling strength in figure 2(c), where we see current reaches a maximum at some optimal coupling. Keeping the coupling constant of the left impurity fixed and varying that of the right impurity, we find that if the impurity is weakly coupled to the BECs, then current cannot be carried from one BEC to the other through the vibrations of these impurities. If on the other hand, this is coupled too strongly, the effect of the noise induced by the baths (BECs) reduces the current that can be transmitted This is in agreement with the findings in figures 2(a) and (b).
In figure 2(d) we see how the heat current varies as a function of the distance between the two impurities. As expected, we see that increasing the distance between the impurities reduces the heat current. In particular, we find k á ñ µ J 2 -in the parameter regime that we study. It is possible that, at shorter distances, another resonance effect occurs between the value of the spring constant, which depends inversely on the distance between the impurities and the trapping frequencies of the impurities. Anyhow, we do not study the system at such small distances because the approximation of independence of the BECs breaks down.
Finally we remark here that one could also consider studying homogeneous gases instead of harmonically trapped ones and the induced heat current could be examined for this case, by studying the system in the spirit of [73]. Experimentally, homogeneous gases could be created as in [74]. From our studies, we observe that as expected one can still have current in this case, but since the SD is different in this case, even though still superohmic, there is a quantitative difference in the amount of current. Nevertheless we focused on the harmonically trapped case which is more conventionally implemented experimentally. From [35] it is known that the two spectral densities result in a different degree of non-Markovianity, and it would be interesting in the future to study the effect of non-Markovianity on the heat current. This however exits the scope of this paper.

Current-current correlations
In figure 3(a) the behavior of the current-current correlations is illustrated as a function of temperature T while keeping the temperature difference ΔT constant. From the figure, we see that, for small ΔT, such that the first term of equation (65) prevailed, the current-current correlations are proportional to T 2 . On the contrary when ΔT is large, and for relatively small T, the current-current correlations depended linearly on the temperature T. Nevertheless the behavior seems to be independent of ΔT as the temperature increases, and appears to depend on the square of T as expected in the classical limit.
In figure 3(b) we study the dependence of current-current correlations on the temperature difference ΔT. At large temperature difference, i.e. beyond the linear limit considered in (64), the correlations depend linearly on ΔT. We comment here that this is not directly evident from the LogLog plot, but we checked that this is indeed the case of a linear relation from the corresponding current-current correlation versus temperature difference plot before taking the logarithms. For small ΔT, where (63) is well approximated by equation (64), indeed we find the saturation on the correlations predicted by the second term of equation (63). This implies that Figure 2. (a) Heat current J against the temperature difference of the two baths ΔT, with fixed T R . As expected current increases linearly with ΔT. Furthermore, the current decreases with increasing distance, and increases with increasing the coupling strength of the impurities to the bath. (b) Current as a function of the trapping frequency Ω R of the right impurity. We observe resonance at Ω R =Ω L . Furthermore, for Ω R larger than the trapping frequency of the bath-which is also the cutoff for the spectral density-the current vanishes quickly. Moreover, current decreases with distance as before. In this regime, current decreases as the impurities couple stronger to their respective baths (note that the coupling strengths are different from panel (a)). (c) Current as a function of the coupling strength of the right impurity. Current reaches a maximum in the range studied. (d) Current as a function of the distance between the impurities -| | d d 1 2 . Current decreases linearly with increasing distance in the regime we were allowed to study, that is under the restriction that is imposed on the lower distance in order to maintain the two BECs spatially independent. Current is found to scale as κ 2 . See table 1 for the parameters that we use here. even as D  T 0-in which case current is zero-the fluctuations are still present. One might expect quantum effects to appear in this regime then, but we studied the entanglement in this system, by means of the logarithmic negativity as in [49], and no entanglement could be detected in this regime.

Driven case: heat rectification
Heat rectification is quantified by the heat rectification, equation (69). We use the particular form of driving In figure 4 we depict the heat rectification coefficient R as a function of W d , the driving frequency. We find that it shows to maxima at for two values of the driving frequency Ω d ä {Ω, 3Ω}. Note that, as shown in [44], heat rectification should be maximum at the following frequencies n n W =  | | ( ) . 7 2 d j i Figure 3. (a) Current-current correlations against temperature T at constant DT . Here we observe that for small ΔT and at large temperatures T as expected from equation (65), current-current correlations scale as ∝T 2 . At sufficiently large T this is expected to happen for large ΔT as well, however we were restricted on the range of maximum temperatures we could consider. (b) Current-current correlations as a function of the temperature difference of the two baths ΔT. For small values of ΔT, i.e. in the linear regime, the correlations have a nonzero value as predicted in equation (64). This implies that even as D  T 0, where current also vanishes, correlations still persist. For large enough values of ΔT the correlations increase linearly with ΔT. See table 1 for the parameters that we use here. Table 2. List of default parameters-unless otherwise mentioned-used in the dynamic case, corresponding to figure 4.

Left BEC and impurity
Right BEC and impurity Other parameters   figure 4. Note that we also studied dependence of the rectification coefficient on the other parameters of the system, apart from the driving frequency. In particular, we find that in the regime of parameters we study, maximum rectification decreases when the impurities couple more strongly to their respective baths. On the contrary, decreasing the number of atoms significantly increases R, in fact reaching R>1. What is more, rectification could also be optimized with respect to the detuning between the trapping frequencies of the impurities as in [44].

Conclusions
In this work, we studied in detail the heat transport control and heat current rectification among two BEC. To this aim, we took an open quantum system approach, and focused on experimentally realistic conditions and parameter regimes. In particular we considered a system composed of two harmonically trapped interacting impurities immersed in two independent harmonically trapped 1D BECs kept at different temperatures. The impurities interact through a long range interaction, in particular a dipole-dipole coupling-that under suitable conditions we were able to treat as a spring-like interaction. In this work we considered the particular case of a fixed angle, motivated by the results in [39], but the results should also hold under an angle-averaged scenario, by taking advantage of the results in( [41]). We showed the dynamics of these impurities can be described within the framework of quantum Brownian motion, where the excitation modes of the gas play the role of the bath. In this analogy, the SD of the bath is not postulated, but it is rather derived exactly from the Hamiltonian of the BEC, which turns out to be superohmic. By solving the relevant generalized Langevin equations, we find the steady state covariance matrix of the impurities, which contains all the information describing our Gaussian system. In particular, we use such information to study the heat currents and current-current correlations and their dependence on the controllable parameters of the system. We find that, the heat current scales linearly with temperature difference among the two BECs. Furthermore, we observe that heat current is maximum when the trapping frequencies of the impurities are at resonance. Finally, we showed the existence of an optimal coupling strength of the impurities on their respective baths.
What is more, by periodically driving one of the impurities, we can conduct heat asymmetrically, i.e. we achieve heat rectification-which is in full agreement with the recent proposal of [44]. In particular, we see that one can achieve heat rectification at the driving frequencies predicted in [33], even though our bath is superohmic.
Motivated by recent developments on the usage of BECs as platforms for quantum information processing, as e.g. in [49], our work offers an alternative possibility to use this versatile setting for information transfer and processing, within the context of phononics. The possibility of quantum advantages using many-body impurities in our platform remains an interesting open question (see [27] too). Another future direction is to study heat control in 2D and 3D BECs. In principle this gives rise to a different SD, which opens a new window for further Figure 4. Rectification coefficient against the driving frequency. As predicted analytically, we observe nonzero rectification at the vicinity of frequencies Ω d ={Ω, 3Ω} corresponding to n n  | | j i . Furthermore, we see that rectification decreases with increasing coupling to the baths in this regime of parameters. See table 2 for the parameters that we use here. manipulation of heat current. Finally, it is desirable to investigate scenarios where the system is nonlinear, which raises difficulties in solving the problem analytically, nonetheless, it offers the opportunity to rectify heat even without periodic derivation. Moreover, motivated by the results in [75,76] where the squeezing in position of a single impurity embedded in a BEC was used to measure the temperature of the BEC in the sub-nano-Kelvin regime, one may study if the present two-particle set-up can be used for applications in quantum thermometry.

PX PP
where = ( ( ) ( )) p t p t P , We emphasize that the state of the system-bath is assumed to be a product state as in (35). Hence the average is taken over the thermal state of the bath, while the state of the system is assumed to have reached its unique equilibrium state by considering the long-time, steady state limit  ¥ t which is equivalent to consider that w L L  , where ν − is the minimum standard eigenvalue of~( ) C 0 (ÿ is assumed to be equal to 1 in this case). Static covariance matrix. The elements of the covariance matrix in the static case, in terms of the phonons propagator functions, can be expressed as