Quantum duets working as autonomous thermal motors

TL;DR

This paper investigates a quantum thermal machine composed of two particles on sinusoidal tracks, revealing how broken symmetry induces directed motion and enabling work extraction without cyclic protocols.

Contribution

It provides an exact expression for the velocity of the quantum motor in the small undulation limit and explores the effects of external forces and particle interactions.

Findings

Non-zero velocity arises from broken spatial symmetry.

Work can be extracted without cyclic protocols.

Results validate quantum molecular dynamics algorithms in non-equilibrium regimes.

Abstract

We study the dynamic properties of a thermal autonomous machine made up of two quantum Brownian particles, each of which is in contact with an environment at different temperature and moves on a periodic sinusoidal track. When such tracks are shifted, the center of mass of the system exhibits a non-vanishing velocity, for which we provide an exact expression in the limit of small track undulations. We discuss the role of the broken spatial symmetry in the emergence of directed motion in thermal machines. We then consider the case in which external deterministic forces are applied to the system, and characterize its steady state velocity. If the applied external force opposes the system motion, work can be extracted from such a steady state thermal machine, without any external cyclic protocol. When the two particles are not interacting, our results reduce to those of refs. [1,2] for a single particle moving in a periodic tilted potential. We finally use our results for the motor velocity to check the validity of the quantum molecular dynamics algorithm in the non--linear, non--equilibrium regime.

Equations

$$H_0(\{P_i\},\{Q_i\})$$

$$V(Q_1,Q_2)$$

$$V_1(Q_1)$$

$$V_2(Q_2)$$

$$U(Q_1-Q_2)=\frac{k}{2}(Q_1-Q_2{)}^2,$$

$$H=H_0+H_1+H_2$$

$$H_i=\sum _k\left(\frac{{\mathcal{P}}_{k,i}^2}{2m_{k,i}}+\frac{1}{2}{m}_{k,i}{\omega }_{k,i}^2({\mathcal{X}}_{k,i}-Q_i{)}^2\right)$$

$$N_i(\omega )=2\pi \sum _k{m}_{k,i}{\omega }_{k,i}^2\delta (\omega -{\omega }_{k,i}),$$

$${\eta }_i(t)=\theta (t)\int \frac{d\omega }{2\pi }{N}_i(\omega )\cos \omega t.$$

$$m{\ddot{Q}}_i=-{\partial }_{i}V-{\int }_{t_0}^t{\eta }_i(t-t^{\prime }){\dot{Q}}_i(t^{\prime })\mathrm{d}{t}^{\prime }+{\xi }_i(t),$$

$$m{\ddot{Q}}_1$$

$$m{\ddot{Q}}_2$$

$$[{\xi }_i(t),{\xi }_j(t^{\prime })]={\delta }_{ij}\int \frac{d\omega }{2\pi }{N}_i(\omega )\hslash \omega {\mathrm{e}}^{-\text{mathbbm}i\omega (t-t^{\prime })},$$

$$⟨\{{\xi }_i(t),{\xi }_j(t^{\prime })\}⟩={\delta }_{ij}\int \frac{d\omega }{2\pi }{N}_i(\omega )\hslash \omega {\mathrm{e}}^{-\text{mathbbm}i\omega (t-t^{\prime })}\coth \frac{\hslash \omega }{2T_i}.$$

$$⟨{\xi }_i(t){\xi }_j(t^{\prime })⟩={\delta }_{ij}\int \frac{d\omega }{2\pi }{N}_i(\omega )\frac{\hslash \omega }{2}{\mathrm{e}}^{-\text{mathbbm}i\omega (t-t^{\prime })}\left(1+\coth \frac{\hslash \omega }{2T_i}\right).$$

$$⟨{\tilde{\xi }}_i(\omega ){\tilde{\xi }}_j({\omega }^{\prime })⟩={\delta }_{ij}2\pi \delta (\omega +{\omega }^{\prime })N_i(\omega ){\tilde{F}}_i(\omega ),$$

$${\tilde{F}}_i(\omega )=\frac{\hslash \omega }{2}\left(1+\coth \frac{\hslash \omega }{2T_i}\right).$$

$$N(\omega )=2{\eta }_{0}f(\omega ),$$

$$q_i=q_i^{(0)}+V_0{q}_i^{(1)}+V_0^2{q}_i^{(2)}+O(V_0^3).$$

$$m{\ddot{q}}_1^{(0)}$$

$$m{\ddot{q}}_2^{(0)}$$

$$x(t)=q_1+q_2,y(t)=q_1-q_2,$$

$$m{\ddot{x}}^{(0)}$$

$$m{\ddot{y}}_2^{(0)}$$

$$x^{(0)}(t)$$

$$y^{(0)}(t)$$

$${\tilde{G}}_x(\omega )$$

$${\tilde{G}}_y(\omega )$$

$$q_1^{(0)}$$

$$q_2^{(0)}$$

$${\tilde{q}}_1^{(0)}(\omega )$$

$${\tilde{q}}_2^{(0)}(\omega )$$

$$\displaystyle\big{[}{\tilde{q}_{1}^{(0)}(\omega),\tilde{q}_{1}^{(0)}(\omega^{\prime})}\big{]}=\big{[}{\tilde{q}_{2}^{(0)}(\omega),\tilde{q}_{2}^{(0)}(\omega^{\prime})}\big{]}$$

$$=b^2\pi N(\omega )\hslash \omega \delta (\omega +{\omega }^{\prime })\times \left(|{\tilde{G}}_x(\omega ){|}^2+|{\tilde{G}}_y(\omega ){|}^2\right)$$

$$\displaystyle\big{[}{\tilde{q}_{1}^{(0)}(\omega),\tilde{q}_{2}^{(0)}(\omega^{\prime})}\big{]}=$$

$$=b^2\pi N(\omega )\hslash \omega \delta (\omega +{\omega }^{\prime })\times \left(|{\tilde{G}}_x(\omega ){|}^2-|{\tilde{G}}_y(\omega ){|}^2\right)$$

$$\displaystyle\big{[}{q^{(0)}_{1}(t),q^{(0)}_{1}(t^{\prime})}\big{]}=\big{[}{q^{(0)}_{2}(t),q^{(0)}_{2}(t^{\prime})}\big{]}=$$

Full text

Emerging technologies have been used to engineer and realise quantum devices which are the nanoscopic counterparts of heat engines Roßnagel et al. (2016) or thermoelectric transducers Thierschmann et al. (2015). Similarly to their macroscopic analogues, quantum engines are open systems that exchange energy with the surrounding environments in form of work and heat, and the study of their dynamic and thermodynamic properties has attracted considerable interest Devoret et al. (2014); Kosloff (2013) in an effort to extend the concepts of the classical thermodynamics to the quantum regime Vinjanampathy and Anders (2016); Ali .

The properties of quantum reciprocating motors have been extensively investigated. Such motors are characterized by a time dependent Hamiltonian, and interact cyclically with baths at different temperatures, so as to perform Carnot, Otto or Stirling cycles Kosloff (2013); Gelbwaser-Klimovsky et al. (2013a, b); Insinga et al. (2016); Ali from which work can be extracted.

Autonomous motors are based on a fundamentally different design: these machines can operate in steady state conditions without any external agent changing their Hamiltonian, or their heath bath. Within this framework, a few works have recently appeared where the dynamic and thermodynamic properties of quantum rotors have been studied and where the sole driving force is a thermal gradient Mari et al. (2015); Roulet et al. (2017); Seah et al. (2018); Fogedby and Imparato (2018); Roulet et al. (2018); Hovhannisyan and Imparato (2019). In such models the heat current between different baths is converted into rotational motion.

The properties of autonomous motors exhibiting directed particle transport on rectilinear tracks have not received much attention in the past, with the noticeable exception of refs. Fisher and Zwerger (1985); Aslangul, C. et al. (1987), whose authors studied the steady state velocity of a quantum Brownian particle in a periodic potential under the effect of a constant external force that explicitly breaks the spatial symmetry. More recently ref. Bissbort et al. (2017) studied the conversion of energy current flowing between different baths into motion of a particle on a ring.

In this paper we put forward a rigorous approach to investigate the dynamic properties of a quantum autonomous thermal motor fuelled by a temperature gradient, and by spatial broken symmetry. Specifically, we propose a minimal model of quantum motor, which is based on the directed transport emerging in systems with both broken spatial symmetry and thermal equilibrium. This model motor, is the quantum counterpart of the classical system considered in Fogedby and Imparato (2017), and consists of two particles sitting in two periodic shifted potentials, interacting through a third potential, and kept at different temperatures.

We use the Caldeira-Legget model Caldeira and Leggett (1983a, b); Ford et al. (1988) for the two heat baths, which allow us to write the dynamic equations for the particles’ coordinates in the Heisenberg representation. In the limit of small undulations in the periodic potentials, the problem can be solved analytically, and the system steady state velocity can be obtained by studying the motion of a free quantum Brownian particle, and of a Brownian quantum oscillator. This choice for the baths does not require any assumption on the strength of the system-bath interaction, and allows us to derive our results for baths with an arbitrary distribution of internal modes. This is relevant in the context of systems interacting with non-Markovian baths, that has attracted considerable attention in the physics community working on dynamics of open quantum systems de Vega and Alonso (2017); Naz . Indeed, while it is commonly believed that non-Markovian behaviour emerges in systems which interact strongly with baths with a structured spectral density de Vega and Alonso (2017), it has been recently shown that non-Markovian behaviour can emerge in systems as simple as a single harmonic oscillator in contact with a Ohmic bath, provided that the bath and the system are initially correlated Strasberg and Esposito (2018). Furthermore, the study of quantum heat machines in the strong coupling regime has received considerable attention in recent years Gelbwaser-Klimovsky and Aspuru-Guzik (2015); Uzdin et al. (2016); Newman et al. (2017).

We discuss the case in which deterministic external forces are applied to the autonomous motor. The case of constant external force is particularly interesting, as a force counteracting the system center-of-mass motion can extract work from the motor while it operates in steady state conditions without any cyclic manipulation protocol, at variance with the setups characterizing the reciprocating engines, as described, e.g., in Ali .

We finally use the results obtained for the motor steady state velocity to check the validity of the quantum molecular dynamics (QMD) algorithm in the non-linear, out-of-equilibrium, multi-bath regime.

I The model

In this section we introduce the autonomous motor model as the quantum conterpart of the classical model introduced in Fogedby and Imparato (2017). The system Hamiltonian reads

[TABLE]

with $V_{1}$ and $V_{2}$ two periodic potentials and $U$ an interaction potential. In the following we will take

[TABLE]

As discussed in Fogedby and Imparato (2017), when one takes the interaction potential to be $U(x_{1}-x_{2})=-k\cos[b(x_{1}-x_{2})]$, eq. (1) becomes the Hamiltonian for the $xy$-model in an external field, describing the elastic free energy in ferromagnetic or liquid-crystal systems Chaikin and Lubensky (1995). Throughout this paper we will use the expression

[TABLE]

as it will allow us to obtain closed results for the dynamics when $V_{0}=0$, as detailed below. The Hamiltonian eq. (1) breaks the spatial symmetry in the sense that when $\varphi\neq 0,\pi$ there is no unilateral translation $\Delta$ such that $V(-Q_{1},-Q_{2})=V(Q_{1}+\Delta,Q_{2})$ $\forall Q_{1},Q_{2}$, (or similarly for $Q_{2}$, $V(-Q_{1},-Q_{2})=V(Q_{1},Q_{2}+\Delta)$). This broken symmetry is the key ingredient for the system to exhibit directed transport, as discussed below.

We model the thermal baths as ensembles of coupled harmonic oscillators Caldeira and Leggett (1983a, b); Ford et al. (1988), and thus the total Hamiltonian (system+bahts) reads

[TABLE]

with

[TABLE]

and $i=1,\,2$, and where $\mathcal{X}_{k,i}$ and $\mathcal{P}_{k,i}$ are the bath variables. While every harmonic oscillator in each bath is characterized by the frequency $\omega_{k,i}$ and the mass $m_{k,i}$, the details of the coupling between the baths and the system are embedded in these two quantities Ford et al. (1988). By introducing the density of states

[TABLE]

we can treat the spectrum for each bath as continuous. In this case, the memory function $\eta_{i}(t)$ is defined as

[TABLE]

Following the procedure described in Caldeira and Leggett (1983a, b), one obtains the quantum Langevin equations

[TABLE]

which are the dynamical equations for the operators $Q_{i}$ in the Heisenberg picture. The interaction with the baths are now embodied by the quantum noise operators $\xi_{i}(t)$. For the specific choice of the potential (2)–(4), one obtains

[TABLE]

The commutators and anti-commutators of the quantum noise operators are given by

[TABLE]

with two time correlation

[TABLE]

which in Fourier space reads

[TABLE]

and where we have introduced

[TABLE]

In the previous equations we have taken $k_{B}=1$, a simplification that we will keep in the following.

In order to simplify the notation in the following we will introduce the dimensionless space coordinates $q_{i}=bQ_{i}$.

Furthermore we will assume that the density of states for the baths are equal $N_{1}(\omega)=N_{2}(\omega)=N(\omega)$, and thus $\eta_{1}(t)=\eta_{2}(t)=\eta(t)$. In principle $N(\omega)$ is non-zero for $\omega\geq 0$, but in order to make the integrals in the following sections symmetric we can take $N(\omega)=N(-\omega)$, with $\omega$ ranging over the entire real axis. The density of states takes the form

[TABLE]

where typical choices for the cutoff function $f(\omega)$ are $f(\omega)=1$ (Ohmic bath), or with a soft cutoff $f(\omega)=\Lambda^{2}/(\omega^{2}+\Lambda^{2})$ or $f(\omega)=\exp(-|\omega|/\Lambda)$.

II perturbation expansion

A simple symmetry argument indicates that, if the system exhibits a non-vanishing velocity, then it is invariant under sign inversion of the undulation amplitude $V_{0}\to-V_{0}$. Indeed inspection of the Hamiltonian (1) suggests that changing sign to $V_{0}$ corresponds to a translation of the coordinates $Q_{i}\to Q_{i}\pm\pi/b$ which does not change the phase shift between the potentials $V_{1}$ and $V_{2}$ in eq. (2). In other words, the Heisenberg equations (11)–(12) would be the same for the shifted coordinates under sign inversion of $V_{0}$. Thus we conclude that any systematic non-vanishing velocity must be even in $V_{0}$.

Therefore, following Aslangul, C. et al. (1987) we expand the dimensionless coordinates $q_{1},\,q_{2}$ up to second order in the particle potential amplitude $V_{0}$,

[TABLE]

To zeroth order the quantum Langevin equations (10) read

[TABLE]

It is furthermore convenient to introduce the coordinates

[TABLE]

and their corresponding power series. To the zero-th order, the coordinate $x^{(0)}(t)$ describes the motion of a free Brownian particle, while $y^{(0)}(t)$ describes a particle in a harmonic potential with strength $2k$, and their quantum Langevin equations read

[TABLE]

Since we are interested in the long-time behaviour of the system, we will set $t_{0}=-\infty$ so as the one time averages will be time independent, and the two-time correlation functions will be time translationally invariant. To zero-th order we have thus that the Green’s function solutions of eq. (10) in terms of the coordinate $x$ and $y$ read

[TABLE]

where the Fourier transform of $G_{x}(\tau)$ and $G_{x}(\tau)$ read

[TABLE]

Thus one finally obtains the zeroth-order solutions

[TABLE]

It is straightforward to show that the commutators \big{[}{q^{(0)}_{i}(t_{1}),q^{(0)}_{j}(t_{2})}\big{]} are complex-valued functions of the time lapse $t_{1}-t_{2}$. We have indeed

[TABLE]

By taking into account that $\tilde{G}_{x,y}(-\omega)=\tilde{G}^{*}_{x,y}(\omega)$, and by using eq. (13), we can then calculate

[TABLE]

obtaining thus

[TABLE]

In order to shorten the notation in the following we set

[TABLE]

and

[TABLE]

which are $c$-numbers. We notice for later use that $a_{12}(t-t^{\prime})=a_{21}(t-t^{\prime})$, which can be obtained from eq. (34) by taking into account that $N(\omega),\,|\tilde{G}_{x}(\omega)|^{2}$ and $|\tilde{G}_{y}(\omega)|^{2}$ are all even functions of $\omega$.

Similarly, from eqs. (16) and (29)-(30) we can calculate the correlations $\langle{q^{(0)}_{i}(t)q^{(0)}_{j}(t^{\prime})}\rangle$ (see appendix A) and from them one can obtain the quantities

[TABLE]

Similarly one finds

[TABLE]

II.1 First order

In order to simplify the calculations in the following we will take the phase appearing in eq. (4) to be $\varphi=\pi/2$. We will then generalise our results to the case of arbitrary $\varphi$. To first order in $V_{0}$, eqs. (10) become

[TABLE]

with solutions

[TABLE]

where

[TABLE]

We can now check that the first order term in $V_{0}$ of the velocity vanishes. One possible way to proceed is to differentiate eq. (LABEL:X1), and to perform a statistical average over the quantum noise operators $\xi_{i}(t)$. One thus obtains

[TABLE]

The calculation of the averages can be performed by expressing the trigonometric functions in their exponential forms. Thus one has to evaluate terms of the form $\langle{\exp(\pm\mathbbm{i}q_{i}(t))}\rangle$. This can be done by noticing that the position operators $q^{(0)}_{i}(t)$ are linear combinations of the noise operators $\xi_{i}$, eqs. (23)–(24), which in turn are normally distributed with zero mean. If $X$ is a normally distributed variable, with zero mean, one easily finds that $\langle{\exp(X)}\rangle=\exp(\langle{X^{2}}\rangle/2$). By applying this equality to the trigonometric functions on the right hand side (rhs) of eq. (45), one finds $\langle{\sin q^{(0)}_{1}(t)}\rangle=0$, and given that $\langle{(q^{(0)}_{2}(t))^{2}}\rangle$ is independent of the time, (see appendix A), also $\langle{\cos q^{(0)}_{2}(t)}\rangle$ it time independent. Thus we conclude that the rhs of eq. (45) vanishes.

II.2 Second order

To second order in $V_{0}$, eqs. (11)–(12) become

[TABLE]

where the notation \big{\{}{\sin(q^{(0)}_{1}+V_{0}q^{(1)}_{1})}\big{\}}^{(1)} and \big{\{}{\cos(q^{(0)}_{2}+V_{0}q^{(1)}_{2})}\big{\}}^{(1)} indicates the first order term in the expansion. However the expansion of the trigonometric functions must be taken carefully, given that the operators $q^{(0)}_{i}(t)$ and $q^{(1)}_{i}(t)$ do not commute.

Assuming that the system reaches a constant average center-of-mass velocity in the steady state (which can well be vanishing) is equivalent to take $\langle{\ddot{x}^{(2)}}\rangle=0$, where the average is taken over the quantum noise operators $\xi_{i}$. Thus summing up eqs.(46)-(47), and averaging over the quantum noise, we obtain the following equation for $\dot{x}^{(2)}(t)$

[TABLE]

The next step is to expand the terms on the right hand side of the last equation up to the first order in $V_{0}$. The derivation is long but quite straightforward, and is presented in appendix B.

For the second order contribution to the steady velocity one thus obtains

[TABLE]

where

[TABLE]

Thus, up to second order in $V_{0}$ the velocity of the center of mass reads

[TABLE]

We now consider the case where the potentials $V_{i}$ have an arbitrary phase shift $\varphi$: $V_{2}(Q_{2})=-V_{0}\cos(bQ_{2}+\varphi)$. By introducing the shifted variable $\bar{q}^{(0)}_{2}=q^{(0)}_{2}+\varphi$ and retracing the previous steps, where we have taken the specific value $\varphi=\pi/2$, one obtains for the steady state velocity up to the second order in $V_{0}$

[TABLE]

In this section we have derived the result (49) by assuming that $\langle{\ddot{x}^{(2)}}\rangle=0$. While it is reasonable to assume that in the steady state the average center-of-mass velocity is constant, in appendix D we provide an exact proof of the result $\langle{\ddot{x}^{(2)}}\rangle=0$.

III Second order velocity term: analysis

We now analyze the results for the constant second order velocity, eqs. (49)–(52). Plots of the system center of mass velocity as a function of the interaction strength $k$, of the wavenumber $b$, and of the temperature scale are shown in fig. 1 for the soft cutoff function $f(\omega)=\Lambda^{2}/(\omega^{2}+\Lambda^{2})$. The physical parameter space considered in Fig. 1 has been chosen in order to analyze a potential experimental realization of the model where the two particles are represented by two sideband laser-cooled atomic ions. Confinement of the ions by a common harmonic trap potential will lead to an effective harmonic binding force between them Grimm et al. (2000), while spatially periodical potentials along the axis defined by the two ions can be realized through off-resonant electrical dipole forces induced by standing wave light fields Linnet et al. (2012); Enderlein et al. (2012); Bylinskii et al. (2015); Lauprêtre et al. (2019). By either choosing two identical ion species initialized in different internal states or two different ion isotopes, different polarization states and/or longitudinal modes of an Fabry-Perot cavity could enable particle dependent periodical potentials with a constant, but tunable phase relation between them Linnet et al. (2014).

When does the right hand side term of equation (50), and thus $\bar{v}^{(2)}$ vanish? From eq. (17) we see that $\tilde{F}_{1}(\omega)=\tilde{F}_{2}(\omega)$ when $T_{1}=T_{2}$. By inspecting eqs. (37)–(38) we find that $c_{12}(t-t^{\prime})=c_{21}(t-t^{\prime})$ for equal temperatures. Therefore we conclude that the velocity $\bar{v}^{(2)}$ vanishes at thermal equilibrium, as expected: it is the heat current flowing between the two baths that sustains a non-vanishing velocity.

Furthermore the the center-of-mass velocity vanishes in the trivial case of vanishing undulation amplitude $V_{0}=0$, i.e., the motor requires gear racks in order to work, which are represented by the two periodic potentials $V_{1}$ and $V_{2}$.

Most importantly, we notice that the second order term of the velocity vanishes for a phase shift between the potentials $\varphi=l\pi$ with $l$ an integer number. This is consistent with the results for the classical counterpart of the present model found in Fogedby and Imparato (2017), where it was numerically shown that the velocity vanishes for $\varphi=l\pi$ independently of the order of $V_{0}$. In this case the system does not break the spatial symmetry discussed in section I. Such a broken symmetry has been shown to be a prerequisite in order for directed motion in classical and quantum duets to arise Fogedby and Imparato (2017); Hovhannisyan and Imparato (2019), as well as in many-body systems Suñé and Imparato (2019) . We recall that the required broken symmetry reads as follows: there is no translation distance $\Delta$ such that $V(-Q_{1},-Q_{2})=V(Q_{1}+\Delta,Q_{2})$. Such a broken symmetry is the 2D counterpart of the broken symmetry discussed in Reimann (2002) for 1D non–autonomous Brownian motors. In 1D this broken spatial symmetry amounts to require that the potential is, e.g., saw-tooth shaped.

One should also expect that the velocity vanishes in the limit $k\to 0$ (decoupled system) and in the limit of tight coupling $k\to\infty$. In the latter case the system behaves as a single particle in contact with two environments at different temperatures, for which the fluctuations of the relative coordinate are suppressed $\langle{y}\rangle\to 0$, leading to $\tilde{G}_{y}(\omega)\to 0$, which in turn implies that the asymmetric term in (37) vanishes. In simpler terms, in the limit $k\to\infty$ the sinusoidal tracks, which are responsible for the broken spatial symmetry and thus for the motor effect, become negligible with respect to the interaction potential (5), as $V_{0}/k\to 0$.

Therefore one should expect an optimal coupling strength between these two regimes, as confirmed by inspection of fig. 1-(a), where we plot the second order steady state velocity $\bar{v}^{(2)}$ as a function of the particle-particle interaction strength for different values of the cutoff frequency. Furthermore from eq. (52) we find that the velocity, up to the second order in $V_{0}$, vanishes in the limit of $b\to 0$ (infinite period). Inspection of eqs. (11)–(12) suggests that the center of mass velocity must also vanish in the limit $b\to\infty$ as the force exerted by the particles’ potentials will prevail on the bath forces $\xi_{i}$ which drive the directed motion. Thus one must expect an optimal value of $b$ for which the velocity achieves a maximum. This is confirmed by inspection of fig. 1-(b) where we plot $\bar{v}^{(2)}$ as a function of the wavenumber $b$, for different values of the cutoff frequency.

The velocity $\bar{v}^{(2)}$ is expected to vanish in both the large and the low temperature regime. The limit of large $T_{1}$ and $T_{2}$ corresponds trivially to the limit $V_{0}\to 0$, where the gears represented by the periodic tracks are flattened. In the limit of low temperatures it is the heat currents that fuels the motor that vanishes. These considerations are confirmed by inspection of fig. (1)-(c), where the second order steady state velocity is plotted as a function of the temperature scale.

The friction coefficient appearing in the expression of the density of states $N(\omega)$, eq. (18), plays the role of an effective coupling strength with the baths’ oscillators, see also eq. (8). It is thus interesting to evaluate the velocity $\bar{v}^{(2)}$ as a function of such a parameter. One expects that for large $\eta_{0}$, the motor is so strongly coupled with the baths’ degrees of freedom, that its velocity should vanish. Similarly for vanishing $\eta_{0}$, the motor velocity should vanish, as the propelling effect of the baths would also vanish. However, one should keep in mind that for very small values of $\eta_{0}$ the integrand in eq. (50) becomes highly oscillating and its numerical evaluation becomes unfeasible. We therefore consider here the lower bond $\eta_{0}/m>1$ Hz. The results for the velocity as function of $\eta_{0}$ are shown in fig. (1)-(d), and we find that in the range $\eta_{0}/m>1$ Hz the velocity decreases with different rates for increasing $\eta_{0}$, depending on the cutoff frequency $\Lambda$.

Inspections of the four panels in fig. 1 indicates that in general the velocity increases for larger cutoff $\Lambda$. This can be easily understood by noticing that it is the modes in the baths with density of states $N(\omega)$ as given by eq. (18) that propel the motor. The larger the cutoff $\Lambda$ the more modes are present in the bath, with higher frequency, and thus larger average thermal energy. This is no longer true in the limit of larger temperatures, see fig. 1-(c). This can be understood by considering that, as discussed above, the limit of large temperatures correspond to the limit of $V_{0}\to 0$. So having fewer modes in the baths (with lower average energy) is beneficial for the motor as the flattening of the tracks is reduced with respect to the case of large $\Lambda$. It is important to notice that the results obtained in this section for $\bar{v}^{(2)}$ are exact up to the second order in $V_{0}$. Thus $V_{0}$ has to be small compared to the other energy scales in the system: $k_{B}T_{i}$, $k/b^{2}$ and $\hbar\sqrt{k/m}$. In this regard, if one considers the limit where one of the two temperatures vanishes, e.g. $T_{1}\to 0$, while the second is finite $T_{2}>0$, the potential amplitude must also go to zero $V_{0}\to 0$, and in this limit the motor effect vanishes, as discussed above.

The above results for the velocity, and in particular eqs. (50)–(51), hold true also in the classical limit of large temperatures ($\hbar\to 0$), provided that in the expression of the correlation $\langle(q^{(0)}_{1}(\tau)-q^{(0)}_{2}(t))^{2}\rangle$ eq. (37), one takes $\tilde{F}_{i}(\omega)\to T_{i}$ for $\hbar\to 0$, and by noticing that from eq. (35) one finds

[TABLE]

IV External forces

We now consider the case where two external forces $F_{i}(t)$ are applied on the two particles. In terms of the rescaled coordinates the quantum Langevin equations (10) become

[TABLE]

To zeroth order the solutions of eqs. (54)–(55) consist now of two contributions, one arising from the random forces $\xi_{i}(t)$ and one from the systematic forces $F_{i}(t)$, and read

[TABLE]

[TABLE]

with $q^{(0)}_{1,2}=(x^{(0)}\pm y^{(0)})/2$. By comparing the equations eqs. (56)-(57) with eqs. (23)-(24) one finds that, as expected, the random part of the solutions $x^{(0)}_{R}(t)$, $y^{(0)}_{R}(t)$ are identical to the solutions $x^{(0)}(t)$, $y^{(0)}(t)$ discussed in sec. II. Similarly, to the first order, one finds

[TABLE]

which are identical to the first order solutions (LABEL:X1)–(LABEL:Y1) in sec. II.

Finally, to second order in $V_{0}$ one obtains

[TABLE]

which are identical to the second order quantum Langevin equations (46)–(47) in absence of external force. Thus, up to this point, the only difference when one applies the external force $F_{i}(t)$ is in the zeroth order solutions (56)–(57). In particular the calculation of the first order expansion of the rhs of equations (60)–(61) can be carried out as in appendix B obtaining identical expressions to eqs. (75)–(76). As in section  II and in appendix B we need to calculate the average of the rhs of eqs. (60)–(61) over the bath variables $\xi_{i}(t)$. However, when averaging the rhs of eqs. (75)–(76), we have to take into account that the averages $\langle\sin q^{(0)}_{i}(\tau)\cos q^{(0)}_{i}(t)\rangle$ do not vanish, as the zeroth order solution contains now a contribution from the deterministic forces $F_{i}(t)$, see eqs. (56)–(57). With these considerations in mind we can thus calculate the average of the rhs of (60)–(61) and obtain

[TABLE]

IV.1 Constant forces

We now consider the case in which the deterministic forces are constant $F_{i}(t)=\bar{F}_{i}$. The system will thus reach a steady state with a constant average velocity. By summing eqs. (62)–(63), and exploiting the fact that the two time commutators for $q^{(0)}_{1}$ and $q^{(0)}_{2}$ are identical (eq. (33)), we obtain for the center of mass

[TABLE]

where

[TABLE]

with the systematic part of the zeroth order solution given by

[TABLE]

and where the expressions for the correlations $c_{11}(t)$ and $c_{22}(t)$ are given by eq. (73) in appendix A.

Thus, up to the second order in $V_{0}$ the velocity of the center of mass reads

[TABLE]

If the forces $\bar{F}_{1,2}$ have opposite sign with respect to the system velocity, and are not large enough in modulus to invert the direction of the motion, one can thus extract work from the thermal machine by doing work against such external forces. If one takes the case $\bar{F}_{1}=\bar{F}_{2}=\bar{F}$, the output work rate reads

[TABLE]

When one consider the non interacting case $k=0$, the two particles move independently under the effect of the two forces $\bar{F}_{1}$ and $\bar{F}_{2}$. In this limit the Green’s function $G_{y}(t)$ becomes equal to $G_{x}(t)$ (see eqs. (25) and (26)) and the second term in eq. (65) vanishes. This is the term responsible for the thermal propulsion even in absence of external forces, see eqs. (50)–(51). On the other hand the first term in eq. (65) is non zero when $\bar{F}_{1}$ or $\bar{F}_{2}\neq 0$. Thus for $k=0$, eqs. (65)–(67) reduce to the results of Fisher and Zwerger (1985); Aslangul, C. et al. (1987) for the case of a single quantum Brownian particle in a sinusoidal potential under the effect of a constant external force.

V Simulations

In the Quantum Molecular Dynamics (QMD) algorithm, as introduced in Dammak et al. (2009) for systems at equilibrium with a single bath, a quantum thermal bath is replaced by a classical bath that accounts for quantum statistics in the framework of a standard MD algorithm. Specifically, in the QMD algorithm the Heisenberg equations of motion for the quantum operators $Q_{i}$ and $P_{i}$ are replaced by a classical Langevin equation, of the type (10), where the power spectral density of the stochastic noise is given by the quantum mechanical fluctuation-dissipation relation (17). Thus the QMD algorithm neglects one basic quantum feature, namely the noncommuting character of the system variables, while it retains the power spectral density of the bath variables.

In order to avoid any possible misinterpretation, we remark once more that the formalism in sections I–IV is purely quantum, as eq. (10) is a dynamical equation for the operators $Q_{i}$ in the Heisenberg picture. Although the Ehrenfest theorem prescribes the equivalence of classical and quantum dynamic equations only for harmonic potentials, the QMD algorithm has provided accurate results for different types of systems with various degrees of anharmonicity Dammak et al. (2012); Calvo et al. (2012a, b); Qi and Reed (2012); Bronstein et al. (2014). For example the validity of the approach has been tested by reproducing several equilibrium experimental data at low temperatures in a regime where quantum statistical effects cannot be neglected, for MgO crystal or nonsuperfluid liquid ${}^{4}\mathrm{He}$, characterized by anharmonic potentials Dammak et al. (2009).

The QMD has been extended to the non–equilibrium case, with two heat reservoirs at different $T_{i}$ in ref. Fogedby and Imparato (2018), but only for the case of a parabolic potential.

Given that the results contained in section III are exact, up to the second order in $V_{0}$, the model discussed in the present paper represents an excellent test-bed to check whether the QMD can be used to evaluate the dynamic properties of a system that is both non–linear and out-of-equilibrium. We will consider in the following Ohmic baths ($\Lambda\to\infty$), while the numerics have been performed by using dimensionless quantities for the system parameters. Similarly to ref. Dammak et al. (2009), we use the numerical method discussed in Maradudin et al. (1990) to generate the correlated noise with fluctuation-dissipation relation (17).

The results are reported in fig. 2. We notice that the agreement between the second order velocity $\bar{v}^{2}$ and the results obtained from the QMD is quite good, in the regime where the corrugation amplitude $V_{0}$ is smaller that the other energy scale (i.e. the thermal energy $k_{B}T$ and the energy scale associated with the harmonic interaction $\hbar\sqrt{k/m}$). For completeness we also plot the steady state velocity as obtained from the classical molecular dynamics (MD) algorithm.

As discussed in section IV.1, when the interaction strength $k$ is set to zero, the model described by eqs. (54)–(55) reduces to two independent particles moving under the effect of the external forces $F_{i}$. This corresponds to the quantum Brownian particle in a tilted sinusoidal potential $U(x)=-V_{0}\cos(bx)-Fx$, discussed in Fisher and Zwerger (1985); Aslangul, C. et al. (1987). In those references the second order velocity of the single particle has been calculated as a function of the applied force. It is thus interesting to compare that exact results with the outcomes of the QMD algorithm in presence of an external constant force. Such a comparison is shown in fig. 3 for two different choices of the parameter set: the agreement between the expected curve and the velocity predicted by the QMD algorithm is quite good.

The agreement between the predicted second order velocity and the QMD algorithm worsens if one, e.g., reduces the potential period, see fig. 4. This can be understood by noticing that the equilibrium QMD algorithm has already been reported to provide approximated results when tunnelling between neighbour wells become predominant Brieuc et al. (2016).

VI Conclusions

In this paper we have introduced and discussed the properties of a steady state quantum motor that can continuously convert heat flows into motion and thus work. As such the motor is different from the reciprocating motors performing thermodynamics cycles with ”moving parts” which play the role of ”pistons” in the classical picture, as discussed, e.g, in Ali .

The present model relies only on the broken spatial symmetry of the underlying potentials and on the temperature difference between the two baths. Had one considered the interaction potential $U(x_{1}-x_{2})=-k\cos[b(x_{1}-x_{2})]$ discussed in section I, the system total potential (2) would be periodic, and so the expected velocity would be the same if one considers periodic tracks (the two potentials $V_{i}(Q)$) with periodic boundary conditions, where the system only can move in the interval $[0,2\pi/b]$, or periodic tracks with open boundary conditions, where the system can move in the range $\left]-\infty,+\infty\right[$. With such a cosine potential our model would be similar to the quantum rotating gears discussed in, e.g., MacKinnon (2002); Liu et al. (2019). The quadratic potential used through this paper, eq. (5), represents the large coupling limit (large $k$) of such a cosine interaction potential. Therefore in this limit, characterized by small fluctuations of the relative coordinate $y$, one should also expect that the steady state velocity is the same for the two different boundary conditions.

Experimentally, there are quite a few examples in literature of systems moving in periodic potentials and subject to temperature gradients. In the context of classical stochastic thermodynamics, experiments where a single particle is constrained on periodic circular tracks have been performed, e.g, in Blickle et al. (2007) where a 3D toroidal laser trap was used to force a colloidal particle to perform circular trajectories along a periodic potential. Brownian systems with 2 degrees of freedom, where the local temperature can be controlled, have been studied, for example, in Ciliberto et al. (2013a, b); Bérut et al. (2016a, b).

The quantum regime considered here should potentially be realizable with two sideband laser-cooled atomic ions commonly confined by a harmonic potential, while individually interacting with each their periodic dipole-induced potential with variable spatial phase relation Linnet et al. (2012); Enderlein et al. (2012); Bylinskii et al. (2015); Lauprêtre et al. (2019); Linnet et al. (2014). The individual temperatures of the two ions can be controlled by addressing each of the ions with specific laser cooling beams. For single ions, temperatures of few $\mu$K have been achieved Poulsen et al. (2012), and even tens of mK deep periodic dipole-induced potentials have been applied to smaller strings of ions Linnet et al. (2012). In this ion-scenario the common harmonic potential will eventual impede the ions’ motion, and instead lead to a stationary offset in the positions of the ions as compared to the thermalized case without the periodical potential. Moving the center of the common trapping potential with a constant velocity corresponding to the relevant one in fig. 1 should, however, not lead to any displacement of the ions in the moving trap frame.

Our results on the QMD provide a solid evidence that such an algorithm can be successfully used to predict dynamical properties, such as particle currents, in quantum, out-of-equilibrium systems in contact with multiple reservoirs. It would be interesting to test the algorithm on other systems whose dynamic properties are known exactly.

Finally, it would be interesting to compare the exact results discussed in this paper with those obtained through a master equation approach, corresponding to a weak-coupling description of the interaction with the baths Breuer and Petruccione (2002). Indeed, it has been found that the local or global master equations may give different dynamics in systems in contact with baths at different temperatures Hovhannisyan and Imparato (2019).

Acknowledgements.

This work was supported by the Danish Council for Independent Research and the Villum Foundation. The numerical results presented in this work were obtained at the Centre for Scientific Computing, Aarhus http://phys.au.dk/forskning/cscaa. We thank F. Barra for a critical reading of the manuscript, and helpful comments. AI is grateful to Hans C. Fogedby for several useful discussions.

Appendix A Two time correlation functions

In this appendix we calculate the two time correlations for the zeroth order operators $\langle{q^{(0)}_{i}(t)q^{(0)}_{j}(t^{\prime})}\rangle$. From eqs. (29)–(30) in the main text we obtain

[TABLE]

We notice that the same time correlation functions $\langle{q^{(0)}_{i}(t)q^{(0)}_{j}(t)}\rangle$ are independent of the time.

One can thus calculate the two-time two-particle correlation function (37) in the main text, and the single particle correlation functions $c_{ii}(\tau)$ appearing in eqs. (62)–(65)

[TABLE]

Appendix B Taylor expansion of the sinusoidal functions

In this appendix we derive the expansion of the trigonometric functions on the rhs of eq. (48) up to the first order in $V_{0}$. By writing the trigonometric functions in their exponential forms, and properly differentiating the exponential operators (see appendix  C) one obtains for the sine

[TABLE]

By setting $\bar{q}^{(0)}_{2}=q^{(0)}_{2}+\pi/2$, one obtains the expansion for the cosine function in eq. (48).

We have now to evaluate the series of nested commutators in (74). One can do this by noticing that the operators $q^{(1)}_{1}$ and $q^{(1)}_{2}$ are linear functionals of $\sin q^{(0)}_{1}(t)$ and $\cos\bar{q}^{0}_{2}(t)$, see eqs. (41)-(LABEL:Y1) and appendix (C). Thus one obtains

[TABLE]

where

[TABLE]

and where the function $a_{ij}(t-t^{\prime})$ has been defined in eq. (35), and is proportional to the commutator of the position operators at different times, which are c-numbers as discussed above. Therefore the functions $\Sigma_{ij}(t-t^{\prime})$ and $\Sigma^{\prime}_{ij}(t-t^{\prime})$ are also complex-valued functions of the time lapse $t-t^{\prime}$. We furthermore notice that the following equality holds

[TABLE]

Our goal is to calculate the average on the right hand side of eq. (48), and we recall that the average has to be taken over the bath variables, represented by the noise operators $\xi_{i}$: such operators appear in the zeroth-order solutions for the position operators $q^{(0)}_{i}(t)$, see eqs. (23)-(24). Thus we have to calculate the average over the bath variables of the products of trigonometric functions appearing on the rhs of eqs. (75)–(76). It is convenient to express those sine and cosine functions in their complex exponential form. Furthermore we notice once more that the commutators \big{[}{q_{i}(t),q_{j}(t^{\prime})}\big{]} are c-numbers, so we can use the Glauber formula

[TABLE]

with $\alpha,\,\alpha^{\prime}=\pm 1$. The zero-th order particle positions $q^{(0)}_{i}(t)$ are linear combinations of the bath variables, which are Gaussian variables with zero average. We can therefore use the following equality

[TABLE]

Thus a straightforward calculation gives

[TABLE]

[TABLE]

[TABLE]

We notice that the ”even” terms $\langle\sin q^{(0)}_{i}(\tau)\cos q^{(0)}_{i}(t)\rangle$ in eqs. (75)–(76) vanish, as can be also checked by a direct calculation.

The variables $q^{(0)}_{1,2}=(x^{(0)}\pm y^{(0)})/2$ depend on the variable $x^{(0)}$ which describes a free Brownian motion. As such the terms $\langle(q^{(0)}_{i}(\tau)+q^{(0)}_{j}(t))^{2}\rangle$ are divergent. Indeed they contain the term $\langle(x^{(0)}(\tau)+x^{(0)}(t))^{2}\rangle$ and by using eqs. (23), (29)–(30) one finds

[TABLE]

and the integrand with the plus sign on the right hand side diverges as $\sim 1/\omega$ for $\omega\to 0$, given that $N(\omega=0)\neq 0$ as discussed in section I. Thus the terms \exp\big{[}{-\langle(q^{(0)}_{i}(\tau)+q^{(0)}_{j}(t))^{2}\rangle/2}\big{]} in eqs. (80)-(83) vanish.

By using the results in eqs. (75)–(76) and (80)–(83), the right hand side of equation (48) can thus be written as

[TABLE]

where

[TABLE]

and where we have used eq. (77).

By noticing that the last term in eq. (84) is time-independent, and by recalling the definition of the memory function (9), one can invert equation (84) and obtain the expression for the long time steady state velocity eq. (49) in the main text.

Appendix C Taylor expansion of eq. (74)

In order to prove eq. (74) in the previous appendix we need to calculate the first order contribution to the expression

[TABLE]

The calculation is more easily performed when the sine is expressed in terms of complex exponentials $\exp(\pm\mathbbm{i}(q^{(0)}(t)+V_{0}q^{(1)}(t))$. Where in general the operators $q^{(0)}(t)$ and $q^{(1)}(t)$ do not commute.

Let $X$ and $Y$ be two non commuting operators, and let us introduce the operator $A(\epsilon)=X+\epsilon Y$, for which the operator identity holds Snider (1964)

[TABLE]

Thus one has Feynman (1951); Snider (1964)

[TABLE]

By taking $A(\epsilon)=\pm\mathbbm{i}(q^{(0)}(t)+\epsilon q^{(1)}(t))$ and by applying the above result (88), one obtains the first order term of eq. (86), and thus eq. (74) in the main text.

Appendix D Proof of the equality $\langle{\ddot{x}^{(2)}}\rangle=0$

In order to prove that in the steady state the average center-of-mass velocity is constant we proceed as follows. By summing eqs. (46) and (47), we obtain

[TABLE]

whose solution reads

[TABLE]

and where the Fourier transform of $G_{x}(t)$ is given by eq. (25). The second derivative reads

[TABLE]

By averaging over the quantum noise, we thus obtain

[TABLE]

The average $-b^{2}\langle{h(t)}\rangle$ is given by the rhs of eq. (84), see also eq. (48). Therefore we have $-b^{2}\langle{h(t)}\rangle=b^{2}\mathcal{I}$, where $\mathcal{I}$ is given by eq. (85) and is time independent. Thus, as the systems is in the steady state $\langle{\partial_{t}h(t)}\rangle=\partial_{t}\langle{h(t)}\rangle=0$, and the rhs of eq. (92) vanishes, which is the desired result.