2023

[1,2]Luca Razzoli These authors contributed equally to this work. These
authors contributed equally to this work. These authors contributed equally to
this work. These authors contributed equally to this work. [1]Center for
Nonlinear and Complex Systems, Dipartimento di Scienza e Alta Tecnologia,
Università degli Studi dell’Insubria, Via Valleggio 11, Como, 22100, Italy
2]Istituto Nazionale di Fisica Nucleare, Sezione di Milano, Via Celoria 16,
Milano, 20133, Italy 3]Dipartimento di Fisica, Università di Genova, Via
Dodecaneso 33, Genova, 16146, Italy 4]CNR-SPIN, Via Dodecaneso 33, Genova,
16146, Italy 5]NEST, Istituto Nanoscienze-CNR, Piazza S. Silvestro 12, Pisa,
56127, Italy

# Efficiency and thermodynamic uncertainty relations of a dynamical quantum
heat engine

<EMAIL_ADDRESS>Fabio Cavaliere<EMAIL_ADDRESS>Matteo
Carrega<EMAIL_ADDRESS>Maura Sassetti<EMAIL_ADDRESS>Giuliano Benenti<EMAIL_ADDRESS>*  [  [  [  [

###### Abstract

In the quest for high-performance quantum thermal machines, looking for an
optimal thermodynamic efficiency is only part of the issue. Indeed, at the
level of quantum devices, fluctuations become extremely relevant and need to
be taken into account. In this paper we study the thermodynamic uncertainty
relations for a quantum thermal machine with a quantum harmonic oscillator as
a working medium, connected to two thermal baths, one of which is dynamically
coupled. We show that parameters can be found such that the machine operates
both as a quantum engine or refrigerator, with both sizeable efficiency and
small fluctuations.

###### keywords:

Quantum thermodynamics, Thermodynamic uncertainty relations, Out of
equilibrium quantum systems

##  1 Introduction

The rapid development of quantum technologies calls for a deeper understanding
of thermodynamics and energetics at a microscopic level, where unavoidable
quantum effects have to be taken into account. The extension of classical
concepts of thermodynamics to the quantum realm is not only at the forefront
of fundamental theoretical research campisi11 ; kosloff12 ; niedenzu ; anders
; talkner20 ; benenti17 , but also relevant for new nanodevice applications
bjorn ; giazotto ; martinez ; pekola ; arrachea . Mastering the thermodynamics
of quantum systems far from equilibrium, from energy storage to energy
transfer and transduction and heat–to–work conversion, is of great importance
for emerging technologies with applications in quantum computing arrachea ;
krantz19 ; calzona , quantum communication aufeve ; qcbook , and quantum
sensing golokolenov .

Spurred by the rapid emergence of new quantum technology platforms, in the
last few years the first prototype realizations of quantum thermal machines
pekola ; spin ; cangemireview ; sheng ; vischi  and quantum batteries quach ;
maillette  have been reported, calling for new theoretical and experimental
investigations. For instance, a great interest revolves around heat nano-
engines and refrigerators where the working medium (WM) is a quantum system
coupled to several thermal reservoirs. Indeed, it has been shown that a multi-
terminal configuration can improve both the output power and the efficiency at
maximum power levy12 ; huber19 ; imry12 ; mazza04 ; mazza15 ; erdman17 ;
sanchez2 ; cavaliere_iscience ; lopez23 ; lu23 . Also, the impact of quantum
coherence on the performance of thermal machines has been inspected karimi16 ;
pekola19 ; brandner17 ; kosloff02 ; rezek06 ; fried17 ; korze16 ; hammam22 .

Quantum devices are highly sensitive to fluctuations that may limit device
performances anders ; cangemireview ; paternostro ; gomez0 ; gomez , in
contrast to macroscopic devices, where usually fluctuations can be safely
neglected. To assess fluctuations, the thermodynamic uncertainty relations
(TURs) have been recently introduced barato ; shiraishi16 ; timpanaro ;
guarnieri ; asegava ; pal ; shiraishi18 ; cangemi ; menczel ; koyuk ; potanina
. These TURs combine steady-state currents $J$, their fluctuations $D_{J}$ and
dissipation (measured by the entropy production rate $\dot{S}$), giving limits
on the precision of currents for a given dissipation. Indeed, in Ref. barato
it was reported that for certain classical Markovian systems, the signal-to-
noise ratio satisfies a TUR bound, showing that relative current fluctuations
are lower bounded in terms of the inverse entropy production rate. By
introducing a dimensionless trade-off parameter $Q_{J}$, TURs can be expressed
as

$Q_{J}\equiv\dot{S}\frac{D_{J}}{J^{2}}\geq 2k_{B},$ (1)

where $k_{B}$ is the Boltzmann constant (from now on we will set $k_{B}=1$).
After that, several bounds for TURs horowitz17 ; pietzonka18 ; koyuk18 ;
brandner18 ; coherence18 ; cangemi20 ; jiao20 ; singh , and their possible
violations, have been put forward under different assumptions, starting from
classical Markovian dynamics to quantum systems, including time-dependent
forces and the role of time-reversal symmetry. For instance, when the model
preserves time–reversal–symmetry it has been shown that, at least in the
linear regime cangemi ; brandner18 , the trade-off parameter $Q$ never
achieves values lower than 2. Recently, experimental results on TURs have
appeared pal . They thus represent new tools to inspect the performances of
so-far unexplored resources based on non-thermal states, correlations and
quantum coherences in non-equilibrium quantum systems paternostro .

In this work, we consider a single quantum harmonic oscillator, a building
block of many quantum technologies krantz19 ; arrachea , as the WM coupled to
two thermal reservoirs. One of them is statically linked, while the other one
is modulated at a frequency $\Omega$. Driven system-bath couplings offer
enhanced design flexibility, such as the possibility to switch between regimes
of quantum heat engine or refrigerator cavaliere_prr ; cavaliere_iscience .
Here we are interested in assessing the precision of the currents in such
regimes, by studying the TURs related to the heat current and the work current
(i.e., the total exchanged power). Via a systematic perturbative approach, we
obtain expressions for the thermodynamic quantities and the TURs, that can be
numerically integrated with standard techniques. Even though we report no
violation of TURs, i.e. we do not observe values of $Q<2$, we show that $Q$
can attain minimal values (close to 2) both in the engine and refrigerator
regimes, accompanied by sizeable efficiency. In the limit of very small
damping, simple analytical expressions for the TURs are obtained, which
support our findings.

In Sec. 2 the model is illustrated, along with the definition of all
thermodynamical quantities and corresponding correlators describing their
fluctuations. We then proceed to Sec. 3 evaluating their quantum and time
averaged values in the dynamical steady state, to lowest order in the strength
of the dynamical coupling. Results are reported and discussed in Sec. 4, and
conclusions are drawn in Sec. 5.

##  2 Model and general setting

###  2.1 Model

Figure 1: Setup. Schematic depiction of the setup under investigation. The WM
– a quantum harmonic oscillator with characteristic frequency $\omega_{0}$ –
is in contact with two thermal reservoirs at temperatures $T_{\nu}$, with
$\nu=1,2$. The WM can exchange heat currents $J_{\nu}$ with the reservoirs and
total power $P$ with an external source. The purple dashed lines represent the
coupling between the WM and the thermal reservoirs. The coupling to the
$\nu=1$ reservoir is modulated in time with frequency $\Omega$, while that to
the $\nu=2$ reservoir is static.

We consider a minimal model for an operating heat engine in the quantum
regime. A single quantum harmonic oscillator acts as the WM and is coupled to
two thermal reservoirs, respectively kept at temperatures $T_{1}$ and $T_{2}$
–see Fig. 1. The total Hamiltonian is ($\hbar=k_{B}=1$)

$H^{(t)}=H_{{\rm WM}}+\sum_{\nu=1}^{2}\left[H_{\nu}+H_{{\rm
int},\nu}^{(t)}\right]\,,$ (2)

where the Hamiltonian of the WM is

$H_{{\rm WM}}=\frac{p^{2}}{2m}+\frac{1}{2}m\omega_{0}^{2}x^{2}\,,$ (3)

with $m$ and $\omega_{0}$ the mass and the characteristic frequency. The
Hamiltonian of the $\nu=1,2$ reservoir is described in the usual Caldeira-
Legget framework weiss ; CL83  in terms of an infinite set of independent
harmonic oscillators as

$H_{{\nu}}=\sum_{k=1}^{+\infty}\left[\frac{P^{2}_{k,\nu}}{2m_{k,\nu}}+\frac{1}{2}m_{k,\nu}\omega^{2}_{k,\nu}X^{2}_{k,\nu}\right]\,,$
(4)

and the interaction between the WM and the $\nu$–th reservoir is

$\\!\\!\\!\\!\\!\\!\\!\\!H_{{\rm
int},\nu}^{(t)}\\!\\!=\\!\\!\sum_{k=1}^{+\infty}\left[-g_{\nu}(t)c_{k,\nu}xX_{k,\nu}+\frac{g^{2}_{\nu}(t)c^{2}_{k,\nu}}{2m_{k,\nu}\omega_{k,\nu}^{2}}x^{2}\right]\\!,$
(5)

where $c_{k,\nu}$ describes the interaction strength and the counter–term
$\propto x^{2}$ prevents the renormalization of the WM potential. Following
Refs. carrega_prxquantum ; cavaliere_prr ; cavaliere_iscience , we focus on a
dynamical coupling situation to establish a working thermal machine in a
periodic steady state regime. In particular, we assume here that the
dimensionless coupling with the $\nu=1$ reservoir is modulated in time as
$g_{1}(t)=\cos(\Omega t)$, with $\Omega$ the driving frequency, while the
$\nu=2$ coupling is kept constant $g_{2}(t)=1$ – see the sketch in Fig. 1. At
the initial time $t_{0}{\to-\infty}$, the reservoirs are assumed in their
thermal equilibrium at temperatures $T_{\nu}$, with the total density matrix
written in a factorized form as

$\rho(t_{0})=\rho_{{\rm
WM}}(t_{0})\otimes\rho_{1}(t_{0})\otimes\rho_{2}(t_{0})\,,$ (6)

with

$\rho_{\nu}(t_{0})=\frac{e^{-H_{\nu}/T_{\nu}}}{{\rm
Tr}\left\\{e^{-H_{\nu}/T_{\nu}}\right\\}}$ (7)

and with $\rho_{{\rm WM}}(t_{0})$ the initial density matrix of the WM.

In the Heisenberg picture the equations of motion for the WM and reservoir
degrees of freedom carrega_prxquantum  are, respectively,

$\displaystyle\dot{x}(t)$ $\displaystyle=$ $\displaystyle\dfrac{p(t)}{m},$
$\displaystyle\dot{p}(t)$ $\displaystyle=$
$\displaystyle-m\omega_{0}^{2}x(t)+{\displaystyle\sum_{\nu=1}^{2}}{\displaystyle\sum_{k=1}^{+\infty}}\bigg{[}g_{\nu}(t)c_{k,\nu}X_{k,\nu}(t)$
(8) $\displaystyle-$
$\displaystyle\dfrac{g^{2}_{\nu}(t)c^{2}_{k,\nu}}{m_{k,\nu}\omega_{k,\nu}^{2}}x(t)\bigg{]},$

and

$\displaystyle\\!\\!\\!\\!\dot{X}_{k,\nu}(t)$ $\displaystyle=$
$\displaystyle\frac{P_{k,\nu}(t)}{m_{k,\nu}},$
$\displaystyle\\!\\!\\!\\!\dot{P}_{k,\nu}(t)$ $\displaystyle=$
$\displaystyle\\!\\!-m_{k,\nu}\omega_{k,\nu}^{2}X_{k,\nu}(t)+g_{\nu}(t)c_{k,\nu}x(t),$
(9)

where overdots denote time derivatives.

The solution for the position operator of the $k$–th oscillator of the
$\nu$–th reservoir are

$\displaystyle\\!\\!\\!\\!\\!\\!X_{k,\nu}(t)$ $\displaystyle=$
$\displaystyle\xi_{k,\nu}(t)+\frac{c_{k,\nu}}{m_{k,\nu}\omega_{k,\nu}}\cdot$
(10) $\displaystyle\cdot$ $\displaystyle\int_{t_{0}}^{t}\mathrm{d}s\
g_{\nu}(s)x(s)\sin[\omega_{k,\nu}(t-s)]\,,$
$\displaystyle\\!\\!\\!\\!\\!\\!\xi_{k,\nu}(t)$ $\displaystyle=$
$\displaystyle X_{k,\nu}(t_{0})\cos[\omega_{k,\nu}(t-t_{0})]$ (11)
$\displaystyle+$
$\displaystyle\frac{P_{k,\nu}(t_{0})}{m_{k,\nu}\omega_{k,\nu}}\sin[\omega_{k,\nu}(t-t_{0})]\,.$

The usual fluctuating force operator is defined as

$\xi_{\nu}(t)\equiv\sum_{k=1}^{+\infty}c_{k,\nu}\xi_{k,\nu}(t),$ (12)

with zero quantum average $\langle\xi_{\nu}(t)\rangle={\rm
Tr}\\{\xi_{\nu}(t)\rho(t_{0})\\}=0$ and correlation function
$\langle\xi_{\nu}(t)\xi_{\nu^{\prime}}(t^{\prime})\rangle\equiv{\cal
L}_{\nu}(t-t^{\prime})\delta_{\nu,\nu^{\prime}}$, with

${\cal L}_{\nu}(t-t^{\prime})={\cal L}_{\nu,{\rm s}}(t-t^{\prime})+{\cal
L}_{\nu,{\rm a}}(t-t^{\prime})\,.$ (13)

Here, the symmetric (s) and anti-symmetric (a) contributions with respect to
the time argument are

$\displaystyle\\!\\!\\!\\!\\!\\!{\cal L}_{\nu,{\rm s}}(t)$
$\displaystyle\\!\\!=$
$\displaystyle\\!\\!\\!\\!\int_{0}^{+\infty}\\!\\!\frac{\mathrm{d}\omega}{\pi}{\cal
J}_{\nu}(\omega)\coth\left(\frac{\omega}{2T_{\nu}}\right)\cos(\omega t),$ (14)
$\displaystyle\\!\\!\\!\\!\\!\\!{\cal L}_{\nu,{\rm a}}(t)$
$\displaystyle\\!\\!=$
$\displaystyle\\!\\!-i\int_{0}^{+\infty}\\!\\!\frac{\mathrm{d}\omega}{\pi}{\cal
J}_{\nu}(\omega)\sin(\omega t).$ (15)

In the above expressions, the so-called spectral density weiss

${\cal
J}_{\nu}(\omega)=\frac{\pi}{2}\sum_{k=1}^{+\infty}\frac{c^{2}_{k,\nu}}{m_{k,\nu}\omega_{k,\nu}}\delta(\omega-\omega_{k,\nu})$
(16)

has been introduced. It governs the properties of the $\nu$–th reservoir, e.g.
its possible non–Markovian behaviour nm1 ; nm2 ; nm3 ; nm4 ; nm5 . The precise
shape of ${\cal J}_{\nu}(\omega)$ will be specified later.

###  2.2 Thermodynamic quantities

In this work we focus on thermodynamic quantities in the long time limit, when
a periodic steady state has been reached. Therefore, all quantities of
interest can be averaged over one period of the drive ${\cal T}=2\pi/\Omega$,
and are therefore well-defined both for weak and strong coupling jurgen ;
esposito19 ; liu22 ; brandner166 . The average heat current $J_{\nu}$ of the
$\nu$–th reservoir and the total power $P$ can be written as cavaliere_prr

$\displaystyle P$ $\displaystyle=$ $\displaystyle\frac{1}{\cal
T}\int_{\bar{t}}^{\bar{t}+{\cal
T}}\\!\\!\\!\\!\mathrm{d}t^{\prime}\\!\\!\sum_{\nu=1}^{2}{\rm
Tr}\left[\frac{\partial H_{{\rm int},\nu}^{(t^{\prime})}}{\partial
t^{\prime}}\rho(t^{\prime})\right]\,,$ (17) $\displaystyle J_{\nu}$
$\displaystyle=$ $\displaystyle-\frac{1}{\cal T}\int_{\bar{t}}^{\bar{t}+{\cal
T}}\\!\\!\\!\\!\mathrm{d}t^{\prime}{\rm
Tr}\left[H_{\nu}\frac{\mathrm{d}}{\mathrm{d}t^{\prime}}\rho(t^{\prime})\right]\,,$
(18)

where here and in all the following $\bar{t}$ denotes a large positive time.
It is worth noting that a positive sign in the above quantities indicates an
energy flow towards the WM. Since the power contribution is associated to the
temporal variation of the interaction term, it is only due to the dynamical
coupling of the WM to the $\nu=1$ bath. An energy balance relation holds true,

$P+J_{1}+J_{2}=0,$ (19)

reflecting the first law of thermodynamics cavaliere_prr ; benenti17 ;
paternostro . An important quantity to assess thermodynamic performances is
the time–averaged entropy production rate, related to the currents by
benenti17 ; paternostro ; esposito19 )

$\dot{S}=-\sum_{\nu=1}^{2}\frac{J_{\nu}}{T_{\nu}}\,,$ (20)

with $\dot{S}\geq 0$ in accordance to the second law of thermodynamics
paternostro .

Due to dissipation, all quantities undergo fluctuations that affect the device
performances. During the whole time interval $t-t_{0}$, fluctuations can be
characterized by auto–correlation functions weiss ; cangemi : given an
observable $O(t)$ we define

$D_{O}(t)=\frac{1}{t-t_{0}}\int_{t_{0}}^{t}\mathrm{d}s\int_{t_{0}}^{t}\mathrm{d}s^{\prime}\langle{O(s)O(s^{\prime})}\rangle,$
(21)

where $t-t_{0}$ is a large time interval, $t-t_{0}\to+\infty$ and where
$\langle A(t)\rangle\equiv\mathrm{Tr}\\{A(t)\rho(t_{0})\\}$ denotes the
quantum average. In this limit, the above auto–correlation function reduces to
a single integral, that can be written as cangemi

$D_{O}(t)=\int_{0}^{+\infty}\mathrm{d}\tau\langle\\{O(t),O(t-\tau)\\}\rangle\,,$
(22)

with $\\{A,B\\}=AB+BA$ the anti–commutator, leading to the period–averaged
fluctuation

$D_{O}=\frac{1}{{\cal T}}\int_{\bar{t}}^{\bar{t}+{\cal T}}\ \mathrm{d}s\
D_{O}(s).$ (23)

##  3 Perturbative approach

The dynamics of the driven dissipative quantum system can be solved by
resorting to non-equilibrium Green function formalism paz1 ; paz2 ; arrachea12
, as described in Ref. carrega_prxquantum . In this work, however, we will
only focus on the case in which the modulated coupling ($\nu=1$) is much
weaker than the static one ($\nu=2$). Under these circumstances, a
perturbative expansion can be used to evaluate all the quantities introduced
above. To do so, it is convenient to rewrite the total Hamiltonian in Eq. (2)
as

$H^{(t)}=H^{(0)}+\Delta H^{(t)}\,,$ (24)

as the sum of an unperturbed Hamiltonian, $H^{(0)}$, and $\Delta H^{(t)}$, the
time–dependent interaction with the reservoir $\nu=1$:

$H^{(0)}=H_{{\rm WM}}+\sum_{\nu=1}^{2}H_{\nu}+H_{{\rm int},2}\,\mathchar
24635\relax\;\,\Delta H^{(t)}=H_{{\rm int},1}^{(t)}.$

In the Heisenberg picture, the position operator of the WM can be written as

$x(t)=x^{(0)}(t)+\Delta x(t),$ (25)

where $x^{(0)}(t)=e^{iH^{(0)}(t-t_{0})}x_{\mathrm{S}}e^{-iH^{(0)}(t-t_{0})}$,
with the label S denoting the Schrödinger picture, and where the perturbative
correction is

$\\!\\!\\!\\!\\!\\!\Delta x(t)=-i\int_{t_{0}}^{t}{\rm
d}s\\!\left[g_{1}(s)x^{(0)}(s)\xi_{1}(s),x^{(0)}(t)\right],$ (26)

neglecting terms $O\left({\Delta H^{(t)}}^{2}\right)$.

###  3.1 Heat current and total power

We start considering the average heat current $J_{1}$. From Eq. (18) we have
the general result carrega_prxquantum

$\displaystyle\\!\\!\\!\\!J_{1}$ $\displaystyle=$
$\displaystyle-\int_{\bar{t}}^{\bar{t}+{\cal T}}\frac{\mathrm{d}s}{{\cal
T}}g_{1}(s)\left\langle x(s)\left[\dot{\xi}_{1}(s)+\right.\right.$ (27)
$\displaystyle+$ $\displaystyle\left.\left.\int_{t_{0}}^{s}{\rm
d}s^{\prime}\,g_{1}(s^{\prime})x(s^{\prime})F_{1}(s-s^{\prime})\right]\right\rangle\,,$

where we have introduced

$F_{1}(t)=\sum_{k=1}^{+\infty}\frac{c_{k,1}^{2}}{m_{k,1}}\cos(\omega_{k,1}t)$
(28)

with Fourier transform $F_{1}(\omega)=2\omega{\cal J}_{1}(\omega)$. Expanding
$x(t)$ as in Eq. (25) we obtain

$\displaystyle J_{1}\\!\\!$ $\displaystyle=$
$\displaystyle\\!\\!\\!-\\!\\!\\!\int_{\bar{t}}^{\bar{t}+{\cal
T}}\\!\\!\frac{\mathrm{d}s}{{\cal T}}g_{1}(s)\bigg{\\{}\langle
x^{(0)}(s)\dot{\xi}_{1}(s)\rangle+\langle\Delta
x(s)\dot{\xi}_{1}(s)\rangle\\!\\!+$ (29) $\displaystyle+$
$\displaystyle\int_{t_{0}}^{s}{\rm d}s^{\prime}\,g_{1}(s^{\prime})\langle
x^{(0)}(s)x^{(0)}(s^{\prime})\rangle F_{1}(s-s^{\prime})\bigg{\\}}.$

We now perform the quantum average on the initial state of Eq. (6), exploiting
the fact that $\langle x^{(0)}(t)\dot{\xi}_{1}(t)\rangle=0$ and that

$\displaystyle\\!\\!\\!\\!\langle\Delta x(t)\dot{\xi}_{1}(t)\rangle$
$\displaystyle=$ $\displaystyle-i\int_{t_{0}}^{t}{\rm d}s\,g_{1}(s)\cdot$ (30)
$\displaystyle\cdot$ $\displaystyle\langle[x^{(0)}(s),x^{(0)}(t)]\rangle
Z_{1}(s-t),$

where we have defined

$Z_{1}(t-t^{\prime})\equiv\langle\xi_{1}(t)\dot{\xi}_{1}(t^{\prime})\rangle=\frac{\mathrm{d}}{\mathrm{d}t^{\prime}}\mathcal{L}(t-t^{\prime})\,,$
(31)

where the decomposition in symmetric and antisymmetric contributions

$Z_{1}(t-t^{\prime})=Z_{1,\mathrm{s}}(t-t^{\prime})+Z_{1,\mathrm{a}}(t-t^{\prime})\,,$
(32)

has been introduced, see Eq. (13). Finally, we arrive at the expression for
the heat current

$\displaystyle J_{1}$ $\displaystyle=$
$\displaystyle-\int_{\bar{t}}^{\bar{t}+{\cal T}}\frac{\mathrm{d}s}{{\cal
T}}\bigg{\\{}\int_{t_{0}}^{s}{\rm
d}s^{\prime}g_{1}(s)g_{1}(s^{\prime})F_{1}(s-s^{\prime})\cdot$ (33)
$\displaystyle\cdot$ $\displaystyle\langle
x^{(0)}(s)x^{(0)}(s^{\prime})\rangle-i\int_{t_{0}}^{s}{\rm
d}s^{\prime}\,g_{1}(s^{\prime})g_{1}(s)\cdot$ $\displaystyle\cdot$
$\displaystyle
Z_{1}(s^{\prime}-s)\langle[x^{(0)}(s^{\prime}),x^{(0)}(s)]\rangle\bigg{\\}}.$

To perform the average over ${\cal T}$ we observe that the correlation
function

$C(t,t^{\prime})=\langle x^{(0)}(t)x^{(0)}(t^{\prime})\rangle$ (34)

only depends, in the long time limit $t\gg t_{0}$, on the time difference
$t-t^{\prime}$, i.e. in the long time limit $C(t,t^{\prime})=C(t-t^{\prime})$.
This allows to change variable $\tau=s-s^{\prime}$, to perform the limit
$t_{0}\to-\infty$, and thus to obtain

$\displaystyle J_{1}$ $\displaystyle=$
$\displaystyle-\frac{1}{2}\int_{0}^{+\infty}\\!\\!\\!{\rm
d}\tau\left\\{{F_{1}(\tau)}\cos(\Omega\tau)C(\tau)\right.$ (35)
$\displaystyle\\!\\!\\!\\!-$
$\displaystyle\\!\\!\\!\\!\left.iZ_{1}(-\tau)\cos(\Omega\tau)\left[C(-\tau)-C(\tau)\right]\right\\}\,.~{}$

By direct inspection one finds that ${F_{1}(t)=-2iZ_{1,{\rm s}}(t)}$ and
introducing the response function $\chi_{0}(t)\equiv
im\theta(t)\langle[x^{(0)}(t),x^{(0)}(0)]\rangle$ (see vignale ) with
$\theta(t)$ the Heaviside step function, whose Fourier transform is

$\chi_{0}(\omega)=-\frac{1}{\omega^{2}-\omega_{0}^{2}+i\omega\gamma_{2}(\omega)}\,,$
(36)

where $\gamma_{2}(\omega)=\int_{-\infty}^{+\infty}\mathrm{d}te^{i\omega
t}\gamma_{2}(t)$ is the damping kernel of the $\nu=2$ bath in Fourier space
with

$\gamma_{\nu}(t)\equiv\sum_{k=1}^{+\infty}\frac{c_{k,{\nu}}^{2}}{m_{k,{\nu}}\omega_{k,{\nu}}^{2}}\cos(\omega_{k,{\nu}}t)\,,$
(37)

it is then possible to rewrite the average heat current as

$\\!\\!\\!\\!\\!\\!J_{1}=\int_{-\infty}^{+\infty}\frac{\mathrm{d}\omega}{4\pi
m}(\omega+\Omega){\cal
J}_{1}(\omega+\Omega)N(\omega)\chi^{\prime\prime}_{0}(\omega)\,,$ (38)

where here and in the following, the double prime denotes the imaginary part,
and

$N(\omega)=\coth\left(\frac{\omega+\Omega}{2T_{1}}\right)-\coth\left(\frac{\omega}{2T_{2}}\right)\,.$
(39)

Turning to the average total power of Eq. (17) we have the general expression
carrega_prxquantum

$\displaystyle P$ $\displaystyle=$ $\displaystyle\int_{\bar{t}}^{\bar{t}+{\cal
T}}\mathrm{d}s\bigg{\\{}-\dot{g}_{1}(s)\langle x(s)\xi_{1}(s)\rangle+$
$\displaystyle\\!\\!\\!\\!+$
$\displaystyle\\!\\!\\!\\!\dot{g}_{1}(s)\left\langle x(s)\int_{t_{0}}^{s}{\rm
d}s^{\prime}\,\gamma_{1}(s-s^{\prime})\frac{{\rm d}}{{\rm
d}s^{\prime}}[g_{1}(s^{\prime})x(s^{\prime})]\right\rangle\bigg{\\}}\,.$

Following the same steps of the perturbative expansion detailed above, after
the quantum average we get

$\displaystyle P$ $\displaystyle=$ $\displaystyle\int_{\bar{t}}^{\bar{t}+{\cal
T}}\frac{ds}{{\cal T}}\Big{\\{}2\int_{t_{0}}^{s}{\rm
d}s^{\prime}\,\Big{[}\mathcal{L}_{1}(s^{\prime}-s)\dot{g}_{1}(s)g_{1}(s^{\prime})\cdot$
(40) $\displaystyle\cdot$ $\displaystyle
C^{\prime\prime}(s-s^{\prime})\Big{]}-\int_{t_{0}}^{s}{\rm
d}s^{\prime}\,\Big{[}\dot{g}_{1}(s)g_{1}(s^{\prime})C(s-s^{\prime})\cdot$
$\displaystyle\cdot$ $\displaystyle\frac{{\rm d}}{{\rm
d}s^{\prime}}\gamma_{1}(s-s^{\prime})\Big{]}+\gamma_{1}(0)\dot{g}_{1}(s)g_{1}(s)C(0)\Big{\\}}\,.$

Taking the average over the period of the drive, recalling Eqs. (13), (14),
(15) and exploiting the identity $\dot{\gamma}_{1}(t)=-2i{\cal L}_{1,{\rm
a}}(t)$ one obtains

$P=-\Omega\int_{-\infty}^{+\infty}\frac{\mathrm{d}\omega}{4\pi m}{\cal
J}_{1}(\omega+\Omega)N(\omega)\chi^{\prime\prime}_{0}(\omega).$ (41)

Notice that the above results are fully consistent with those obtained using a
fully non-equilibrium Green function formalism in the perturbative regime (see
Appendix E of Ref. cavaliere_prr ). In the following we will focus only on
$J_{1}$ and $P$, but $J_{2}$ can be easily obtained from the energy balance
expressed by Eq. (19).

###  3.2 Correlation functions

The auto–correlation function given in Sec. 2.2 can be evaluated in the
perturbative regime following steps similar to those described above. Here we
briefly outline the procedure to obtain $D_{J_{1}}$: starting from the
definition of Eq. (22), to lowest order it is sufficient to consider only the
zero-th order term $x^{(0)}(t)$ of the position operator $x(t)$ in Eq. (25).
Then, the fluctuations of $J_{1}$ can be written as

$\displaystyle D_{J_{1}}(t)$ $\displaystyle=$
$\displaystyle\\!\\!\int_{0}^{+\infty}{\mathrm{d}}\tau
g_{1}(t)g_{1}(t-\tau)\cdot$ $\displaystyle\\!\\!\cdot$
$\displaystyle\\!\\!\left\langle\\{x^{(0)}(t)\dot{\xi}_{1}(t),x^{(0)}(t-\tau)\dot{\xi}_{1}(t-\tau)\\}\right\rangle\,.$

The quantum averages in the above expression decouple in terms of the form
$\langle x^{(0)}(t)x^{(0)}(t^{\prime})\rangle\equiv C(t,t^{\prime})$ – see Eq.
(34) – and in terms of the form

$\langle\dot{\xi}_{1}(t)\dot{\xi}_{1}(t^{\prime})\rangle=\frac{{\mathrm{d}}}{{\mathrm{d}}t}\frac{{\mathrm{d}}}{{\mathrm{d}}t^{\prime}}{\cal
L}_{1}(t-t^{\prime})\,.$ (42)

Exploiting again the time translational invariance
$C(t,t^{\prime})=C(t-t^{\prime})$ and following steps analogous to those
discussed above, the fluctuation of $J_{1}$ averaged over ${\cal T}$ – see Eq.
(23) – turn out to be

$\\!\\!\\!\\!\\!D_{J_{1}}\\!\\!=\\!\\!\int_{-\infty}^{+\infty}\\!\\!\frac{\mathrm{d}\omega}{4\pi
m}\\!(\omega+\Omega)^{2}{\cal
J}_{1}(\omega+\Omega)R(\omega)\chi^{\prime\prime}_{{0}}(\omega),$ (43)

with

$\\!\\!R(\omega)=\coth\left(\frac{\omega+\Omega}{2T_{1}}\right)\coth\left(\frac{\omega}{2T_{2}}\right)-1\,.$
(44)

In complete analogy, we derive the period–averaged fluctuations of the power
$P$ that read

$D_{P}=\Omega^{2}\int_{-\infty}^{+\infty}\frac{\mathrm{d}\omega}{4\pi
m}~{}{\cal J}_{1}(\omega+\Omega)R(\omega)\chi^{\prime\prime}_{{0}}(\omega)\,.$
(45)

##  4 Results

Hereafter, we focus on the specific case of a system with a structured,
Lorentzian spectral density for the dynamically coupled bath

${\mathcal{J}}_{1}(\omega)=\frac{d_{1}m\gamma_{1}\omega}{(\omega^{2}-\omega_{1}^{2})^{2}+\gamma_{1}^{2}\omega^{2}}\,,$
(46)

peaked at $\omega\approx\omega_{1}$, with typical broadening $\gamma_{1}$ and
coupling strength parameterized by $d_{1}$. Such a spectral density, which is
known to induce non–Markovian effects cavaliere_prr ; nm1 ; nm2 ; nm3 ; nm4 ;
nm5 , can be either obtained in cavity–optomechanical systems cavop1 ; cavop2
or in circuit-QED setups cqed1 ; cqed2 , with the latter representing a
particularly convenient solid–state platform which allows to tune dynamical
couplings up to a high degree of precision cqed3 . For the statically coupled
contact we choose a Ohmic spectral density
${\mathcal{J}}_{2}(\omega)=m\gamma_{2}\omega$, which implies
$\gamma_{2}(\omega)\equiv\gamma_{2}$ in Eq. (36).

###  4.1 Average thermodynamic quantities

The quantum-thermodynamical properties of this setup have been investigated
recently cavaliere_prr . It has been shown that it can either operate as a
quantum engine or a quantum refrigerator depending on the parameters. Here we
focus on the case $\gamma_{1,2}\ll\omega_{0}$, favourable to obtain better
performances and smaller fluctuations, and on temperatures in the quantum
regime $T_{\nu}<\omega_{0}$. In particular, we show here the case
$T_{2}=2T_{1}$, with $T_{1}=0.4\omega_{0}$.111To be consistent with the
perturbative approach, $d_{1}$ has been chosen small enough so that the peak
value of $\mathcal{J}_{1}(\omega)$ satisfies $\mathcal{J}_{1}(\omega_{1})\ll
m\omega_{0}^{2}$.

Figure 2: Energy currents and operating modes. Panels (a) and (b) show, as a
density plot, the average power $P$ and heat current $J_{1}$ as a function of
the driving frequency $\Omega$ and the resonance frequency $\omega_{1}$. Both
quantities are here normalized to $\gamma_{2}^{2}$. The dashed line marks the
resonance condition $\omega_{1}=\omega_{0}-\Omega$. In all panels
$T_{1}=0.4\omega_{0}$, $T_{2}=0.8\omega_{0}$, $\gamma_{1}=0.05\omega_{0}$,
$\gamma_{2}=0.01\omega_{0}$ and $d_{1}=10^{-3}\omega_{0}^{4}$.

Results for the total power $P$ and for the heat current $J_{1}$, obtained
integrating numerically Eqs. (41),(38) are shown as density plots in Figs.
2(a,b) as a function of the driving frequency $\Omega$ and $\omega_{1}$.
Depending on the sign of $J_{\nu}$ and $P$ different operating modes are
identified cavaliere_iscience . Here we focus on the most relevant ones: heat
engine with $P<0$ and $J_{1}<0$ and refrigerator with $P>0$ and $J_{1}>0$. The
best performance for these operating modes occurs, for $T_{1}<T_{2}$, when
$\omega_{1}=\omega_{0}-\Omega$ is satisfied cavaliere_prr , see the dashed
line.222The situation for $T_{1}>T_{2}$ is essentially mirrored along the axis
$\omega_{1}=\omega_{0}$. In this case the best engine and refrigerator
operating regimes occur along the resonance condition
$\omega_{1}=\omega_{0}+\Omega$ cavaliere_prr .

###  4.2 Thermodynamic uncertainty relations

It is now interesting to assess the impact of fluctuations on the performance
of this setup. To this end, we use TURs barato ; cangemi ; horowitz17 ;
brandner18 ; agarwalla . We are particularly interested in the engine and
refrigerator regimes and thus in the TURs for $P$ and $J_{1}$, expressed in
terms of the trade-off parameter $Q_{P}$ and $Q_{J_{1}}$ as cangemi ;
brandner18 ; agarwalla

$Q_{P}=\dot{S}\frac{D_{P}}{P^{2}}\geq 2\ \mathchar 24635\relax\;\
Q_{J_{1}}=\dot{S}\frac{D_{J_{1}}}{J_{1}^{2}}\geq 2\,.$ (47)

These TURs combine energy flows, their fluctuations, and the entropy
production rate in a dimensionless quantity expressing the trade-off between
how the system fluctuates with respect to the quality (in terms of magnitude
and degree of dissipation) of the energy flow. The lower $Q_{\mu}$ (with
$\mu\in\\{J_{1},P\\}$), the better the operation of the thermal machine
pietzonka18 .

Figure 3: Quantifying fluctuations via TURs. Density plot of $Q_{P}$ (panel a)
and $Q_{1}$ (panel b) as a function of $\Omega$ and $\omega_{1}$. Parameters
as in Fig. 2.

Figure 3 shows density plots of the numerically evaluated $Q_{P}$ (a) and
$Q_{J_{1}}$ (b) as a function of $\Omega$ and $\omega_{1}$ and it is clear
that $Q_{\mu}>2$ everywhere. Via extensive numerical investigations we can
report that regardless of all parameter configurations, $Q_{\mu}<2$ is never
found in this model, and thus no violation of TUR bounds barato ; cangemi ;
brandner18  occurs. The trade-off parameters can attain very large values, see
the white regions in the plots, and even diverge – e.g. at the crossover
between the engine and the refrigerator regimes along
$\omega_{1}=\omega_{0}-\Omega$ (due to $J_{1}$, $P$ crossing zero).
Interestingly, however, they are instead particularly low (even approaching
$Q_{\mu}\approx 2$) around the latter resonance line, where the best
performances occur. Let us then relate $Q_{\mu}$ to the efficiency
$\eta=-P/J_{2}$ (engine mode) or to the coefficient of performance
$\mathrm{COP}=J_{1}/P$ (refrigerator mode), normalized to the respective
Carnot limits $\eta_{C}=1-\frac{T_{1}}{T_{2}}$ and
$\mathrm{COP}_{C}=\frac{T_{1}}{T_{2}-T_{1}}$. 333For the parameters chosen
here, $\eta_{C}=0.5$ and COPC=1.

Figure 4: Efficiency vs. fluctuations. Panel (a) shows as a red curve the
efficiency, normalized to the Carnot limit, of the thermal machine operating
as an engine, and $Q_{P}$ as a green line, as a function of $\Omega$ along the
resonance line $\omega_{1}=\omega_{0}-\Omega$ – see Fig. 2. Panel (b) shows
the same quantities along the same resonance, but in the regime where the
machine operates as a refrigerator. Parameters as in Fig. 2.

Figure 4(a) shows $\eta/\eta_{C}$ (red line) and $Q_{P}$ (green line) as a
function of $\Omega$ within the engine region along the dashed line of Fig. 2.
As is clear, there is a large range of driving frequencies where $Q_{P}$ is
very close to 2. Furthermore, when it attains its minimum value, for
$\Omega\approx 0.35\omega_{0}$, the efficiency of the engine is $\approx
0.7\eta_{C}$. On the other hand, when the engine operates at its maximum
efficiency of $\eta\approx 0.95\eta_{C}$, $Q_{P}$ almost doubles. Notice also
that $Q_{P}$ diverges at the boundaries of the engine regime and that the
maximum of $\eta$ is located in close proximity of one of these edges. The
situation in the refrigerator regime is qualitatively similar, and is reported
in Fig. 4(b).

To gain further insight, we can exploit the fact that in the limit
$\gamma_{2}\to 0$ one can approximate cavaliere_prr
$\chi^{\prime\prime}_{0}(\omega)\approx\frac{\pi}{2\omega_{0}}\sum_{p=\pm
1}p\delta(\omega-p\omega_{0})$, with $\delta(\omega)$ the Dirac delta. The
integrals in Eqs. (38),(41),(43),(45) are then easily solved, leading to

$\displaystyle J_{1}$ $\displaystyle=$
$\displaystyle\frac{1}{8m\omega_{0}}\sum_{p=\pm
1}(\omega_{0}+p\Omega)\mathcal{J}_{1}(\omega_{0}+p\Omega)N_{p},$
$\displaystyle P$ $\displaystyle=$
$\displaystyle-\frac{\Omega}{8m\omega_{0}}\sum_{p=\pm
1}p\mathcal{J}_{1}(\omega_{0}+p\Omega)N_{p},$ $\displaystyle D_{J_{1}}$
$\displaystyle=$ $\displaystyle\frac{1}{8m\omega_{0}}\sum_{p=\pm
1}(\omega_{0}+p\Omega)^{2}\mathcal{J}_{1}(\omega_{0}+p\Omega)R_{p},$
$\displaystyle D_{P}$ $\displaystyle=$
$\displaystyle\frac{\Omega^{2}}{8m\omega_{0}}\sum_{p=\pm
1}\mathcal{J}_{1}(\omega_{0}+p\Omega)R_{p}\,,$

where $N_{p}=N(p\omega_{0})$ and $R_{p}=R(p\omega_{0})$ – see Eqs. (39),(44).
In the regime $\gamma_{1}\ll\omega_{0}$ considered in this paper, restricting
along the resonance line $\omega_{1}=\omega_{0}-\Omega$ one can safely drop
all terms with $p=+1$ in the above expressions cavaliere_prr . Then, after
simple rearrangements one can see that in this limit $Q_{P}\equiv Q_{J_{1}}$,
with

$Q_{\mu}\\!=\\!2\left(\frac{\omega_{0}}{2T_{2}}-\frac{\omega_{0}-\Omega}{2T_{1}}\right)\coth\left(\frac{\omega_{0}}{2T_{2}}-\frac{\omega_{0}-\Omega}{2T_{1}}\right).$

This proves that, at least in the above limits, the TURs are never violated
and $Q_{\mu}\geq 2$ in agreement with our numerical findings. Indeed, in this
limit $Q_{\mu}=2$ precisely when $\Omega=\Omega^{*}=\omega_{0}\eta_{C}$, which
is the “turning point” between the operation as an engine and that as a
refrigerator along the resonance cavaliere_prr . There, at
$\Omega=\Omega^{*}$, we have $\eta\to\eta_{C}$ (with $P\to 0$) and
COP$\to$COPC (with $J_{1}\to 0$). The net coalescence of engine and
refrigerator modes is lifted at finite values of $\gamma_{\nu}$, which results
in $\eta<\eta_{C}$ (COP$<$COPC), with the maximum for $\eta$ (COP) occurring
at $\Omega\lesssim\Omega^{*}$ ($\Omega\gtrsim\Omega^{*}$). Also, when
$\gamma_{\nu}$ is small but non zero the trade-off parameters are no longer
identical, and develop the divergence near $\Omega^{*}$ seen in Fig. 4.

##  5 Conclusions

In this paper we have analyzed the fluctuations of the heat current and power
exchanged between a quantum thermal machine with a single quantum harmonic
oscillator as the WM and two thermal baths at uneven temperature, in the
presence of a dynamical coupling between the WM and one of the baths. Such a
machine can operate either as a quantum heat engine or as refrigerator
depending on the frequency of the drive and the other parameters. To
understand the impact of fluctuations on such operating modes we have
developed expressions for the fluctuations to leading order in the amplitude
of the dynamical coupling. Such expressions allow us to numerically evaluate
the trade-off parameters $Q_{P}$ and $Q_{J_{1}}$ associated to total power (in
the engine case) or heat current from the cold bath (in the refrigerator
case), whose magnitude is associated to the strength and impact of
fluctuations. We report no violation of the lower bound proposed for these
quantities in the literature. However, we have found that in typical operating
regimes the trade-off parameters can reach values very close to the lower
bound, with an efficiency of about 70% of the Carnot limit in the case of a
quantum engine. In the limit of weak damping, appropriate for the regime of
parameters considered in this paper, we have also developed an analytical
approximation for the trade-off parameters, which supports our findings and
shows that TUR bounds cannot be violated in our model, at least for weak
damping.

Despite no violation of TURs has been found, showing that the trade-off
parameters can attain small values in sensible operating regimes is a first
step towards the development of optimal protocols for the operation of quantum
thermodynamic machines. Future developments of this work may lead to consider
more complex driving protocols for the thermal machine, which could help to
even increase the efficiency at the lowest values of the trade-off parameters.
Also, it will be interesting to look for more complex quantum thermal
machines, either in multi–terminal or multi–WM configurations. Due to their
ability to perform multiple thermodynamical tasks at once cavaliere_iscience ,
assessing the impact of fluctuations can be a key issue towards their
optimization.

Acknowledgments L.R. and G.B. acknowledge financial support by the Julian
Schwinger Foundation (Grant JSF-21-04-0001) and by INFN through the project
‘QUANTUM’. F.C. and M.S. acknowledge support by the “Dipartimento di
Eccellenza MIUR 2018-2022.”

Data availability All data generated during this study have been obtained by
numerically implementing the equations presented in the text using standard
techniques.

## References

  * (1) M. Campisi, P. Hänggi, and P. Talkner, Colloquium: Quantum fluctuation relations: Foundations and applications, Rev. Mod. Phys. 83, 771 (2011).
  * (2) R. Kosloff, Quantum thermodynamics: A dynamical Viewpoint, Entropy 15, 2100 (2013).
  * (3) D. Gelbwaser-Klimovsky, W. Niedenzu, and G. Kurizki, Thermodynamics of quantum systems under dynamical control, Adv. At. Mol. Opt. Phys. 64, 329 (2015).
  * (4) S. Vinjanampathy and J. Anders, Quantum thermodynamics, Contemp. Phys. 57, 545 (2016).
  * (5) P. Talkner and P. Hänggi, Colloquium: Statistical mechanics and thermodynamics at strong coupling: Quantum and classical, Rev. Mod. Phys. 92, 041002 (2020).
  * (6) G. Benenti, G. Casati, K. Saito, and R. S. Whitney, Fundamental aspects of steady-state conversion of heat to work at the nanoscale, Phys. Rep. 694, 1 (2017).
  * (7) B. Sothmann, R. Sánchez, and A. N. Jordan, Thermoelectric energy harvesting with quantum dots, Nanotechnology 26, 032001 (2015).
  * (8) F. Giazotto, T. T. Heikkilä, A. Luukanen, A. M. Savin, and J. P. Pekola, Opportunities for mesoscopics in thermometry and refrigeration: Physics and applications, Rev. Mod. Phys. 78, 217 (2006).
  * (9) M. J. Martínez-Pérez and F. Giazotto, A quantum diffractor for thermal flux, Nat. Commun. 5, 3579 (2014).
  * (10) J. P. Pekola and B. Karimi, Colloquium: Quantum heat transport in condensed matter systems, Rev. Mod. Phys. 93, 041001 (2021).
  * (11) L. Arrachea, Energy dynamics, heat production and heat–work conversion with qubits: toward the development of quantum machines, Rep. Prog. Phys. 86, 036501 (2023).
  * (12) P. Krantz, M. Kjaergaard, F. Yan, T. P. Orlando, S. Gustavsson, and W. D. Oliver, A quantum engineer’s guide to superconducting qubits, Appl. Phys. Rev. 6, 021318 (2019).
  * (13) A. Calzona and M. Carrega, Multi-mode architectures for noise-resilient superconducting qubits, Supercond. Sci. Technol. 36 023001 (2023).
  * (14) A. Auffèves, Quantum Technologies Need a Quantum Energy Initiative, PRX Quantum 3, 020101 (2022).
  * (15) G. Benenti, G. Casati, D. Rossini, and G. Strini, Principles of quantum computation and information (A comprehensive textbook) (World Scientific, Singapore, 2019).
  * (16) I. Golokolenov et al., Thermodynamics of a single mesoscopic phononic mode, Phys. Rev. Research 5, 013046 (2023).
  * (17) D. von Lindenfels et al., Spin Heat Engine Coupled to a Harmonic-Oscillator Flywheel, Phys. Rev. Lett. 123, 080602 (2019).
  * (18) L. M. Cangemi, C. Bhadra, and A. Levy, Quantum engines and refrigerators, Arxiv:2302.00726.
  * (19) J. Sheng, C. Yang, and H. Wu, Realization of a coupled-mode heat engine with cavity-mediated nanoresonators, Sci. Adv. 7, eabl7740 (2021).
  * (20) F. Vischi, M. Carrega, P. Virtanen, E. Strambini, A. Braggio, and F. Giazotto, Thermodynamic cycles in Josephson junctions, Sci. Rep. 9, 3238 (2019).
  * (21) J. Q. Quach et al., Superabsorption in an organic microcavity: toward a quantum battery, Sci. Adv. 8, eabk3160 (2022).
  * (22) I. Maillette de Buy-Wenniger et al., Coherence-powered work exchanges between a solid-state qubit and light fields, Arxiv:2202.01109
  * (23) A. Levy, R. Alicki, and R. Kosloff, Quantum refrigerators and the third law of thermodynamics, Phys. Rev. E 85, 061126 (2012).
  * (24) F. Clivaz, R. Silva, G. Haack, J.B. Brask, N. Brunner, and M. Huber, Unifying Paradigms of Quantum Refrigeration: A Universal and Attainable Bound on Cooling, Phys. Rev. Lett. 123, 170605 (2019).
  * (25) J.-H. Jiang, O. Entin-Wohlman, and Y. Imry, Thermoelectric three-terminal hopping transport through one-dimensional nanosystems, Phys. Rev. B 85, 075412 (2012).
  * (26) F. Mazza, R. Bosisio, G. Benenti, V. Giovannetti, R. Fazio, and F. Taddei, Thermoelectric efficiency of three-terminal quantum thermal machines, New J. Phys. 16, 085001 (2014).
  * (27) F. Mazza, S. Valentini, R. Bosisio, G. Benenti, V. Giovannetti, R. Fazio, and F. Taddei, Separation of heat and charge currents for boosted thermoelectric conversion, Phys. Rev. B. 91, 245435 (2015).
  * (28) P. A. Erdman, F. Mazza, R. Bosisio, G. Benenti, R. Fazio, and F. Taddei, Thermoelectric properties of an interacting quantum dot based heat engine, Phys. Rev. B 95, 245432 (2017).
  * (29) J. Yang, C. Elouard, J. Splettstoesser, B. Sothmann, R. Sánchez, and A. N. Jordan, Thermal transistor and thermometer based on Coulomb-coupled conductors, Phys. Rev. B 100, 045418 (2019).
  * (30) F. Cavaliere, L. Razzoli, M. Carrega, G. Benenti, and M. Sassetti, Hybrid quantum thermal machines with dynamical couplings, iScience 26, 106235 (2023).
  * (31) R. López, J.S. Lim, and K.W. Kim, Optimal superconducting hybrid machine, Phys. Rev. Research 5, 013038 (2023).
  * (32) J. Lu, Z. Wang, R. Wang, J. Peng, C. Wang, and J.-H. Jiang, Multitask quantum thermal machines and cooperative effects, Phys. Rev. B 107, 075428 (2023).
  * (33) B. Karimi and J. P. Pekola, Otto refrigerator based on a superconducting qubit: Classical and quantum performance, Phys. Rev. B 94, 184503 (2016).
  * (34) J. P. Pekola, B. Karimi, G. Thomas, and D. V. Averin, Supremacy of incoherent sudden cycles, Phys. Rev. B 100, 085405 (2019).
  * (35) K. Brandner, M. Bauer, and U. Seifert, Universal Coherence-Induced Power Losses of Quantum Heat Engines in Linear Response, Phys. Rev. Lett. 119, 170602 (2017).
  * (36) R. Kosloff and T. Feldmann, Discrete four-stroke quantum heat engine exploring the origin of friction, Phys. Rev. E 65, 055102(R) (2002).
  * (37) Y. Rezek and R. Kosloff, Irreversible performance of a quantum harmonic heat engine, New J. Phys. 8, 83 (2006).
  * (38) A. Friedenberger and E. Lutz, When is a quantum heat engine quantum?, Europhys. Lett. 120, 10002 (2017).
  * (39) K. Korzekwa, M. Lostaglio, J. Oppenheim, and D. Jennings, The extraction of work from quantum coherence, New J. Phys. 18, 023045 (2016).
  * (40) K. Hammam, H. Leitch, Y. Hassouni, and G. De Chiara, Exploiting coherence for quantum thermodynamic advantage, New J. Phys. 24, 113053 (2022).
  * (41) G.T. Landi, and M. Paternostro, Irreversible entropy production: From classical to quantum, Rev. Mod. Phys. 93, 035008 (2021).
  * (42) S. Hernández-Gómez, et al., Experimental test of exchange fluctuation relations in an open quantum system, Phys. Rev. Research 2, 023327 (2020).
  * (43) S. Hernández-Gómez, N. Staudenmaier, M. Campisi, and N. Fabbri, Experimental test of fluctuation relations for driven open quantum systems with an NV center, New J. Phys. 23, 065004 (2021).
  * (44) A. C. Barato and U. Seifert, Thermodynamic Uncertainty Relation for Biomolecular Processes, Phys. Rev. Lett. 114, 158101 (2015).
  * (45) N. Shiraishi, K. Saito, and H. Tasaki, Universal Trade-Off Relation between Power and Efficiency for Heat Engines, Phys. Rev. Lett. 117, 190601 (2016).
  * (46) A. M. Timpanaro, G. Guarnieri, J. Goold, and G. T. Landi, Thermodynamic Uncertainty Relations from Exchange Fluctuation Theorems, Phys. Rev. Lett. 123, 090604 (2019).
  * (47) G. Guarnieri, G. T. Landi, S. R. Clark, and J. Goold, Thermodynamics of precision in quantum nonequilibrium steady states, Phys. Rev. Research 1, 033021 (2019).
  * (48) Y. Hasegawa, Thermodynamic uncertainty relation for general open quantum systems, Phys. Rev. Lett. 126, 010602 (2021).
  * (49) S. Pal et al., Experimental study of the thermodynamic uncertainty relation, Phys. Rev. Research 2, 022044 (2020).
  * (50) N. Shiraishi and K. Saito, Fundamental Relation Between Entropy Production and Heat Current, J. Stat. Phys. 174, 433 (2019).
  * (51) L. M. Cangemi, M. Carrega, A. De Candia, V. Cataudella, G. De Filippis, M. Sassetti and G. Benenti, Optimal energy conversion through antiadiabatic driving breaking time-reversal symmetry, Phys. Rev. Research 3, 013237 (2021).
  * (52) P. Menczel, E. Loisa, K. Brandner, and C. Flindt, Thermodynamic uncertainty relations for coherently driven open quantum systems, J. Phys. A: Math. Theor. 54, 314002 (2021).
  * (53) T. Koyuk and U. Seifert, Thermodynamic Uncertainty Relation for Time-Dependent Driving, Phys. Rev. Lett. 125, 260604 (2020).
  * (54) E. Potanina, C. Flindt, M. Moskalets, and K. Brandner, Thermodynamic bounds on coherent transport in periodically driven conductors, Phys. Rev. X 11, 021013 (2021).
  * (55) J. M. Horowitz and T. R. Gingrich, Proof of the finite-time thermodynamic uncertainty relation for steady-state currents, Phys. Rev. E 96, 020103 (2017).
  * (56) P. Pietzonka and U. Seifert, Universal Trade-Off between Power, Efficiency, and Constancy in Steady-State Heat Engines, Phys. Rev. Lett. 120, 190602 (2018).
  * (57) T. Koyuk, U. Seifert, and P. Pietzonka, A generalization of the thermodynamic uncertainty relation to periodically driven systems, J. Phys. A: Math. Theor. 52, 02LT02 (2019).
  * (58) K. Brandner, T. Hanazato, and K. Saito, Thermodynamic Bounds on Precision in Ballistic Multiterminal Transport, Phys. Rev. Lett. 120, 090601 (2018).
  * (59) K. Ptaszyński, Coherence-enhanced constancy of a quantum thermoelectric generator, Phys. Rev. B 98, 085425 (2018).
  * (60) L. M. Cangemi, V. Cataudella, G. Benenti, M. Sassetti, and G. De Filippis, Violation of thermodynamics uncertainty relations in a periodically driven work-to-work converter from weak to strong dissipation, Phys. Rev. B 102, 165418 (2020).
  * (61) G. Benenti, G. Casati, and J. Wang, Power, efficiency, and fluctuations in steady-state heat engines, Phys. Rev. E 102, 040103(R) (2020).
  * (62) V. Singh, V. Shaghaghi, Ö. E. Müstecaplıoǧlu, and D. Rosa, Thermodynamic uncertainty relation in nondegenerate and degenerate maser heat engines, ArXiv:2211.08377.
  * (63) F. Cavaliere, M. Carrega, G. De Filippis, V. Cataudella, G. Benenti, and M. Sassetti, Dynamical heat engines with non-Markovian reservoirs, Phys. Rev. Research 4, 033233 (2022).
  * (64) U. Weiss, Quantum Dissipative Systems - 5th Edition (World Scientific, Singapore, 2021).
  * (65) A. O. Caldeira and A. J. Leggett, Quantum tunnelling in a dissipative system, Ann. Phys. 149, 374 (1983).
  * (66) M. Carrega, L. M. Cangemi, G. De Filippis, V. Cataudella, G. Benenti, and M. Sassetti, Engineering Dynamical Couplings for Quantum Thermodynamic Tasks, PRX Quantum 3, 010323 (2022).
  * (67) M. Thorwart, E. Paladino, and M. Grifoni, Dynamics of the spin-boson model with a structured environment, Chem. Phys. 296, 333 (2004).
  * (68) E. Paladino, A. G. Maugeri, M. Sassetti, G. Falci, and U. Weiss, Structured environments in solid state systems: Crossover from Gaussian to non-Gaussian behavior, Phys. E 40, 198 (2007).
  * (69) J. Iles-Smith, N. Lambert, and A. Nazir, Environmental dynamics, correlations, and the emergence of noncanonical equilibrium states in open quantum systems, Phys. Rev. A 90, 032114 (2014).
  * (70) P. Strasberg, G. Schaller, N. Lambert, and T. Brandes, Nonequilibrium thermodynamics in the strong coupling and non-Markovian regime based on a reaction coordinate mapping, New J. Phys. 18, 073007 (2016).
  * (71) S. Restrepo, J. Cerrillo, P. Strasberg, and G. Schaller, From quantum heat engines to laser cooling: Floquet theory beyond the Born-Markov approximation, New J. Phys. 20, 053063 (2018).
  * (72) M. Wiedmann, J. T. Stockburger, and J. Ankerhold, Non-Markovian dynamics of a quantum heat engine: out-of-equilibrium operation and thermal coupling control, New J. Phys. 22, 033007 (2020).
  * (73) K. Ptaszyński and M. Esposito, Entropy Production in Open Systems: The Predominant Role of Intraenvironment Correlations, Phys. Rev. Lett. 123, 200603 (2019).
  * (74) J. Liu, K. A. Jung, and D. Segal, Periodically Driven Quantum Thermal Machines from Warming up to Limit Cycle, Phys. Rev. Lett. 127, 200602 (2021).
  * (75) K. Brandner and U. Seifert, Periodic thermodynamics of open quantum systems, Phys. Rev. E 93, 062134 (2016).
  * (76) N. Freitas and J. P. Paz, Fundamental limits for cooling of linear quantum refrigerators, Phys. Rev. E 95, 012146 (2017).
  * (77) N. Freitas and J. P. Paz, Cooling a quantum oscillator: A useful analogy to understand laser cooling as a thermodynamical process, Phys. Rev. A 97, 032104 (2018).
  * (78) L. Arrachea, E. R. Mucciolo, C. Chamon, and R. B. Capaz, Microscopic model of a phononic refrigerator, Phys. Rev. B 86, 125424 (2012).
  * (79) G.F. Giuliani and G. Vignale, Quantum Theory of the Electron Liquid (Cambridge University Press, Cambridge, 2005).
  * (80) M. Aspelmeyer, T. J. Kippenberg, and F. Marquardt, Cavity optomechanics, Rev. Mod. Phys. 86, 1391 (2014).
  * (81) S. Barzanjeh, A. Xuereb, S. Gröblacher, M. Paternostro, C. A. Regal, and E. M. Weig, Optomechanics for quantum technologies, Nat. Phys. 18, 15 (2022).
  * (82) A. Cottet et al., Cavity QED with hybrid nanocircuits: from atomic-like physics to condensed matter phenomena, J. Phys.: Condens. Matter 29, 433002 (2017).
  * (83) B. Peropadre, D. Zueco, F. Wulschner, F. Deppe, A. Marx, R. Gross, and J.J. García-Ripoll, Tunable coupling engineering between superconducting resonators: from sidebands to effective gauge fields, Phys. Rev. B 87, 134504 (2013).
  * (84) F. Wulschner et al., Tunable coupling of transmission-line microwave resonators mediated by an rf SQUID, EPJ Quant. Technol. 3, 10 (2016).
  * (85) B. K. Agarwalla and D. Segal, Assessing the validity of the thermodynamic uncertainty relation in quantum systems, Phys. Rev. B 98, 155438 (2018).


# Quantum-enhanced bosonic learning machine

Chi-Huan Nguyen  These authors contributed equally to this work  Ko-Wei Tseng
These authors contributed equally to this work  Gleb Maslennikov222Present
address: NKT Photonics, Bregnerødvej 144, 3460 Birkerød, Denmark H. C. J. Gan
Centre for Quantum Technologies, National University of Singapore, 3 Science
Dr 2, 117543, Singapore  Dzmitry Matsukevich  Centre for Quantum Technologies,
National University of Singapore, 3 Science Dr 2, 117543, Singapore
Department of Physics, National University of Singapore, 2 Science Dr 3,
117551, Singapore

(August 27, 2024)

###### Abstract

Quantum processors enable computational speedups for machine learning through
parallel manipulation of high-dimensional vectors Biamonte  _et al._ (2017).
Early demonstrations of quantum machine learning have focused on processing
information with qubits Cai  _et al._ (2015); Trávníček  _et al._ (2019);
Havlíček  _et al._ (2019); Johri  _et al._ (2020); Bartkiewicz  _et al._
(2020); Xin  _et al._ (2021); Peters  _et al._ (2021). In such systems, a
larger computational space is provided by the collective space of multiple
physical qubits. Alternatively, we can encode and process information in the
infinite dimensional Hilbert space of bosonic systems such as quantum harmonic
oscillators Braunstein and van Loock (2005); Ortiz-Gutiérrez  _et al._ (2017);
Flühmann  _et al._ (2019). This approach offers a hardware-efficient solution
with potential quantum speedups to practical machine learning problems. Here
we demonstrate a quantum-enhanced bosonic learning machine operating on
quantum data with a system of trapped ions. Core elements of the learning
processor are the universal feature-embedding circuit that encodes data into
the motional states of ions, and the constant-depth circuit that estimates
overlap between two quantum states. We implement the unsupervised $K$-means
algorithm to recognize a pattern in a set of high-dimensional quantum states
and use the discovered knowledge to classify unknown quantum states with the
supervised $k$-NN algorithm. These results provide building blocks for
exploring machine learning with bosonic processors.

Machine learning techniques excel at finding patterns in data that are usually
encoded as high dimensional feature vectors. As the size of these vectors
grows, processing them classically requires enormous computational resources
Lloyd  _et al._ (2013). Recently proposed quantum algorithms show evidences of
computational speedups, that comes from the exploitation of large Hilbert
space of quantum systems Lloyd  _et al._ (2013, 2014, 2016); Rebentrost  _et
al._ (2014). Demonstrations of quantum machine learning algorithms implemented
on discrete variables platforms have shown tremendous potential Cai  _et al._
(2015); Trávníček  _et al._ (2019); Havlíček  _et al._ (2019); Johri  _et al._
(2020); Bartkiewicz  _et al._ (2020); Xin  _et al._ (2021). However, the
number of gates and qubits must be increased for solving practical problems on
near-term quantum processors.

An alternative approach is to utilize the infinite-dimensional Hilbert space
of bosonic systems as the feature space for quantum machine learning Lau  _et
al._ (2017); Schuld and Killoran (2019). One advantage of the bosonic
architecture is that a single operation, such as squeezing or displacement, is
formally equivalent to processing infinitely-dimensional feature vectors
Schuld and Killoran (2019). Furthermore, it allows for efficient use of
physical resources. For discrete-variable schemes, approximately $\log(n)$
physical qubits are required to encode an $n$-bit data string, whereas in the
alternate approach, one bosonic mode is sufficient. As an example, in the
ideal case the motion of a single trapped ion suffices to encode 3 infinitely
dimensional feature vectors, whereas the internal state of an ion typically
only encodes a single qubit Ortiz-Gutiérrez  _et al._ (2017).

Figure 1:  Experimental implementation. a. Schematic of the linear rf trap
with two trapped 171Yb+ ions. One ion (black) is prepared in the
${}^{2}F_{7/2}$ metastable state such that it does not interact with optical
fields applied during the experiment. Two Raman beams (R1,R2) couple the
motion in the radial directions with the internal states
$|g\rangle,|e\rangle$. The $x$ and $y$ axes point along the diagonally
opposite trap electrodes (Methods). b. Energy level diagram of 171Yb+. The
internal states used in the experiment are
$|g\rangle=|^{2}S_{1/2},F=0,m_{F}=0\rangle$ and
$|e\rangle=|^{2}S_{1/2},F=1,m_{F}=0\rangle$. The arrows represent the Raman
beams. c. Schematic of radial motional modes. $\omega_{A}$, $\omega_{B}$, and
$\omega_{C}$ are frequencies of the motional modes $A$, $B$, and $C$ which
denotes the out-of-phase motion in the $y$ direction, out-of-phase and in-
phase motion in the $x$ direction respectively.
$\omega_{AB}=|\omega_{A}-\omega_{B}|,\,\omega_{AC}=|\omega_{A}-\omega_{C}|$ d.
The experimental sequence to implement the bosonic learning machine. The
motional modes are initialized in the ground state. The embedding circuit
prepares arbitrary $n$-dimensional states in mode $B$ and $C$ by applying
$n-1$ pairs of carrier ($U_{\textrm{c}}$) and red sideband pulses
($U_{\textrm{rsb}}$). The gate sequence for the SWAP test consists of two
controlled beam splitters $U_{\textrm{CBS}}(\theta,\psi)$ followed by a
projective measurement of the internal state.

In this letter, we report a bosonic learning machine that consists of two
building blocks for encoding and processing either quantum or classical data
(see Fig. 1d). First, a universal embedding circuit implements a feature map
$\phi$ of classical data $\boldsymbol{x}=(x_{1},x_{2},...,x_{n})$ to a quantum
state $|\phi(\boldsymbol{x})\rangle$. The map chosen in this work is amplitude
encoding
$|\phi_{\textrm{amp}}(\boldsymbol{x})\rangle=1/||\boldsymbol{x}||_{2}\sum\limits_{j=1}^{n}x_{j}|j-1\rangle$,
where $\\{|j\rangle\\}$ are the Fock states, $n$ denotes the dimension of the
data, and $||\boldsymbol{x}||_{2}=\sqrt{\sum\limits_{j=1}^{n}x^{2}_{j}}$ is
the Euclidean norm. The universal embedding circuit $\phi$ can be associated
with a state preparation gate sequence $U_{\phi}(\boldsymbol{x})$, that
prepares an arbitrary state $|\phi(\boldsymbol{x})\rangle$ satisfying
$U_{\phi}(\boldsymbol{x})|0\rangle=|\phi(\boldsymbol{x})\rangle$ Schuld and
Killoran (2019); Lloyd  _et al._ (2020).

Machine learning algorithms frequently rely on distance measures that can
quantify the similarity between feature vectors Hastie  _et al._ (2009). Here,
we employ the Hilbert-Schmidt distance defined as
$D_{\textrm{HS}}(|\phi(\boldsymbol{x})\rangle,|\phi(\boldsymbol{x^{\prime}})\rangle)=\sqrt{2-2|\langle\phi(\boldsymbol{x})|\phi(\boldsymbol{x^{\prime}})\rangle|^{2}}$,
where $|\langle\phi(\boldsymbol{x})|\phi(\boldsymbol{x^{\prime}})\rangle|^{2}$
is the overlap between two quantum states $|\phi(\boldsymbol{x})\rangle$ and
$|\phi(\boldsymbol{x^{\prime}})\rangle$ Trávníček  _et al._ (2019). The
overlap measurement, reported throughout this work, can be efficiently
realized with a constant-depth circuit implemented by a SWAP test Filip
(2002); Gao  _et al._ (2018); Gan  _et al._ (2020); Nguyen  _et al._ (2021).

Figure 2:  $\boldsymbol{K}$-means clustering results. a. Measurement of
overlap $|\langle d_{m}|{c_{lk}}\rangle|^{2}$ between the data states and the
$k$-th centroid $|{c_{lk}}\rangle$ of the assignment step $l$. The data states
are the squeezed vacuum states ($\mathcal{H}_{\textrm{sqz}}$), superposition
of $|0\rangle$ and $|1\rangle$ Fock states ($\mathcal{H}_{\textrm{Fock}}$),
and coherent states ($\mathcal{H}_{\textrm{coh}}$). Top row shows the step-by-
step assignment of states into clusters, visualized in the $x_{3}$-$x_{4}$
scatter plot. b. The scatter plot matrix of the centroids (crosses) and data
states: $\mathcal{H}_{\textrm{sqz}}$ (circles), $\mathcal{H}_{\textrm{Fock}}$
(squares), and $\mathcal{H}_{\textrm{coh}}$ (diamonds). The classification
results of $K$-means algorithm into three clusters are labeled by different
colors: cluster 1 (red), cluster 2 (blue), and cluster 3 (purple). The
algorithm assigns $\mathcal{H}_{\textrm{sqz}}$ to cluster 1,
$\mathcal{H}_{\textrm{Fock}}$ to cluster 2, and $\mathcal{H}_{\textrm{coh}}$
to cluster 3. Inset shows the process of updating the centroids visualized in
the $x_{1}$-$x_{2}$ scatter plot. The error bar denotes the one $\sigma$ fit
error of the state characterization (see Methods).

In this work, the bosonic learning machine comprises two trapped 171Yb+ ions
in a linear radio-frequency Paul trap, employing three motional modes denoted
$A$, $B$, and $C$. We start the experiment by preparing all the modes in the
ground state of motion. Two of the modes, $B$ and $C$, are used for storing
information while $A$ acts as an auxiliary mode required for implementing the
SWAP test. The hyperfine states $|g\rangle,|e\rangle$ of one ion mediate
interaction between the motional modes, while the other ion is prepared in the
${}^{2}F_{7/2}$ metastable state and does not interact with optical fields
during the experiment (See Fig.1). The Jaynes-Cumming interaction between the
internal state and the motional modes enables deterministic preparation of an
arbitrary state with dimension $n$ in $n-1$ steps, using methods presented in
Law and Eberly (1996); Ben-Kish  _et al._ (2003). Alternatively, we use
displacement and squeezing operations to prepare coherent and squeezed vacuum
states respectively. The SWAP test relies on controlled beam splitter gates of
the form
$U_{\textrm{CBS}}(\theta,\psi)=\exp[-(i/\hbar)\int_{0}^{T}\sigma_{x}H_{\textrm{BS}}(\psi,t)dt)]$
that enacts the Hamiltonian
$H_{\textrm{BS}}(\psi)/\hbar=\frac{g(t)}{2}(a^{\dagger}be^{i\psi}+ab^{\dagger}e^{-i\psi})$,
depending on the internal state of the ions in the $\sigma_{x}$ basis Filip
(2002); Gao  _et al._ (2018); Nguyen  _et al._ (2021). Here $a(a^{\dagger})$
and $b(b^{\dagger})$ are the annihilation (creation) operator of the motional
modes, $g(t)$ is the coupling strength, $\psi$ an experimentally tunable
phase, and $\theta=\int_{0}^{T}g(t)dt$ the effective mixing angle for a gate
duration $T$ (Methods).

Given two quantum states $|\phi(\boldsymbol{x})\rangle_{B}$ and
$|\phi(\boldsymbol{x^{\prime}})\rangle_{C}$ prepared in mode $B$ and $C$
respectively, the SWAP test can be implemented with controlled beam splitter
gates $U_{\textrm{CBS}}(\pi,0)$ acting between modes $A$ and $B$, and
$U_{\textrm{CBS}}(\pi/2,\pi)$ between modes $A$ and $C$, followed by a
projection measurement of the internal state. The overlap between two states
is proportional to the probability $P_{g}$ to detect the internal state in
$|g\rangle$:
$|_{B}\langle\phi(\boldsymbol{x})|\phi(\boldsymbol{x^{\prime}})\rangle_{C}|^{2}=|1-2P_{g}|$
(Fig.1d). The SWAP test duration is independent from the dimension of the data
and is significantly shorter than the coherence times of the motional modes
Nguyen  _et al._ (2021).

To highlight applications of our bosonic learning machine, we demonstrate the
unsupervised $K$-means clustering and supervised $k$-NN algorithms using
quantum states as the input data. We note that the input data may also be
classical as it can be mapped to quantum states by using the universal
embedding circuit. We begin with the description of implementing the
unsupervised $K$-means clustering algorithm. Consider two parties Alice and
Bob. Alice prepares and sends Bob identical copies of $M$ quantum states. The
quantum states can be expressed in the Fock basis
$|d_{m}\rangle=\sum\limits_{j=0}^{n-1}x_{m,j+1}|j\rangle$, where $m$ runs from
$1$ to $M$. In this scenario, Bob is assumed to have no a priori information
about the states that Alice prepares. Without resorting to a full tomography
of the input states, Bob’s task is to sort the quantum states into $K$
clusters based on similarities among the states.

The algorithm begins with Bob making a guess of the number of clusters $K$,
and associates each cluster to a centroid state that are randomly initialized.
The states are assigned to the cluster that has the nearest centroid. In step
$l$ we implement a map $|d_{m}\rangle\rightarrow
r^{m}_{l}=\operatorname*{argmin}_{k}\\{D(|d_{m}\rangle,|c_{lk}\rangle)\\}$,
where $|c_{lk}\rangle$ denotes the centroids corresponding to the $k-$th
cluster. Our quantum advantage comes from the reduction of complexity involved
in estimating the distance of the two states $|d_{m}\rangle$ and
$|c_{lk}\rangle$. This estimation takes approximately $n$ steps to perform
classically, in comparison to the constant-time operation in our experiment.
The centroids are subsequently updated by Alice to be the mean of the states
contained in the corresponding cluster. The assignment and update are repeated
until the clusters cease to change.

Figure 3:  Pattern recognition of quantum states with $k$-nn algorithm. The
overlap between trial states and all members in each cluster were measured
with SWAP test. a. Proportion of members in each cluster in the set of $k=7$
nearest neighbors of the trial states. The trial states are assigned to the
cluster that has the largest proportion. For visualization, the calculated
Wigner functions of trial and cluster states are shown in the left column and
top row, respectively. b. Measured overlap between the training data and the
trial states: $\frac{1}{\sqrt{2}}(|0\rangle+|1\rangle)$ (top) and (up to a
normalization factor)
$\sqrt{0.3}S(r=0.8e^{i\pi})+\sqrt{0.2}(|0\rangle+|1\rangle)+\sqrt{0.5}D(\alpha=1.7)|0\rangle$
(bottom). The error bars correspond to statistical uncertainty and represent
$95\%$ confidence intervals. Other overlap measurements can be found in the
Supplementary Information.

Our data is a set of 15 quantum states and consists of three clusters: the set
$\mathcal{H}_{\textrm{sqz}}$ of squeezed vacuum states $|re^{i\pi}\rangle$,
$\mathcal{H}_{\textrm{Fock}}$ superposition of Fock states
$\cos(\varphi/2)|0\rangle+\sin(\varphi/2)|1\rangle$, and
$\mathcal{H}_{\textrm{coh}}$ coherent states $|\alpha\rangle$. We choose the
parameters as $r\in\\{1.5\,,1.1\,,1.2\,,1.3\,,1.4\\}$,
$\varphi/\pi\in\\{0.37\,,0.45\,,0.47\,,0.72\,,0.55\\}$, and
$\alpha\in\\{1.5\,,1.3\,,1.0\,,1.2\,,1.1\\}$. Our simulations verify that
$K$-means algorithm can partition this data set correctly on a classical
computer. The data states are prepared in the mode $B$. For the first step, we
choose the centroid states in mode C to be $|c_{01}\rangle=|0\rangle$,
$|c_{02}\rangle=|1\rangle$, and $|c_{03}\rangle=|2\rangle$. After each step
the centroids, truncated to dimension 5, are updated and approximately
synthesized using the embedding circuit. We characterize the state generation
through measurements of Rabi oscillations on the $|g\rangle|n+1\rangle$ to
$|e\rangle|n\rangle$ red-sideband transition (see Methods).

Figure 2a shows the overlap measurements between the data states and the
centroids. We used on average 700 experimental runs for each measurement. With
this particular choice of the initial centroids, the algorithm attains a
stationary assignment after four steps. The final classification result is
given in Fig. 2b, where the scatter plot of the amplitudes $x_{j}$ of the data
states are shown for $j$ up to 5. The algorithm associates cluster 1 to
$H_{\textrm{sqz}}$, cluster 2 to $H_{\textrm{Fock}}$, and cluster 3 to
$H_{\textrm{coh}}$. The SWAP test measurements, albeit imperfect, produce
classification results that agrees perfectly with the prior knowledge about
the states sent to Bob. Stationary assignment of the centroids is evident from
the step-by-step evolution shown in the inset of Fig. 2b. The assignment
results after each step can be found in Extended Data.

We next demonstrate the $k$-NN algorithm to classify unknown trial states, by
using the three labelled clusters from the $K$-means results as training
states. The traditional approach to solving the problem of classifying unknown
states is to perform a complete tomography, which requires a number of
measurements that grows exponentially with the state dimension $n$ Haah  _et
al._ (2017). However, if complete information regarding the state is
unnecessary, $k$-NN algorithm could offer an efficient alternative, especially
when $n$ is large. With a constant depth SWAP test, computational complexity
of the $k$-NN becomes independent of $n$.

The $k$-NN algorithm requires estimation of the overlap between each trial
state and all training states. We identify the $k$ states having highest
overlap with the trial state. The trial state is then assigned to the cluster
that labels the majority of $k$ states. Figure 3 shows the experimental
results of $k$-NN classification with $k=7$ for a variety of trial states. The
first three trial states utilize the same operation that prepares the three
clusters, respectively, and are observed to be correctly classified. We
further carry out the classification of trial states that are prepared with
other operations (See Methods). Such examples are the two states
$D(\alpha=1)(|0\rangle-|1\rangle)$ and
$\sqrt{0.3}S(r=0.8e^{i\pi})|0\rangle+\sqrt{0.2}(|0\rangle+|1\rangle)+\sqrt{0.5}D(\alpha=1.7)|0\rangle$
(up to a normalization factor), that are both classified to cluster 2. These
results indicate that the outcome of the algorithm is determined by overlap of
the quantum states, which can be visualized by the similarity of their Wigner
functions. This can be interpreted as a form of Wigner function pattern
recognition, much like applying $k$-NN to handwriting recognition in classical
ML Hastie  _et al._ (2009).

In conclusion, we realize a bosonic learning machine using motional modes in a
system of trapped ions. The core subroutine of the machine is the efficient
SWAP test that estimates the overlap of two motional modes with a constant
gate duration, a feature that is preferred in solving problems with large
dimension such as quantum support vector machines. We expect that the
experimentally accessible dimensions can be further increased by improving
various technical aspects of the setup, such as the motional coherence time
and coupling strength between modes. The same bosonic learning machine can
also be used to implement several quantum kernels proposed in Schuld and
Killoran (2019) and on different experimental platforms such as
superconducting microwave cavities Hofheinz  _et al._ (2009); Gao  _et al._
(2018). These results provide a key primitive for developing advanced quantum-
enhanced machine learning algorithms in infinite-dimensional feature spaces.

###### Acknowledgements.

We thank P.T. Le, H.N. Le, and L.C. Kwek for discussions. This research is
supported by the National Research Foundation, Prime Ministers Office,
Singapore, and the Ministry of Education, Singapore, under the Research
Centers of Excellence program and NRF Quantum Engineering Program (Award
QEP-P4)

Author contributions

C.H.N. and D.M. conceived the experiment. C.H.N. K.W.T, G.M. and H.C.J.G.
built and maintained the experimental apparatus. C.H.N. and K.W.T performed
the experiments and analyzed the experimental data with assistance from G.M.
and H.C.J.G. All authors contributed to the manuscript preparation.

Data availability

The data that support the findings of this study are available from the
corresponding authors on reasonable request.

Author Information

The authors declare no competing financial interests. Correspondence and
requests for materials should be addressed to C.H.N<EMAIL_ADDRESS>or
D.M. (phymd@nus.edu.sg).

Methods

Experimental setup. We trap two 171Yb+ ions in a four-rod linear radio-
frequency Paul trap with single-ion secular trap frequencies
$(\omega_{x},\omega_{y},\omega_{z})=2\pi\times(1.274,0.945,0.519)$ MHz. The
normal motional modes $A$, $B$, and $C$ are the out-of-phase motion in the $y$
direction, out-of-phase motion and in-phase in the $x$ direction with
frequencies $(\omega_{A},\omega_{B},\omega_{C})=2\pi\times(0.782,1.159,1.274)$
MHz, respectively. The transverse trapping frequencies are actively stabilized
and drift less than 100 Hz/hour.

One of the ions, pumped to the metastable state ${}^{2}F_{7/2}$ by driving the
${}^{2}D_{3/2}\rightarrow$ $3[1/2]^{\circ}_{3/2}$ transition at $398.98$ nm,
does not interact with optical fields and is referred to as the “dark” ion. At
the beginning of all experiments, the other ion is Doppler cooled by a 369.53
nm laser red-detuned from the $|^{2}S_{1/2},F=1\rangle$ to
$|^{2}P_{1/2},F=0\rangle$ transition, and sympathetically cools the “dark”
ion. The two hyperfine states used in the experiment are
$|g\rangle=|^{2}S_{1/2},F=0,m_{F}=0\rangle$ and
$|e\rangle=|^{2}S_{1/2},F=1,m_{F}=0\rangle$, which are separated by
$\omega_{0}/2\pi=12.643$ GHz. The ion is initialized via optical pumping in
the state $|g\rangle$. Resonance fluorescence state detection technique is
used to detect the ion’s internal state Olmschenk  _et al._ (2007).

A frequency-doubled, mode-locked Ti:Sapphire laser with a central wavelength
of 372.85 nm is used to drive the stimulated Raman transitions. The mode-
locked laser has a pulse duration of 3 ps and a repetition rate of 76.2 MHz.
Two beams, R1 and R2 (Fig. 1a,b) with an average power of 60mW each, are sent
through acoustic-optical modulators (AOMs) and focused to a beam waist of
$\sim 15\,\mu$m at the ions. The beams R1 and R2 form 45∘ and 135∘ angles with
respect to the $z$-axis (Fig. 1a). The detuning of the Raman beams are
controlled by the AOMs.

We begin all experiments with all the radial motional modes prepared in the
ground state (residual $\bar{n}\leq 0.05$) using 10 ms Sisyphus cooling
followed by 250 cycles of Raman sideband cooling.

Controlled beam splitter gates. We implement controlled beam splitter gates
$U_{\textrm{CBS}}$ between modes $A$ and $B$ by applying a bichromatic Raman
beatnote at $\omega_{0}\pm\omega_{AB}$, where
$\omega_{AB}=|\omega_{A}-\omega_{B}|$. The resulting Hamiltonian is

$H_{I}=\frac{\hbar
g(t)}{2}\sigma_{x}(a^{\dagger}be^{-i\psi}+ab^{\dagger}e^{i\psi})=\sigma_{x}H_{\textrm{BS}},$
(1)

where $\psi=\frac{\psi_{b}-\psi_{r}}{2}$, and $\psi_{r}$ and $\psi_{b}$ are
the phases of the bichromatic fields. We let
$\sigma_{x}=(e^{-i\psi_{S}}\sigma_{+}+e^{i\psi_{S}}\sigma_{-})$, where
$\psi_{\textrm{S}}=-\frac{\psi_{r}+\psi_{b}}{2}$, and
$\sigma_{+}=|e\rangle\langle g|$, $\sigma_{-}=|g\rangle\langle e|$. $g(t)$ is
the coupling strength.

The controlled beam splitter gate corresponds to the unitary evolution
$U_{\textrm{CBS}}(\theta,\psi)=\exp[-(i/\hbar)\int_{0}^{T}\sigma_{x}H_{\textrm{BS}}(\psi,t)dt)]$,
where $\theta=\int_{0}^{T}g(t)dt$. In the experiment, we obtain a coupling
strength $g/2\pi\sim 680$ Hz. This corresponds to implementing a beam splitter
operation $U_{\textrm{CBS}}(\theta=\pi/2)$ in $T\sim 365$ $\mu$s, much shorter
than the motional phase coherence time for the superposition of
$|0\rangle+|1\rangle$ which is approximately $10$ ms.

Preparation and characterization of quantum states. The universal embedding
circuit relies on a series $n-1$ pairs of carrier and red-sideband:

$\prod_{l=1}^{n-1}U_{\textrm{rsb}}(\tau_{l},\Upsilon_{l})U_{c}(t_{l},\upsilon_{l}+\upsilon^{{}^{\prime}}_{l})$
(2)

where $t_{l}$ ($\tau_{l}$) and $\upsilon_{l}$ ($\Upsilon_{l}$) are the
duration and phase for the $l$-th carrier (red-sideband) pulse that are
determined by the Eberly-Law algorithm Law and Eberly (1996). The additional
phases $\upsilon^{{}^{\prime}}_{l}$ compensate for the rotation of the
internal states due to the AC Stark shifts of the red-sideband pulses (See
Supplementary). The embedding circuit is used to prepare the centroid states
of $K$-means, the cluster state $\mathcal{H}_{\textrm{Fock}}$ of $K$-means
experiment, and the trial state (up to a normalization factor)
$\sqrt{0.3}S(r=0.8e^{i\pi})|0\rangle+\sqrt{0.2}(|0\rangle+|1\rangle)+\sqrt{0.5}D(\alpha=1.7)|0\rangle$
of $k$-NN experiment. We use spin-dependent force to prepare all coherent
states and the squeezed vacuum states Leibfried  _et al._ (2003). For the
trial state $D(\alpha=1)(|0\rangle-|1\rangle)$, we first prepare
$|0\rangle-|1\rangle$, followed by a displacement of $\alpha=1$.

The real and positive coefficients $x_{j+1}$ of the data states
$|d_{m}\rangle$, shown in Fig. 2b, are the square root of population
distribution $p(j)$, where $p(j)=|\langle d_{m}|j\rangle|^{2}$. The population
distribution is determined from a Fourier transform of the time evolution
while driving the red-sideband transition

$P_{e}(\tau)=\sum\limits_{j=1}^{\infty}p(j)e^{-\sqrt{j}\gamma\tau}(1-\cos{(\sqrt{j}\Omega_{\textrm{rsb}}\tau}))/2$
(3)

where $P_{e}$ is the probability to detect the ion in the state $|e\rangle$,
$\Omega_{\textrm{rsb}}$ is the Rabi frequency of the red-sideband, and
$\gamma$ is the decoherence rate of motional states. To determine $p(j)$ a fit
constrained to be non-negative and $\sum_{j=0}^{6}p(j)=1$ is applied to Eq. 3
for $j=0$ to $6$.

## References

  * Biamonte  _et al._ (2017) J. Biamonte, P. Wittek, N. Pancotti, P. Rebentrost, N. Wiebe, and S. Lloyd, Nature 549, 195 (2017).
  * Cai  _et al._ (2015) X.-D. Cai, D. Wu, Z.-E. Su, M.-C. Chen, X.-L. Wang, L. Li, N.-L. Liu, C.-Y. Lu, and J.-W. Pan, Physical Review Letters 114, 110504 (2015).
  * Trávníček  _et al._ (2019) V. Trávníček, K. Bartkiewicz, A. Černoch, and K. Lemr, Physical Review Letters 123, 260501 (2019).
  * Havlíček  _et al._ (2019) V. Havlíček, A. D. Córcoles, K. Temme, A. W. Harrow, A. Kandala, J. M. Chow, and J. M. Gambetta, Nature 567, 209 (2019).
  * Johri  _et al._ (2020) S. Johri, S. Debnath, A. Mocherla, A. Singh, A. Prakash, J. Kim, and I. Kerenidis,  (2020), arXiv:2012.04145 .
  * Bartkiewicz  _et al._ (2020) K. Bartkiewicz, C. Gneiting, A. Černoch, K. Jiráková, K. Lemr, and F. Nori, Scientific Reports 10, 12356 (2020).
  * Xin  _et al._ (2021) T. Xin, L. Che, C. Xi, A. Singh, X. Nie, J. Li, Y. Dong, and D. Lu, Phys. Rev. Lett. 126, 110502 (2021).
  * Peters  _et al._ (2021) E. Peters, J. Caldeira, A. Ho, S. Leichenauer, M. Mohseni, H. Neven, P. Spentzouris, D. Strain, and G. N. Perdue, (2021), arXiv:2101.09581 [quant-ph] .
  * Braunstein and van Loock (2005) S. L. Braunstein and P. van Loock, Reviews of Modern Physics 77, 513 (2005).
  * Ortiz-Gutiérrez  _et al._ (2017) L. Ortiz-Gutiérrez, B. Gabrielly, L. F. Muñoz, K. T. Pereira, J. G. Filgueiras, and A. S. Villar, Optics Communications 397, 166 (2017).
  * Flühmann  _et al._ (2019) C. Flühmann, T. L. Nguyen, M. Marinelli, V. Negnevitsky, K. Mehta, and J. P. Home, Nature 566, 513 (2019).
  * Lloyd  _et al._ (2013) S. Lloyd, M. Mohseni, and P. Rebentrost,  (2013), arXiv:1307.0411 .
  * Lloyd  _et al._ (2014) S. Lloyd, M. Mohseni, and P. Rebentrost, Nature Physics 10, 631 (2014).
  * Lloyd  _et al._ (2016) S. Lloyd, S. Garnerone, and P. Zanardi, Nature Communications 7, 10138 (2016).
  * Rebentrost  _et al._ (2014) P. Rebentrost, M. Mohseni, and S. Lloyd, Physical Review Letters 113, 130503 (2014).
  * Lau  _et al._ (2017) H.-K. Lau, R. Pooser, G. Siopsis, and C. Weedbrook, Physical Review Letters 118, 080501 (2017).
  * Schuld and Killoran (2019) M. Schuld and N. Killoran, Physical Review Letters 122, 040504 (2019).
  * Lloyd  _et al._ (2020) S. Lloyd, M. Schuld, A. Ijaz, J. Izaac, and N. Killoran,  (2020), arXiv:2001.03622 .
  * Hastie  _et al._ (2009) T. Hastie, R. Tibshirani, and J. H. Friedman,  _The Elements of Statistical Learning: Data Mining, Inference, and Prediction_ , Springer series in statistics (Springer, 2009).
  * Filip (2002) R. Filip, Physical Review A 65, 062320 (2002).
  * Gao  _et al._ (2018) Y. Y. Gao, B. J. Lester, Y. Zhang, C. Wang, S. Rosenblum, L. Frunzio, L. Jiang, S. M. Girvin, and R. J. Schoelkopf, Physical Review X 8, 021073 (2018).
  * Gan  _et al._ (2020) H. C. J. Gan, G. Maslennikov, K.-W. Tseng, C. Nguyen, and D. Matsukevich, Phys. Rev. Lett. 124, 170502 (2020).
  * Nguyen  _et al._ (2021) C.-H. Nguyen, K.-W. Tseng, G. Maslennikov, H. C. J. Gan, and D. Matsukevich, (2021), arXiv:2103.10219 [quant-ph] .
  * Law and Eberly (1996) C. K. Law and J. H. Eberly, Physical Review Letters 76, 1055 (1996).
  * Ben-Kish  _et al._ (2003) A. Ben-Kish, B. DeMarco, V. Meyer, M. Rowe, J. Britton, W. M. Itano, B. M. Jelenković, C. Langer, D. Leibfried, T. Rosenband, and D. J. Wineland, Phys. Rev. Lett. 90, 037902 (2003).
  * Haah  _et al._ (2017) J. Haah, A. W. Harrow, Z. Ji, X. Wu, and N. Yu, IEEE Transactions on Information Theory 63, 5628 (2017).
  * Hofheinz  _et al._ (2009) M. Hofheinz, H. Wang, M. Ansmann, R. C. Bialczak, E. Lucero, M. Neeley, A. D. O’Connell, D. Sank, J. Wenner, J. M. Martinis, and A. N. Cleland, Nature 459, 546 (2009).
  * Olmschenk  _et al._ (2007) S. Olmschenk, K. C. Younge, D. L. Moehring, D. N. Matsukevich, P. Maunz, and C. Monroe, Physical Review A 76, 052314 (2007).
  * Leibfried  _et al._ (2003) D. Leibfried, R. Blatt, C. Monroe, and D. Wineland, Reviews of Modern Physics 75, 281 (2003).

Supplementary Material: Experimental quantum-enhanced machine learning over
infinite dimensions

##  I Delay time between state preparation and the SWAP test

The two steps required for operation of the bosonic learning machine are the
state preparation and the overlap measurement using the SWAP test. However,
the quantum state of a harmonic oscillator evolves in time as
$|\Psi\rangle=\sum_{j}c_{j}e^{-i\omega jt}|j\rangle$ where $\omega$ is the
motional frequency. If the frequencies of motional modes $\omega_{B}$ and
$\omega_{C}$ are different, care must be taken to calibrate waiting time
between preparation and measurement to ensure that intended states are passed
to the SWAP test. We set the delay time as $nT+\textrm{offset}$. Otherwise,
for example, the overlap of the two coherent states in modes $B$ and $C$ will
oscillate as a function of the delay between the two steps, with a period of
$T=2\pi/|\omega_{C}-\omega_{B}|$.

The trap frequency measurements done by scanning the motional sidebands are
affected by the Stark shifts of the internal states. To obtain a frequency
difference between the modes which is free from AC Stark shifts we measure the
overlap $|\langle\Psi|\Phi_{1}\rangle|^{2}$ as a function of the delay time
$\tau$, where
$|\Psi\rangle=|\Phi_{1}\rangle=1/\sqrt{4}(|0\rangle+|1\rangle+|2\rangle+|3\rangle)$.
Figure S1 shows the results of overlap measurement. The period $T$ is obtained
from the time difference between two peaks in the
$|\langle\Psi|\Phi_{1}\rangle|^{2}$ overlap measurement. We estimate $T\sim
8.2\,\mu$s.

Figure S1:  Measurement of $|\langle\Psi|\Phi_{1}\rangle|^{2}$ (blue circles)
and $|\langle\Psi|\Phi_{2}\rangle|^{2}$ (orange squares) as a function of
waiting time $\tau$, where
$|\Psi\rangle=|\Phi_{1}\rangle=1/\sqrt{4}(|0\rangle+|1\rangle+|2\rangle+|3\rangle)$
and $|\Phi_{2}\rangle=1/\sqrt{4}(|0\rangle-|1\rangle+|2\rangle-|3\rangle)$.

##  II Additional details on state preparation

The Fourier analysis of the time evolution while driving the red sideband can
measure only the diagonal terms (phonon distribution) of the density matrix of
the states. To verify that the universal embedding circuit can prepare states
with desired off-diagonal terms, we carry out the following experiment. We
prepare $|\Phi_{2}\rangle=1/\sqrt{4}(|0\rangle-|1\rangle+|2\rangle-|3\rangle)$
and compare the two overlaps $|\langle\Psi|\Phi_{1}\rangle|^{2}$ and
$|\langle\Psi|\Phi_{2}\rangle|^{2}$ as a function of waiting time $\tau$.
Simulations predict that the second overlap
$|\langle\Psi|\Phi_{2}\rangle|^{2}$ is shifted in time, relative to
$|\langle\Psi|\Phi_{1}\rangle|^{2}$, by $T/2$ if the states with correct off-
diagonal terms are prepared. The experimental results shown in Fig. S1 agree
with calculations.

Figure S2:  a. Circuit diagram for measuring the AC Stark shift induced by the
inserted rsb pulse $U_{\textrm{rsb}}$. $R(\pi/2,\phi)$: $\pi/2$ rotation
pulses with phase $\phi$ of the internal states. b. Probability of detecting
the internal state in $|e\rangle$ as a function of the red-sideband pulse
duration $\tau$ with $\phi=0$ (blue circles) and $\phi=-2\pi\times 5.9\tau$
kHz (orange squares). The dashed line is the calculated probability when the
AC Stark is compensated.

The red-sideband pulse $U_{rsb}(\tau_{l},\Upsilon_{l})$ induces an AC Stark
shift $\hbar\Delta\omega_{\textrm{AC}}$ that gives rise to a rotation
$U=e^{-i\Delta\omega_{\textrm{AC}}\tau_{l}\sigma_{z}}$ of the internal states.
To compensate for this rotation, an additional phase
$\upsilon^{{}^{\prime}}_{l+1}=\sum_{l^{\prime}\leq
l}\Delta\omega_{\textrm{AC}}\tau_{l^{\prime}}$ is included to the next carrier
pulses $U_{c}(t_{l},\upsilon_{l}+\upsilon^{{}^{\prime}}_{l})$. We determine
$\Delta\omega_{\textrm{AC}}$ using the Ramsey experiment with an inserted red-
sideband pulse described in Figure S2a. The second $\pi/2$ rotation pulse has
a phase $\phi$. We vary the red-sideband pulse duration $\tau$ while
monitoring the probability of detecting the internal state in $|e\rangle$.
When $\phi$ is set to $-\Delta\omega_{\textrm{AC}}\tau$ and
$\Delta\omega_{\textrm{AC}}\sim 2\pi\times 5.9\,$KHz, the measurement result
(Fig. S2b) agrees with the expectation when the AC Stark shift is compensated.

##  III Extended Data

Figure S3:  Step-by-step cluster assignment of $K$-mean algorithm visualized
in scatter plot matrix of the centroids (crosses) and data states:
$\mathcal{H}_{\textrm{sqz}}$ (circles), $\mathcal{H}_{\textrm{Fock}}$
(squares), and $\mathcal{H}_{\textrm{coh}}$ (diamonds). The error bar denotes
the one $\sigma$ fit error of the state characterization.

Figure S4:  The overlap between trial states and all members in each cluster
in the $k$-nn experiment. The error bars correspond to statistical uncertainty
and represent $95\%$ confidence intervals.

Figure S5:  Population distribution of the centroid states in the Fock basis.
Blue sideband pulses are used to analyze centroid states of step 0 and red
sideband pulses to analyze the rest of the centroid states. Insets show the
probability to detect the internal state in $|e\rangle$ as a function of the
pulse duration. The red curves are the fits to the population using Eq. 3. The
orange bars show the desired population distribution.