Heat current control in trapped BEC

TL;DR

This paper explores heat transport control between two trapped Bose-Einstein Condensates connected by impurities with dipole interactions, demonstrating heat rectification and potential for phononic quantum information processing.

Contribution

It introduces a novel setup using long-range dipole interactions and periodic driving to achieve heat rectification in coupled BECs, advancing quantum thermodynamics and phononics.

Findings

Heat current depends on system parameters.

Heat rectification occurs with periodic driving.

System can function as a phononic circuit.

Abstract

We investigate the heat transport and the control of heat current among two spatially separated trapped Bose-Einstein Condensates (BEC), 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.

Tables

Table 1. Table 1: List of default parameters—unless otherwise mentioned—used in the static case, corresponding to Figures.  2 and   3 .

Left BEC and impurity	Right BEC and impurity	Other parameters
$T_L=75nK$	$T_R=7.5nK$	$d=9.5a_{\mathrm{ho}}$
$N_L=7.5\times {10}^4$	$N_R=7\times {10}^4$	$a_{\mathrm{ho}}=0.7\mu m$
${\eta }_L=0.5$	${\eta }_R=0.5$	$\Omega =200\pi Hz$
$g_B^{\left(L\right)}=2.5\times {10}^{-40}J\cdot m$	$g_B^{\left(R\right)}=2\times {10}^{-40}J\cdot m$
${\Omega }_L=\Omega$	${\Omega }_R=\Omega$
${\hat{\Omega }}_{B,L}=3\Omega$	${\hat{\Omega }}_{B,R}=3\Omega$

Table 2. Table 2: List of default parameters—unless otherwise mentioned—used in the dynamic case, corresponding to Fig.  4

Left BEC and impurity	Right BEC and impurity	Other parameters
$T_L=15nK$	$T_R=1nK$	$d=35a_{\mathrm{ho}}$
$N_L=7.5\times {10}^4$	$N_R=7\times {10}^4$	$a_{\mathrm{ho}}=0.7\mu m$
${\eta }_L=0.5$	${\eta }_R=0.5$	$\Omega =200\pi Hz$
$g_B^{\left(L\right)}=5\times {10}^{-39}J\cdot m$	$g_B^{\left(R\right)}=4.5\times {10}^{-39}J\cdot m$	$v=0.1\Omega$
${\Omega }_L=2\Omega$	${\Omega }_R=2\Omega$
${\Omega }_{B,L}=2.5\Omega$	${\Omega }_{B,R}=2.5\Omega$

Equations

$$H_{\mathrm{Tot}}=H_{\mathrm{S}}+\sum _{j\in \left\{\mathrm{L},\mathrm{R}\right\}}\left(H_{\mathrm{B},\mathrm{j}}+H_{\mathrm{BB},\mathrm{j}}\right)+\sum _{j\in \left\{\mathrm{L},\mathrm{R}\right\}}{H}_{Int.,\mathrm{j}}+H_{\mathrm{Drive}}\left(t\right),$$

$$H_{\rm{S}}=\sum_{\begin{array}[]{c}i,j\in\left\{\rm{L,R}\right\},\\ i\neq j\end{array}}\frac{\mathbf{p}_{j}^{2}}{2m}+\mathbf{V}_{j}\left(\mathbf{x}_{j}\right)+H_{\rm{int.}},$$

$$H_{\mathrm{B},\mathrm{j}}=\int {\mathbf{\Psi }}_j^{†}\left({\mathbf{x}}_{\mathrm{B}}\right)\left[\frac{{\mathbf{q}}_j^2}{2m_{\mathrm{B}}}+{\mathbf{V}}_{\mathrm{B}}^{\left(j\right)}\left({\mathbf{x}}_{\mathrm{B}}\right)\right]{\mathbf{\Psi }}_j\left({\mathbf{x}}_{\mathrm{B}}\right)\mathrm{d}{\mathbf{x}}_{\mathrm{B}},$$

$$H_{\mathrm{BB},\mathrm{j}}=g_{\mathrm{B}}^{\left(j\right)}\int {\mathbf{\Psi }}_j^{†}\left({\mathbf{x}}_{\mathrm{B}}\right){\mathbf{\Psi }}_j^{†}\left({\mathbf{x}}_{\mathrm{B}}\right){\mathbf{\Psi }}_j\left({\mathbf{x}}_{\mathrm{B}}\right){\mathbf{\Psi }}_j\left({\mathbf{x}}_{\mathrm{B}}\right)\mathrm{d}{\mathbf{x}}_{\mathrm{B}},$$

$$\begin{array}[]{c}H_{\rm{Int.},j}=g_{\rm{IB}}^{(j)}\int\mathbf{\Psi}_{j}^{\dagger}\left(\mathbf{x_{\rm{B}}}\right)\mathbf{\Psi}_{j}\left(\mathbf{x_{\rm{B}}}\right)\delta\left(\mathbf{x}_{j}-\mathbf{x_{\rm{B}}}\right)\rm{d}\mathbf{x_{\rm{B}}},\\ =g_{\rm{IB}}^{(j)}\mathbf{\Psi}_{j}^{\dagger}\left(\mathbf{x}_{j}\right)\mathbf{\Psi}_{j}\left(\mathbf{x}_{j}\right),\end{array}$$

$$x_j\ll R_j,$$

$$H_{\mathrm{Tot}}=H_{\mathrm{S}}\left(t\right)+\sum _{j\in \left\{\mathrm{L},\mathrm{R}\right\}}{\tilde{H}}_{\mathrm{B},\mathrm{j}}+\sum _{j\in \left\{\mathrm{L},\mathrm{R}\right\}}{\tilde{H}}_{Int.,\mathrm{j}},$$

$${\tilde{H}}_{\mathrm{B},\mathrm{j}}=\sum _q{E}_q^{\left(j\right)}{b}_{j,q}^{†}{b}_{j,q},$$

$${\tilde{H}}_{Int.,\mathrm{j}}=\sum _q\hslash g_q^{\left(j\right)}{x}_j\left(b_{j,q}+b_{j,q}^{†}\right),$$

$$E_q^{\left(j\right)}=\hslash {\Omega }_{\mathrm{B},\mathrm{j}}\sqrt{q\left(q+1\right)},$$

$$g_q^{\left(j\right)}=\frac{g_{\mathrm{IB}}^{\left(j\right)}{\mu }_j}{\hslash {\pi }^{3/2}}{\left[\frac{1+2q}{\hslash {\Omega }_{\mathrm{B},\mathrm{j}}{g}_{\mathrm{B}}^{\left(j\right)}{R}_j^3}\right]}^{1/2}\frac{\Gamma \left[\frac{1}{2}\left(1-q\right)\right]\Gamma \left[\frac{1}{2}\left(1+q\right)\right]\sin \left(\pi q\right)}{{\left[q\left(q+1\right)\right]}^{1/4}}.$$

$${\mu }_j={\left(\frac{3}{4\sqrt{2}}{g}_{\mathrm{B}}^{\left(j\right)}{N}_j{\Omega }_{\mathrm{B},\mathrm{j}}\sqrt{m_{\mathrm{B}}}\right)}^{2/3},$$

$$R_j=\sqrt{\frac{2{\mu }_j}{m_{\mathrm{B}}{\left({\Omega }_{\mathrm{B},\mathrm{j}}\right)}^2}}.$$

$$H_{int.}=A\frac{1}{r^3},$$

$$C_{dd}={\mu }^2{\mu }_0,$$

$$H_{int.}=A{d}^{-3}{\left|1+\frac{\left(x_{\mathrm{L}}-x_{\mathrm{R}}\right)}{d}\right|}^{-3},$$

$$H_{int.}=A{d}^{-3}\left(1-3\frac{\left(x_{\mathrm{L}}-x_{\mathrm{R}}\right)}{d}+6\frac{{\left(x_{\mathrm{L}}-x_{\mathrm{R}}\right)}^2}{d^2}\right).$$

$${\hat{H}}_{int.}=6A{d}^{-5}\left(x_{\mathrm{L}}^2-2x_{\mathrm{L}}{x}_{\mathrm{R}}+x_{\mathrm{R}}^2\right),$$

$$H_{\mathrm{S}}^{\mathrm{non}-inter.}=\frac{p_{\mathrm{L}}^2}{2m}+\frac{p_{\mathrm{R}}^2}{2m}+\frac{m}{2}{\hat{\Omega }}_{\mathrm{L}}^2{x}_{\mathrm{L}}^2+\frac{m}{,}2{\hat{\Omega }}_{\mathrm{R}}^2{x}_{\mathrm{R}}^2$$

$$\begin{array}[]{c}\hat{\Omega}_{\rm{L}}^{2}=\Omega_{\rm{L}}^{2}+6Ad^{-5},\\ \hat{\Omega}_{\rm{R}}^{2}=\Omega_{\rm{R}}^{2}+6Ad^{-5}.\end{array}$$

$${\tilde{H}}_{int.}=\kappa x_1{x}_2,$$

$$\kappa =\frac{12C_{dd}\left(1-3{\left(\cos \theta \right)}^2\right)}{4\pi d^5}.$$

$$H_{\mathrm{Drive}}\left(t\right)=\Theta \left(t-t_0\right)\frac{1}{2}{\mathbf{X}}^T\cdot \mathbf{V}\left(t\right)\mathbf{X},$$

$$\mathbf{V}\left(t\right)=\sum _k{\mathbf{V}}_k{e}^{ik{\Omega }_{d}t},$$

$$\frac{d{b}_{k,j}\left(t\right)}{dt}$$

$$\frac{d{b}_{k,j}^{†}\left(t\right)}{dt}$$

$$\frac{d{x}_j\left(t\right)}{dt}=\frac{i}{\hslash }\left[H_{\mathrm{Tot}},x_j\left(t\right)\right]=\frac{p_j\left(t\right)}{m},$$

$$\frac{d{p}_j}{dt}\left(t\right)=\frac{i}{\hslash }\left[H_{\mathrm{Tot}},p_j\left(t\right)\right]$$

$$=-m({\Omega }_j^2+\sum _k{V}_{k,jj}{e}^{ik{\omega }_{d}t})x_j\left(t\right)$$

$$-(\kappa +\sum _k{V}_{k,jj}{e}^{ik{\omega }_{d}t})x_q\left(t\right)-\hslash \sum _k{g}_{k,j}\left(b_{k,j}\left(t\right)+b_{k,j}^{†}\left(t\right)\right),$$

$$b_{kj}\left(t\right)$$

$$b_{kj}^{†}\left(t\right)$$

$${\ddot{x}}_j+({\Omega }_j^2+\sum _k{V}_{k,jj}{e}^{ik{\omega }_{d}t})x_j+(\kappa +\sum _k{V}_{k,jq}{e}^{ik{\omega }_{d}t})x_q-{\int }_0^t{\lambda }_j\left(t-s\right)x_j\left(s\right)\mathrm{d}\mathrm{s}=\frac{1}{\mathrm{m}}{\mathrm{B}}_{\mathrm{j}}\left(\mathrm{t}\right),$$

Full text

1 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, 2, 3, 4, 21, 23, 24, 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 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 Bose-Einstein condensates (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 long-range dipolar interactions (see Fig. 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 out-of-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.

2 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 $\Omega_{\rm{B,L}}$ and $\Omega_{\rm{B,R}}$, respectively, but one trapping potential is centered around a minima located in $\mathbf{d}_{\rm{L}}$ and the other around $\mathbf{d}_{\rm{R}}$, with certain distance among them. Each BEC has, respectively, $N_{\rm{L}}$ and $N_{\rm{R}}$ interacting atoms of the same species with mass $m_{\rm{B}}$. The two impurities are also trapped in a parabolic potential of frequencies $\Omega_{\rm{L}}$ and $\Omega_{\rm{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

[TABLE]

where

[TABLE]

[TABLE]

[TABLE]

[TABLE]

where $\mathbf{x}$ and $\mathbf{x_{\rm{B}}}$ are the three-dimenional 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_{\rm{B}}^{\left(j\right)}$ and $g_{\rm{IB}}^{(j)}$ respectively. Here, $H_{\rm{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_{\rm{L}}=-d/2$ and $d_{\rm{R}}=d/2$.

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 Bogoliubov-de-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}$,

[TABLE]

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

[TABLE]

where

[TABLE]

[TABLE]

with

[TABLE]

and

[TABLE]

The chemical potential $\mu_{j}$ of the $j^{\rm{th}}$ bath is

[TABLE]

and the Thomas-Fermi radius $R_{j}$ of the $j^{\rm{th}}$ bath is

[TABLE]

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\left(1\right)$ 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:

[TABLE]

where $A=\frac{C_{dd}\left(1-3\left(\cos\theta\right)^{2}\right)}{4\pi}$, with

[TABLE]

and $\mu_{0}=4\pi 10^{-7}\frac{J}{A^{2}m}$ being the vacuum permeability, while $\mu$ is the magnetic moment of the dipole. The angle $\theta$ is that formed between the axis of the two dipoles and it can determined in experiments. Furthermore, in (14), $r=\left|\left(d_{\rm{R}}+x_{\rm{R}}\right)-\left(d_{\rm{L}}+x_{\rm{L}}\right)\right|$ is the distance between the two oscillators, with $d_{\rm{L}}$ and $d_{\rm{R}}$ being the centers of the trapping potentials of the two oscillators. We then rewrite (14) as

[TABLE]

and we assume that the distance between the two oscillators centers $d$ is much larger than the fluctuations $x_{\rm{L}},x_{\rm{R}}$ and hence the fluctuations difference $x_{\rm{L}}-x_{\rm{R}}$, i.e. $x_{\rm{L}},x_{\rm{R}}\ll d$. 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_{\rm{L}}+R_{\rm{R}}<d$. With the additional assumption we made before, namely that the impurities oscillations are much smaller than their corresponding Thomas-Fermi radius, Eq. (6), then one concludes that $x_{\rm{L}},x_{\rm{R}}\ll d$. 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 $x_{\rm{L}},x_{\rm{R}}\ll d$ should be made explicitly.

Finally, after expanding the binomial series, Eq. (16) is rewritten as

[TABLE]

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_{\rm{L}},x_{\rm{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

[TABLE]

and we absorb the terms with $x_{\rm{L}}^{2}$ and $x_{\rm{R}}^{2}$ in the non-interacting part of the Hamiltonian by redefining the frequencies as

[TABLE]

where

[TABLE]

From here on we omit the tilde in the frequencies to avoid unnecessary complications in the nomenclature. Therefore we can rewrite the interaction as

[TABLE]

which models a spring-like interaction among the two impurities with

[TABLE]

The angle $\theta$ can be experimentally controlled as in [39]. 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 can not 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 $x_{\rm{L}},x_{\rm{R}}\ll d$, qualitativley the results will be the same, with the difference being that the distance at which the oscillators should be kept will change [40]. 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, 41]. 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 $H_{\rm{Drive}}\left(t\right)=\Theta\left(t-t_{0}\right)\mathbf{f}^{T}\left(t\right)\mathbf{X}\left(t\right)$ where $\mathbf{X}\left(t\right)=\left(x_{1}\left(t\right),x_{1}\left(t\right)\right)^{\intercal}$ and $\mathbf{f}\left(t\right)$ is some periodic function, e.g. $\mathbf{f}\left(t\right)=\mathbf{f}_{0}e^{-i\omega_{d}t}+c.c.$ with $t_{0}$ the time at which the driving begins, $\Omega_{d}$ the driving frequency and $\mathbf{f}_{0}$ a complex valued constant column vector, $\mathbf{f}_{0}=\left(f_{0,\rm{L}},f_{0,\rm{R}}\right)^{\intercal}$. This type of driving was considered in [42], 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.

[TABLE]

where the driving is either on the trapping frequency of the central oscillators or on their in-between coupling. It was recently shown in [43] 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 [44], by introducing a coherently driven coupling between the two oscillators. We also assume the driving to be periodic, $\mathbf{V}\left(t+\tau\right)=\mathbf{V}\left(t\right)$ with $\tau$ being the time period, such that it can be Fourier expanded as

[TABLE]

where $\Omega_{d}=1/\tau$ being the driving frequency. This type of coupling could also be implemented by a laser.

3 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

[TABLE]

and impurity particles

[TABLE]

where $j,q\in\left\{\rm{L,R}\right\}$ and $j\neq q$. We first solve the bath particles equations of motion

[TABLE]

and we replace these in the impurities equations of motion (26)–(27), to obtain

[TABLE]

where $B_{j}\left(t\right)$ plays the role of the stochastic force and reads as

[TABLE]

Here $\lambda_{j}\left(t\right)$ 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

[TABLE]

with

[TABLE]

being the spectral density. Equation (31) has the form of a Generalized Langevin equation which describes the evolution of a system with memory and under the influence of a stochastic force. These terms, $B_{j}\left(t\right)$ and $\gamma_{j}\left(t\right)$, 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,

[TABLE]

where $\rho\left(0\right)$ is the total density state, $\rho_{\rm{S}}\left(0\right)$ is the initial density state of the system and $\rho_{\rm{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_{j}\left(t\right)$ obeys the fluctuation dissipation relation

[TABLE]

where $\lambda_{j}\left(\omega\right)$ is the Fourier transform of $\lambda_{j}\left(t\right)$ and obeys [45]

[TABLE]

Hence one concludes that, upon determination of the spectral density $J_{j}\left(\omega\right)$, one determines the influence of the baths on the impurities. In [35], it was shown that in the continuous frequency limit, the spectral density takes the following form

[TABLE]

where

[TABLE]

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 $j^{th}$ bath. Here $\frac{\tau_{j}}{\Lambda_{j}^{3}}$ plays the role of a relaxation time. Heat transport with a superohmic spectral density was considered for example in the energy transport in the phenomenon of photosynthesis [46]. Note that Eq. (31) can be rewritten in terms of the damping kernel $\gamma_{j}\left(t\right)$ which is related to the susceptibility by $\frac{1}{m}\lambda_{j}\left(t\right)=-\frac{\partial}{\partial t}\gamma_{j}\left(t\right)$ [47], as

[TABLE]

with $\tilde{\Omega}_{j}^{2}=\Omega_{j}^{2}+\sum_{k}V_{k,jj}e^{ik\omega_{d}t}$ and $\tilde{\kappa}_{q}=\kappa+\sum_{k}V_{k,jq}e^{ik\omega_{d}t}$, where the Leibniz rule was used as in  [48]. Moreover, in [35], it was also shown that the form of the damping kernel for such a spectral density reads as

[TABLE]

In the limit $t\rightarrow 0$ this damping kernel becomes

[TABLE]

In vector form the two coupled equations in (31) read as

[TABLE]

where $\mathbb{I}$ is the identity matrix and

[TABLE]

Static case.

From the system of coupled equations (31), one can now identify the condition that guarantees the positivity of the Hamiltonian in the static case, where $H_{\rm{Drive}}\left(t\right)=0$ . 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 $\mathbf{Q}=\mathbf{O}\cdot\mathbf{X}$, that brings (43) into

[TABLE]

where $\mathbf{D}_{D}\left(t\right)=\mathbf{O}\cdot\mathbf{D}\left(t\right)\cdot\mathbf{O}^{T}$, $\mathbf{B}_{D}^{T}\left(t\right)=\mathbf{O}\cdot\mathbf{B}^{T}\left(t\right)$, and the frequency matrix is diagonalized as

[TABLE]

The positivity of the Hamiltonian condition is then guaranteed by requiring that $\{\Omega_{\rm{L}}^{D},\Omega_{\rm{R}}^{D}\}>0$ which in terms of the original frequencies reads as

[TABLE]

This condition guarantees that we do not have negative renormalized normal frequencies in the system and hence stability of the solution is guaranteed in the long-time limit. With this satisfied, we are in a position to safely neglect the effects of $\gamma_{\rm{L}}\left(0\right)$ and $\gamma_{\rm{R}}\left(0\right)$ in the dynamics of the system.

Upon rewriting the coupled equations of motion for the static case in a vector form as above, and considering it expressed in terms of the susceptibilities $\lambda_{\rm{L}}\left(t\right),\lambda_{\rm{R}}\left(t\right)$, one can obtain its solution by taking the Fourier transform of both sides

[TABLE]

where

[TABLE]

which is understood to play the role of a phonon Green function, and $\mathbf{L}\left(\omega\right)$ is the Fourier transform of \mathbf{L}\left(t\right)=\left(\begin{array}[]{cc}\lambda_{\rm{L}}\left(t\right)&0\\ 0&\lambda_{\rm{R}}\left(t\right)\end{array}\right), with diagonal elements

[TABLE]

where the real part of the susceptibility was obtained through the Kramers-Kronig relation $Re[\lambda_{j}\left(\omega^{\prime}\right)]=\mathcal{H}[Im[\lambda_{j}\left(\omega\right)]]\left(\omega^{\prime}\right)=\frac{1}{\pi}P\int_{-\infty}^{\infty}\frac{Im[\lambda_{j}\left(\omega\right)]}{\omega-\omega^{\prime}}\rm{d}\omega$. Here $\mathcal{H}[\cdot]\left(\omega^{\prime}\right)$ denotes the Hilbert transform and $P$ the principal value. In general, for the parameters we consider in our results, it will always hold that $\Omega_{j}^{2}\gg Re[\lambda_{j}\left(\omega\right)]$ such that we safely neglect the effect of $Re[\lambda_{\rm{L}}\left(\omega\right)]$ and $Re[\lambda_{\rm{R}}\left(\omega\right)]$.

In terms of this solution of the static equations of motion, the positivity condition (49), can be interpreted in a different way: it guarantees that the phonon propagator, i.e. the Green function $\mathbf{G}\left(\omega\right)$, 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 [49].

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 [43]. 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\left(t\right)$ might not be periodic, another function defined based on this one as $\mathbf{W}\left(t,\omega\right)=\int_{\mathbb{R}}dt^{\prime}f\left(t-t^{\prime}\right)e^{i\omega\left(t-t^{\prime}\right)}$ will indeed be periodic. Furthermore, such periodic function can always be expressed in terms of its Fourier components as $\mathbf{W}\left(t,\omega\right)=\sum_{k}\mathbf{B}_{k}\left(\omega\right)e^{ik\Omega_{d}t}$. By performing these two transformations on the equations of motion, i.e. $\int_{\mathbb{R}}dt^{\prime}e^{i\omega\left(t-t^{\prime}\right)}$ and Fourier expanding, one obtains a set of equations for the Fourier coefficients $\mathbf{A}_{k}\left(\omega\right)$ of $\mathbf{P}\left(t,\omega\right)=\int_{\mathbb{R}}dt^{\prime}\mathbf{G}\left(t-t^{\prime}\right)e^{i\omega\left(t-t^{\prime}\right)}=\sum_{k}\mathbf{A}_{k}\left(\omega\right)e^{ik\Omega_{d}t}$ that can be selfconsistently solved. By following this procedure, in [43] the authors were able to obtain expressions for these coefficients

[TABLE]

where $\mathbf{V}\left(t\right)=\sum_{k}\mathbf{V}_{k}e^{ik\Omega_{d}t}$. Note that we will assume that the driving strength coefficient is sufficiently small such that we can ignore terms of the order of $\mathcal{O}\left(\mathbf{V}_{j}^{3}\right)$ or higher. Furthermore, since the Fourier coefficients $\mathbf{A}_{k}\left(\omega\right)$ 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 $\omega$ of the phonon by an amount of $k\Omega_{d}$. Finally, the solution of the equations of motion in this case would read as

[TABLE]

Uncertainty relation.

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

[TABLE]

where $\nu_{-}$ is the minimum standard eigenvalue of $\widetilde{\mathbf{C}}\left(0\right)$ ($\hbar$ is assumed to be equal to 1 in this case). In the above expression, $\widetilde{\mathbf{C}}\left(0\right)=i\mathbf{W}\cdot\mathbf{C}\left(0\right)$ with $\mathbf{W}$ the symplectic matrix and $\mathbf{C}\left(0\right)$ 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 spectral density, for both static and driven cases as is shown in the appendix.

4 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 spectral density that describes the baths. As such, the results in next section are obtained by numerically evaluating the integrals.

4.1 Static case

We begin with the static scenario with $H_{\rm{Drive}}\left(t\right)=0$, and we study the behaviour of heat current and the current-current correlations.

The heat current.

When studying a heat engine, a key quantity is the average current $J_{\rm{L}}$ going from the left reservoir (assuming it to be the hot reservoir) to the left oscillator [16, 50] (by conservation of current $J_{\rm{L}}=-J_{\rm{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,

[TABLE]

where $H^{\prime}_{\rm{L}}=\frac{p_{\rm{L}}^{2}}{2m}+V_{\rm{L}}\left(x_{\rm{L}}\right)+\kappa x_{\rm{L}}x_{\rm{R}}$ and the average is over the total density state $\rho$. The second definition is expressed in terms of the rate of decrease of the bath energy [51]

[TABLE]

The second term vanishes under steady state conditions and weak system-bath coupling, where weak is understood in the sense of assumption (35) (and not of Markovian dynamics for the impurity), which is the case in our study. We are further considering that the correlations among the system bath interactions is negligible, i.e. $\left\langle\left[\widetilde{H}_{Int.,L},\widetilde{H}_{Int.,R}\right]\right\rangle_{\rho_{\rm{B}}}=0$, 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 [52], while the case of interactions among the baths was studied in [12]. Therefore the heat current considered here is

[TABLE]

where the second equality is valid in steady state with $\eta_{\rm{L}}\left(t\right)=\int_{0}^{t}\lambda_{\rm{L}}\left(t-s\right)x_{\rm{L}}\left(s\right)\rm{d}s+\frac{1}{m}B_{\rm{L}}\left(t\right)$, and it is obtained after solving the bath particles equations of motion. Note that in (59) 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, 50, 53], by using a direct solution of the equations of motion for a non-interacting system of bath particles, the average current is proven to be equal to

[TABLE]

with $f\left(\omega,T_{j}\right)=[{e^{\frac{\hbar\omega}{k_{\rm{B}}T_{j}}}-1}]^{-1}$, the phonon occupation number for the respective thermal reservoir and

[TABLE]

being the transmission coefficient for phonons at frequency $\omega$. This is called the Landauer formula. Note that $\mathcal{T}\left(\omega\right)$ can be related to the transmittance of plane waves of frequency $\omega$ across the system as in [54]. 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 (60) 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 [55], 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 [56]. 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 [57].

It is interesting to comment on two limits of (60), namely the linear limit $\Delta T:=T_{\rm{L}}-T_{\rm{R}}\ll T$ where $T=\frac{T_{\rm{L}}+T_{\rm{R}}}{2}$ and the classical limit, $\frac{\hbar\omega}{k_{\rm{B}}T}\rightarrow 0$ where $\omega$ refers to the bath frequencies. In the first limit, the current reduces to

[TABLE]

such that once the two baths are at the same temperature there is no current flow. In the second limit, it becomes

[TABLE]

where the current is independent of the temperature of the baths $T_{\rm{L}},T_{\rm{R}}$ but it only depends on their difference $\Delta 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 [43], 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 [58, 59] 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, 50] and the local effective temperature of driven systems [60].

The current-current time correlations $JJ_{\rm{LL}}\left(t,t^{\prime}\right)$, 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

[TABLE]

We are interested in the steady state correlations, for which an expression for this can be obtained using the non-equilibrium Green’s functions as was mentioned above. Hence the correlation function of interest is a function only of the time difference, $JJ\left(t,t^{\prime}\right)=JJ\left(t-t^{\prime}\right)$. Therefore, the noise strength is characterized by the zero frequency component $JJ_{\alpha\beta}=\int_{-\infty}^{\infty}JJ_{\alpha\beta}\left(t\right)\rm{d}t$, which obeys $JJ_{\alpha\beta}\geq 0$ according to Wiener-Khinchine theorem. It is current conserving, i.e. the sum of currents entering the system from all reservoirs is equal to zero at each instant of time, and gauge invariant and hence physically meaningful [61, 62]. Current conservation implies $JJ_{\rm{LL}}=JJ_{\rm{RR}}$.

One way to obtain such an expression is by first deriving the cumulant generating function $\chi\left(\mu\right)$, employing the non-equilibrium Green’s functions technique within the Keldysh formalism, and noting that $JJ_{\alpha\beta}:=\frac{\partial^{2}\chi\left(\mu\right)}{\partial\mu^{2}}\mid_{\mu=0}$ (see e.g.,  [16, 50, 63]). In this case, one can show that the current fluctuations read as [61, 62]

[TABLE]

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 Eq. (70) 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 current-current correlations read as

[TABLE]

In the classical limit it becomes

[TABLE]

Note that (71), contrary to the expression for the current (67) at the same limit, does not vanish when $\Delta T\rightarrow 0$, which results in nonzero fluctuations of the current even in the scenario that no average current flows in the system.

4.2 The dynamic case

For the driven case, the steady-state averaged heat current is given by [43]

[TABLE]

where the new transmission coefficient reads as

[TABLE]

and

[TABLE]

This expression was also obtained in [62, 64, 65, 66], while in [67] a similar expression was obtained for transport through quantum dots. Unlike the static case, these transmission coefficients are not symmetric, i.e., $\widetilde{T}_{jq}\left(\omega\right)\neq\widetilde{T}_{qj}\left(\omega\right)$. Crucially, this symmetry breaking, attributed to the driving that is now expressed in the form of the transmission coefficients, is responsible for the appearance of heat rectification as addressed in [43]. To observe and quantify rectification, it is useful to evaluate the rectification coefficient

[TABLE]

where $J_{j,\rm{r}}^{\left(D\right)}$ 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 $J_{j}^{\left(D\right)}=-J_{j,\rm{r}}^{\left(D\right)}$, 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.

5 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:

Linearization of the impurity-bath coupling, which is achieved by assuming that the impurity is in the middle of its corresponding trap (see Eq. (6)). This in practice imposes a restriction on the maximum temperature we can consider [35]

[TABLE]

Note that in the scenario when the impurity is driven, $\Omega_{j}^{2}$ is replaced by $\min_{t}\left(\Omega_{j}^{2}+\sum_{k}V_{k,jj}e^{ik\Omega_{d}t}\right)$. 2. 2.

BEC independence condition $R_{\rm{L}}+R_{\rm{R}}<d$ 3. 3.

Positivity condition, Eq. (49).

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, $\mu=10\mu_{\rm{B}}$, where $\mu_{\rm{B}}$ is the Bohr magneton.

5.1 Static system

5.1.1 Heat current.

In Fig. 2 (a) we plot the heat current with the temperature difference $\Delta T=T_{\rm{L}}-T_{\rm{R}}$, where we keep fixed the temperature of the left reservoir $T_{\rm{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 Eq. (67). In addition, we studied heat current in the scenario where the temperature difference between the two baths was fixed to some value, in particular $\Delta T=10nK$, with $T_{\rm{L}}=T$ and $T_{\rm{R}}=T+\Delta T$, and we considered the simultaneous variation of the temperatures of both baths, in the temperature regime $T=10{\rm nK}-100{\rm 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 Eq. (68). A figure for this case is omitted since the current was just constant as a function of $T$. From Fig. 2 (a) we also observe that increasing the distance of the two impurities results in decreasing the heat current flow (red vs blue curves), while increasing the impurity-BEC couplings results in an increase of the current (red vs 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., $\Omega_{\rm R}=\Omega_{\rm L}$. This is understood as follows. The current density, which we define as $J_{den.}:=\frac{1}{4\pi}\mathcal{T}\left(\omega\right)\hbar\omega\left[f\left(\omega,T_{\rm{L}}\right)-f\left(\omega,T_{\rm{R}}\right)\right]$ is maximized whenever the denominator of the transmission function $\mathcal{T}\left(\omega\right)$, given by $\left(-\omega^{2}\mathbf{I}+\mathbf{K}-\mathbf{L}\left(\omega\right)\right)^{-1}\left(\left(-\omega^{2}\mathbf{I}+\mathbf{K}-\mathbf{L}\left(-\omega\right)\right)^{T}\right)^{-1}$ is minimized. In the regime we are looking, $\frac{\tau_{j}}{\widehat{\Omega}_{j}^{3}}$ come out to be of the order of $10^{-4}$, while the values of $\kappa$ that are allowed, are of the order of $10^{-5}$. These are much smaller than the trapping frequencies $\Omega_{\rm{L}}^{2}$, $\Omega_{\rm{R}}^{2}$, such that the denominator is minimized whenever $\left(\omega^{2}-\Omega_{L}^{2}\right)^{2}+\left(\omega^{2}-\Omega_{R}^{2}\right)^{2}$ is minimized. This happens when $\omega=\Omega_{L}=\Omega_{R}$. Secondly—for the specific parameters that we choose—contrary to Fig. 2 (a), increasing the impurity-BEC coupling strength results in reducing the current.

We study the dependence on the impurity-BEC coupling strength in Fig. 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 can not 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 Fig. 2 (a) and (b).

In Fig. 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 $\left\langle J\right\rangle\propto\kappa^{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 [68]. Experimentally, homogeneous gases could be created as in [69]. From our studies, we observe that as expected one can still have current in this case, but since the spectral density 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 densitites 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.

5.1.2 Current-Current correlations.

In Fig. 3 (a) the behavior of the current-current correlations is illustrated as a function of temperature $T$ while keeping the temperature difference $\Delta T$ constant. From the figure, we see that, for small $\Delta T$, such that the first term of Eq. (72) prevailed, the current-current correlations are proportional to $T^{2}$. On the contrary when $\Delta 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 $\Delta T$ as the temperature increases, and appears to depend on the square of $T$ as expected in the classical limit.

In Fig. 3 (b) we study the dependence of current-current correlations on the temperature difference $\Delta T$. At large temperature difference, i.e. beyond the linear limit considered in (71), the correlations depend linearly on $\Delta 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 $\Delta T$, where (70) is well approximated by Eq. (71), indeed we find the saturation on the correlations predicted by the second term of Eq. (70). This implies that even as $\Delta T\rightarrow 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 [48], and no entanglement could be detected in this regime.

5.2 Driven case: Heat rectification

Heat rectification is quantified by the heat rectification, Eq. (76). We use the particular form of driving

[TABLE]

i.e. we drive the frequency of the first oscillator. The parameters that we use are given in Table 2. The temperatures we choose are upper bounded by $\{T_{L},T_{R}\}<10^{3}nK$.

In Fig. 4 we depict the heat rectification coefficient $R$ as a function of $\Omega_{d}$, the driving frequency. We find that it shows to maxima at for two values of the driving frequency $\Omega_{d}\in\{\Omega,3\Omega\}$. Note that, as shown in [43], heat rectification should be maximum at the following frequencies

[TABLE]

Here $\{i,j\}\in\left\{\rm{L,R}\right\}$ with $i\neq j$ and $\nu_{i}$’s are the normal modes of the coupled impurities

[TABLE]

with $\Delta=\Omega_{\rm{R}}^{2}-\Omega_{\rm{L}}^{2}$. In our case, Eq. (79) indeed suggests that rectification should be maximum at $\Omega_{d}\in\{\Omega,3\Omega\}$ which explains the results in Fig. 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 [43].

6 Conclusions

In this work, we studied in detail the heat transport control and heat current rectification among two Bose Einstein Condensates (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. 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 spectral density 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 [43]. 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 [48], 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 spectral density, which opens a new window for further 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 [70, 71] 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.

Acknowledgements

We thank fruitful discussions with Andreu Riera-Campeny. We acknowledge the Spanish Ministry MINECO (National Plan 15 Grant: FISICATEAMO No. FIS2016-79508-P, SEVERO OCHOA No. SEV-2015-0522), the Ministry of Education of Spain (FPI Grant BES-2015-071803), European Social Fund, Fundació Cellex, Generalitat de Catalunya (AGAUR Grant No. 2017 SGR 1341 and CERCA/Program), ERC AdG OSYRIS and NOQIA, EU FETPRO QUIC, and the National Science Centre, Poland-Symfonia Grant No. 2016/20/W/ST4/00314. MM acknowledges support from the Spanish MINECO (QIBEQI FIS2016-80773-P and Severo Ochoa SEV-2015-0522), Fundació Privada Cellex, and the Generalitat de Catalunya (CERCA Program and SGR1381).

Appendix A The uncertainty relation

Here, we present how the covariance matrix of our system is obtained, (both for the static and driven cases) which we used in order to verify that our system fulfils the uncertainty principle. This allows us to guarantee that the system under study, for the parameters considered and according with the assumptions that we have imposed, is a physical system. To do so, one needs to obtain the covariance matrix of the system,

[TABLE]

where $\mathbf{P}=\left(p_{1}\left(t\right),p_{2}\left(t\right)\right)^{T}$ and

[TABLE]

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\rightarrow\infty$ which is equivalent to consider that $\omega\ll\Lambda_{\rm{L}},\Lambda_{\rm{R}}$. The $4\times 4$ matrix in (81) is constructed from the vector $\mathbf{Y}\left(t\right)=\left(x_{1}\left(t\right),x_{2}\left(t\right),p_{1}\left(t\right),p_{2}\left(t\right)\right)^{T}$ as the product $\frac{1}{2}\left\langle\left\{\mathbf{Y}^{T}\left(t\right)\mathbf{Y}\left(t^{\prime}\right),\left(\mathbf{Y}^{T}\left(t\right)\mathbf{Y}\left(t^{\prime}\right)\right)^{T}\right\}\right\rangle_{\rho_{\rm{B}}}$.

The uncertainty relation will then be expressed as a condition on the symplectic transform of the covariance matrix $\widetilde{\mathbf{C}}\left(0\right)$, where the latter is obtained as $\widetilde{\mathbf{C}}\left(0\right)=i\mathbf{W}\cdot\mathbf{C}\left(0\right)$ with $\mathbf{W}$ the symplectic matrix \mathbf{W}=\left(\begin{array}[]{cccc}0&0&1&0\\ 0&0&0&1\\ -1&0&0&0\\ 0&-1&0&0\end{array}\right). The uncertainty relation then simply reads as

[TABLE]

where $\nu_{-}$ is the minimum standard eigenvalue of $\widetilde{\mathbf{C}}\left(0\right)$ ($\hbar$ 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

[TABLE]

where

[TABLE]

with \mathbf{T}=\left(\begin{array}[]{cc}T_{\rm{L}}&0\\ 0&T_{\rm{R}}\end{array}\right) being the temperatures of each bath.

Driven case.

In the driven case, since the solution of the equations of motion is periodic, so are the elements of the covariance matrix as well. In particular, in terms of the Fourier components of the expansion of the periodic phonon Green function, the covariance matrix elements read as

[TABLE]

where

[TABLE]