# Powering an autonomous clock with quantum electromechanics

Oisín Culhane<EMAIL_ADDRESS>School of Physics, Trinity College Dublin,
Dublin 2, Ireland  Michael J. Kewming<EMAIL_ADDRESS>School of Physics,
Trinity College Dublin, Dublin 2, Ireland  Alessandro Silva
<EMAIL_ADDRESS>SISSA, Via Bonomea 265, I-34135 Trieste, Italy  John
Goold<EMAIL_ADDRESS>School of Physics, Trinity College Dublin, Dublin 2,
Ireland  Trinity Quantum Alliance, Unit 16, Trinity Technology and Enterprise
Centre, Pearse Street, Dublin 2, D02YN67  Mark T. Mitchison
<EMAIL_ADDRESS>School of Physics, Trinity College Dublin, Dublin 2,
Ireland  Trinity Quantum Alliance, Unit 16, Trinity Technology and Enterprise
Centre, Pearse Street, Dublin 2, D02YN67

###### Abstract

We theoretically analyse an autonomous clock comprising a
nanoelectromechanical system, which undergoes self-oscillations driven by
electron tunnelling. The periodic mechanical motion behaves as the clockwork,
similar to the swinging of a pendulum, while induced oscillations in the
electrical current can be used to read out the ticks. We simulate the dynamics
of the system in the quasi-adiabatic limit of slow mechanical motion, allowing
us to infer statistical properties of the clock’s ticks from the current auto-
correlation function. The distribution of individual ticks exhibits a tradeoff
between accuracy, resolution, and dissipation, as expected from previous
literature. Going beyond the distribution of individual ticks, we investigate
how clock accuracy varies over different integration times by computing the
Allan variance. We observe non-monotonic features in the Allan variance as a
function of time and applied voltage, which can be explained by the presence
of temporal correlations between ticks. These correlations are shown to yield
a precision advantage for timekeeping over the timescales that the
correlations persist. Our results illustrate the non-trivial features of the
tick series produced by nanoscale clocks, and pave the way for experimental
investigation of clock thermodynamics using nanoelectromechanical systems.

##  I Introduction

The measurement of time is perhaps the most primitive and also the most
important of all physical measurements. Not only are clocks ubiquitous in our
daily lives, but they also underpin satellite navigation, magnetometry, and
gravimetry [1, 2], among many other applications in metrology. However,
despite the spectacular precision that is now routinely achieved with atomic
clocks, fundamental gaps in our understanding of time at the quantum level
still remain.

One fruitful and long-standing [3, 4, 5] strategy to address this problem is
to formulate and study models of clocks as dynamical systems. Recent work in
this direction has shown that basic physical considerations, stemming from
information theory on the one hand [6, 7, 8, 9, 10] and thermodynamics on the
other [11, 12, 13, 14], strongly constrain both the internal structure and the
achievable performance of timekeeping devices. For example, it is now
understood that accurate clocks must produce a lot of entropy [11, 13],
quantum clocks can be more accurate than classical clocks with the same number
of accessible states [10, 9, 15], and that there exists a tradeoff between a
clock’s accuracy and its timing resolution [16]. Similar constraints on
temporal precision known as thermodynamic or kinetic uncertainty relations
have been discovered in other, seemingly unrelated contexts, such as
biomolecular systems [17, 18, 19, 20] or mesoscopic transport setups [21, 22,
23]. Beyond their foundational importance, moreover, the physical limits of
timekeeping have consequences for the fidelity of time-dependent control
protocols arising in quantum thermodynamics [24, 25, 26] and quantum
computation [27, 28], among other applications.

Here we focus on ticking clocks, which produce a continuous time reference to
an external observer, as opposed to stopwatch-like systems (e.g. the Larmor
clock [29]) that measure only a single time interval. Any quantum ticking
clock comprises two distinct parts: the clockwork and the register [7]. The
clockwork is a non-equilibrium system whose dynamics generates the ticks,
which are then recorded by the register. Separation of clockwork and register
ensures that the tick readout does not disrupt the clockwork [30]. Although
almost any out-of-equilibrium system can act as a clockwork, accurate
timekeeping requires temporal probability concentration [14], whereby the
irreversible process that generates each tick is permitted to occur only at
well-spaced time intervals. For example, while a simple hourglass exploits
freefalling mass (e.g. sand) to indicate the passage of time, the escapement
of a pendulum clock allows a mass to fall incrementally — thereby generating
the next tick — only after one full period of the swinging pendulum has
elapsed [13, 31].

A natural form of autonomous clockwork is a self-oscillator [32], i.e. a
system that oscillates in the absence of periodic forcing.
Nanoelectromechanical systems [33] are an interesting and experimentally
relevant class of self-oscillators, in which the periodic motion of a
mechanical degree of freedom is driven by an electron tunnelling. This can
occur either by electron shuttling [34, 35, 36, 37, 38, 39], where single-
electron tunnelling resonantly drives the oscillator, or via a subtler
electromechanical instability that emerges from the time-delayed effect of
many non-resonant tunnelling events [40, 41, 42, 43, 44, 45, 46].
Thermodynamic aspects of both kinds of self-oscillations have recently been
investigated [47, 48, 49].

In this work, we consider an autonomous clock powered by electromechanical
self-oscillations of the latter kind, motivated by the experiments reported in
Refs. [45, 46]. In our model, periodic mechanical motion provides the
clockwork while the ticks are registered by monitoring the induced
oscillations in the electronic current. As in previous studies, the autonomous
nature of the clockwork dynamics allows for a self-contained analysis of the
thermodynamic resources needed to measure time [11]. However, our work goes
beyond existing literature on clock thermodynamics in at least two important
respects.

First, our electromechanical clockwork is based on the oscillations of a
continuous degree of freedom, reminiscent of a pendulum. Yet, unlike a
pendulum clock, or indeed any other clock considered in the literature to our
knowledge, these oscillations are driven by measurement backaction associated
with the tick readout itself. By contrast, existing theoretical models of
quantum clocks consider the tick generation and readout to arise from distinct
irreversible processes [11, 14, 10, 9, 15, 16, 30], while experimental work on
the thermodynamics of timekeeping has either lain in a purely classical regime
where measurement backaction is negligible [50] or focussed purely on the role
of measurement backaction in a non-autonomous setting [51]. Our setup is also
distinct from Ref. [52], which considered a quantum clock powered by
measurements that are unrelated to tick readout.

Second, we go beyond the performance measures introduced in previous works
[11, 10] by evaluating the Allan variance [53], which quantifies clock
accuracy over different integration times. At intermediate timescales, our
results corroborate the general expectation that higher accuracy entails more
dissipation. However, at longer timescales this trend reverses. To explain
this crossover in the Allan variance, we show that the non-trivial dynamics of
the clockwork generates persistent temporal correlations between the clock’s
ticks. Our results thus call for further research to understand both the
origin and consequences of such temporal correlations, which have been largely
neglected in previous studies of autonomous quantum clocks.

Figure 1: Schematic of the setup. A central system, here a quantum dot, is
coupled to two electrodes each at a temperature $T$, with the left electrode
having a chemical potential of $\mu_{L}$ and the right a chemical potential of
$\mu_{R}$. The quantum dot is coupled to simple harmonic oscillator mode.
Under an appropriate static voltage bias, the current flowing through the
system drives self-sustained mechanical oscillations, which in turn induces an
oscillatory component in the current. A series of ticks are recorded by
counting oscillations in the output current.

An outline of the manuscript is as follows. We first introduce our model for
an electromechanical clock (Sec. II.1), and explain how we model the
clockwork’s nonequilibrium dynamics (Sec. II.2) and the readout of ticks from
oscillations induced in the electrical current (Sec. II.3). Our modelling is
based on the quasi-adiabatic assumption, i.e. that the mechanical motion is
slow compared to electronic relaxation, which is well satisfied in current
experiments [45, 46]. Then we examine the induced current oscillations in
detail, identifying an interesting separation of timescales between amplitude
and phase fluctuations (Sec. III.1). Next, we analyse the distribution of
individual ticks of the clock (Sec. III.2) and find a tradeoff between
accuracy, resolution, and entropy production, as expected from previous work
[11]. Finally, we show that the Allan variance exhibits a non-trivial
dependence on integration time and applied voltage (Sec. III.3), which can be
understood in terms of correlations between the clock’s ticks that we quantify
explicitly (Sec. III.4). We conclude with a summary and a brief outlook on
future work (Sec. IV). The reduced Planck constant, the Boltzmann constant,
and the electron charge are all set to unity throughout this work.

##  II Model

###  II.1 Description of the setup

We consider an electromechanical system similar to the model outlined in [49].
A schematic of the setup is shown in Fig. 1. The clockwork consists of a
resonant-level model coupled to a harmonic oscillator. The Hamiltonian of the
clockwork is given by a sum of the following terms:

$\displaystyle\hat{H}_{S}=\epsilon\hat{c}^{\dagger}\hat{c},$ (1)
$\displaystyle\hat{H}_{B}=\sum_{k}\left(\Omega_{kL}\hat{d}^{\dagger}_{kL}\hat{d}_{kL}+\Omega_{kR}\hat{d}^{\dagger}_{kR}\hat{d}_{kR}\right),$
(2)
$\displaystyle\hat{H}_{T}=\sum_{k}\left(g_{kL}\left[\hat{c}^{\dagger}\hat{d}_{kL}+\hat{d}^{\dagger}_{kL}\hat{c}\right]+g_{kR}\left[\hat{c}^{\dagger}\hat{d}_{kR}+\hat{d}^{\dagger}_{kR}\hat{c}\right]\right),$
(3)
$\displaystyle\hat{H}_{V}=\frac{\hat{p}^{2}}{2m}+\frac{m\omega_{0}^{2}\hat{x}^{2}}{2}-F\hat{n}\hat{x}.$
(4)

Here, $\hat{H}_{S}$ is the Hamiltonian of a single-level fermionic quantum
dot, described by its energy $\epsilon$ and the annihilation operator
$\hat{c}$. The leads coupled to the quantum dot are defined by the Hamiltonian
$\hat{H}_{B}$. Each lead is described by an infinite sum of noninteracting
Fermi modes with energy $\Omega_{k\alpha}$ and annihilation operator
$\hat{d}_{k\alpha}$, where $\alpha=L,R$ denotes the left or right lead. The
tunnelling of electrons between the dot and the leads is described by
$\hat{H}_{T}$, with $g_{k\alpha}$ the coupling strength. The vibrational
Hamiltonian $\hat{H}_{V}$ describes a single-mode mechanical oscillator with
momentum operator $\hat{p}$, position operator $\hat{x}$, natural oscillation
frequency $\omega_{0}$, and mass $m$. The final term in $\hat{H}_{V}$
describes the electromechanical interaction between the quantum dot and the
harmonic oscillator, where $F$ is the force per unit charge and $\hat{n}$ is
the excess charge present on the quantum dot, i.e.
$\hat{n}=\hat{c}^{\dagger}\hat{c}-N_{0}$, with
$N_{0}=\langle\hat{c}^{\dagger}\hat{c}\rangle_{F=0}$ being the average charge
on the oscillator in the absence of the electromechanical coupling.

The leads are initially at equilibrium with equal inverse temperature $\beta$
but different chemical potentials, $\mu_{L}\neq\mu_{R}$, with
$V=\mu_{L}-\mu_{R}$ denoting the voltage bias. The tunnelling rate for each
lead is determined by its spectral density

$\kappa_{\alpha}(E)=2\pi\sum_{k}g_{k\alpha}^{2}\delta(E-\Omega_{k\alpha}).$
(5)

Self-sustained oscillations require an energy-dependent spectral density [42].
As in our previous work [49], we assume the simple Lorentzian form

$\kappa_{\alpha}(E)=\frac{\Gamma\delta^{2}}{(E-\omega_{\alpha})^{2}+\delta^{2}},$
(6)

which models a single band of states with bandwidth $\delta$ centred on energy
$\omega_{\alpha}$, with $\Gamma$ representing the overall tunnelling rate. We
focus on the quasi-adiabatic limit, where the relaxation rate of the
electronic degrees of freedom are on a much faster timescale than the
characteristic timescale of the oscillator. This is tantamount to the
conditions

$\min(\beta^{-1},V),\Gamma\gg\omega_{0},\lambda,$ (7)

where $\lambda=F/\sqrt{m\omega_{0}}$ is the effective electromechanical
coupling strength. These conditions are well satisfied in typical experiments
[45, 46] and also by our specific parameter choices below. Note that our
minimal model neglects charging effects due to Coulomb interaction as well as
intrinsic mechanical damping of the oscillator, since neither of these effects
are crucial for the functioning of the clock.

###  II.2 Langevin equation for the clockwork dynamics

When the quasi-adiabatic conditions (7) hold, the electronic degrees of
freedom can be adiabatically eliminated to obtain a description for the
oscillator alone. Full details can be found in the Supplemental Material of
Ref. [49], but we briefly sketch the derivation here because it is serves as a
precursor for the modelling of tick readout described in Sec. II.3.

We divide the Hamiltonian into two parts as
$\hat{H}=\hat{H}_{\text{el}}+\hat{H}_{V}$, where
$\hat{H}_{\text{el}}=\hat{H}_{S}+\hat{H}_{B}+\hat{H}_{T}$ describes the
electronic degrees of freedom and $\hat{H}_{V}$ describes the oscillator
dynamics including the electromechanical coupling. We now introduce the Wigner
transform of the density matrix

$\hat{\rho}(x,p)=\frac{1}{\pi}\int\text{d}y\bra{x+y}\hat{\rho}\ket{x+y}\text{e}^{-i2yp},$
(8)

where $\ket{x}$ is a position eigenstate for the oscillator, so that
$\hat{\rho}(x,p)$ acts as a quantum operator on the electronic Hilbert space.
The Wigner function of the oscillator follows from
$W(x,p)=\text{tr}[\hat{\rho}(x,p)]$, while the electronic state is given by
$\hat{\rho}_{\text{el}}=\int\text{d}x\int\text{d}p\hat{\rho}(x,p)$. Using the
operator correspondences for the Wigner function [54], one finds that the
quantum Liouville equation
$\text{d}\hat{\rho}/\text{d}t=-i[\hat{H},\hat{\rho}]$ is equivalent to

$\frac{\partial}{\partial
t}\hat{\rho}(x,p)=\left(\mathcal{L}_{x}+\mathcal{H}+\mathcal{F}\right)\hat{\rho}(x,p),$
(9)

where we define the following superoperators:

$\displaystyle\mathcal{L}_{x}\bullet$
$\displaystyle=-i\left[\hat{H}_{\text{el}}-Fx\hat{n},\bullet\right],$ (10)
$\displaystyle\mathcal{H}\bullet$
$\displaystyle=\left(m\omega_{0}^{2}x-F\langle\hat{n}\rangle_{x}\right)\frac{\partial}{\partial
p}\bullet-\frac{p}{m}\frac{\partial}{\partial x}\bullet,$ (11)
$\displaystyle\mathcal{F}\bullet$
$\displaystyle=-\frac{F}{2}\frac{\partial}{\partial
p}\left\\{\hat{n}-\langle\hat{n}\rangle_{x},\bullet\right\\}.$ (12)

Here, $\langle\bullet\rangle_{x}=\text{tr}\left[\bullet\hat{\rho}_{x}\right]$
denotes an average with respect to the electronic stationary state,
$\hat{\rho}_{x}$, which depends on $x$ but not on $p$, and satisfies
$\mathcal{L}_{x}\hat{\rho}_{x}=0$ and
$\text{tr}\left[\hat{\rho}_{x}\right]=1$. Physically, $\mathcal{L}_{x}$
represents the evolution of the electronic degrees of freedom conditioned on a
particular oscillator position $x$, $\mathcal{H}$ represents the free
mechanical evolution including the mean-field electromechanical force, and
$\mathcal{F}$ represents the fluctuating force on the oscillator due to
deviations of the electronic charge from its steady-state mean value.

By assumption (i.e. because of the conditions (7)), we have
$\mathcal{L}_{x}\gg\mathcal{H},\mathcal{F}$ so that the electronic degrees of
freedom are only weakly perturbed from their instantaneous steady state
$\hat{\rho}_{x}$. This is formalised by defining, for any operator-valued
function $\hat{A}(x,p)$ of the oscillator’s phase-space coordinates, the
projector

$\mathcal{P}\hat{A}(x,p)=\text{tr}\left[\hat{A}(x,p)\right]\hat{\rho}_{x},$
(13)

whose orthogonal complement is denoted by $\mathcal{Q}=1-\mathcal{P}$. The
“relevant” part of the density matrix $\mathcal{P}\hat{\rho}(x,p)$ determines
the Wigner function, since $W(x,p)=\tr\mathcal{P}\hat{\rho}(x,p)$. Thus,
following the standard procedure [55], we can perturbatively eliminate the
irrelevant part, $\mathcal{Q}\hat{\rho}(x,p)$, taking into account
contributions up to second order in $\mathcal{F}$ and $\mathcal{H}$ within a
Born-Markov approximation (see Ref. [49] for details). This results in a
Fokker-Planck equation for the Wigner function, which is equivalent to the
classical Langevin equation [42, 43, 44]

$m\ddot{x}+m\gamma_{x}\dot{x}+m\omega_{0}^{2}x=F\langle\hat{n}\rangle_{x}+\sqrt{D_{x}}\xi(t),$
(14)

where $\xi(t)$ is a white-noise force with zero mean, $\mathbb{E}[\xi(t)]=0$,
and autocorrelation function
$\mathbb{E}[\xi(t)\xi(t^{\prime})]=\delta(t-t^{\prime})$, while $\gamma_{x}$
and $D_{x}$ are the damping and diffusion terms respectively, given by

$D_{x}=S_{x}(0),\qquad
m\gamma_{x}=\frac{\text{d}S_{x}(\omega)}{\text{d}\omega}\bigg{\rvert}_{\omega=0},$
(15)

with $S_{x}(\omega)$ the noise spectrum of the fluctuating force:

$S_{x}(\omega)=F^{2}\int^{\infty}_{-\infty}\text{d}t\,e^{i\omega
t}\left[\langle\hat{n}(t)\hat{n}(0)\rangle_{x}-\langle\hat{n}\rangle_{x}^{2}\right].$
(16)

In particular, since the oscillator’s Wigner function is positive within the
quasi-adiabatic approximation, it can be reconstructed by stochastically
sampling the solutions of Eq. (14).

The dynamics of the clockwork depends qualitatively on the voltage bias. For
zero or very low voltage bias, the leads behave effectively like an
equilibrium bath for the oscillator, such that the fluctuation-dissipation
relation $\beta D_{x}=2m\gamma_{x}$ is approximately satisfied. The Langevin
equation (14) then predicts a steady state for the oscillator that is almost
thermal. Conversely, when the voltage exceeds a certain threshold value,
electron tunnelling drives self-oscillations of the mechanical motion. In this
regime, the damping rate $\gamma_{x}$ entering the Langevin equation (14)
becomes negative for small $x$. As a result, the long-time state of the
oscillator predicted by Eq. (14) is a limit cycle that oscillates with period
$\sim 2\pi/\omega_{0}$ [45, 49].

###  II.3 Electrical current and current noise

Under a voltage bias, electrical current flows continuously through the
system. Mechanical oscillations induce variations of this current that can be
monitored to generate a series of ticks, as depicted in Fig. 1 and explained
in more detail in Sec. III.1. A key assumption of our work is that the current
signal is effectively classical, i.e. it can be read out without additional
measurement backaction beyond that which is already contained in the damping
and diffusion terms of Eq. (14). This is justified in the quasi-adiabatic
limit, because a very large number of electrons tunnel through the device
within a given oscillation period. In practice, oscillations in the current
can be continuously monitored by using a resonant tank circuit to filter
frequency components near the relevant frequency, together with quantum-
limited amplification [45, 56].

To find how the oscillator motion modifies the current, we use the fact that
$\hat{\rho}(x,p)\approx\mathcal{P}\hat{\rho}(x,p)$ to lowest order in small
quantities [49]. The electronic state at long times is therefore given to
lowest order by

$\hat{\rho}_{\text{el}}\approx\int\text{d}x\int\text{d}p\,W(x,p)\hat{\rho}_{x}=\int\text{d}x\,P(x)\hat{\rho}_{x}.$
(17)

where $W(x,p)$ is the steady-state Wigner function and
$P(x)=\int\text{d}pW(x,p)$ is the corresponding marginal probability
distribution for the oscillator’s position. Equation (17) can be used to
compute observables of the electronic system. In particular, the average of
the electronic current operator, $\hat{I}$, is given to lowest order by

$\langle\hat{I}\rangle=\text{tr}\left[\hat{I}\hat{\rho}_{\text{el}}\right]\approx\int\text{d}x\,P(x)\text{tr}\left[\hat{I}\hat{\rho}_{x}\right]=\int\text{d}x\,P(x)\braket{\hat{I}}_{x}.$
(18)

Here we have defined
$\braket{\hat{I}}_{x}\equiv\tr\left[\hat{I}\hat{\rho}_{x}\right]$ as the
conditional expectation value of the current operator, given the oscillator’s
position $x$.

While the ensemble-averaged current tends to the stationary value (18) at long
times, individual trajectories of the measured current will show periodic
oscillations. This is because fluctuations of the oscillator position are
transduced into current fluctuations, as characterised by the symmetrised
autocorrelation function [57]

$\displaystyle C(t)$
$\displaystyle=\tfrac{1}{2}\langle\\{\hat{I}(t),\hat{I}(0)\\}\rangle-\langle\hat{I}\rangle^{2}$
$\displaystyle=\real\left(\tr\left[\text{e}^{i\hat{H}t}\hat{I}\text{e}^{-i\hat{H}t}\hat{I}\hat{\rho}\right]\right)-\braket{\hat{I}}^{2},$
(19)

where all averages are taken with respect to the steady-state density matrix,
$\hat{\rho}$. Let us define an auxiliary density operator

$\hat{K}(t)=\text{e}^{-i\hat{H}t}\hat{I}\hat{\rho}\text{e}^{i\hat{H}t},$ (20)

such that the correlation function can be found from the real part of

$\langle\hat{I}(t)\hat{I}(0)\rangle=\tr\left[\hat{I}\hat{K}(t)\right]=\int\text{d}x\int\text{d}p\tr\left[\hat{I}\hat{K}(x,p,t)\right],$
(21)

where $\hat{K}(x,p,t)$ is the Wigner transform of of $\hat{K}(t)$, defined
analogously to Eq. (8). Since $\hat{K}(t)$ obeys the Liouville-von Neumann
equation

$\frac{\text{d}\hat{K}}{\text{d}t}=-i\left[\hat{H},\hat{K}\right],$ (22)

we have that $\hat{K}(x,p,t)$ obeys the same equation of motion as
$\hat{\rho}(x,p))$: namely,

$\frac{\partial}{\partial
t}\hat{K}(x,p,t)=(\mathcal{L}_{x}+\mathcal{H}+\mathcal{F})\hat{K}(x,p,t),$
(23)

where the superoperators on the right-hand side are defined in Eqs. (10)–(12).
We can therefore approximate the solution to this equation following an
adiabatic elimination procedure analogous to Sec. II.2, but where now both
$\mathcal{P}$ and $\mathcal{Q}$ projections are needed to reconstruct the
correlation function.

As shown in Appendix A, this procedure yields the current autocorrelation
function $C(t)=C_{\mathcal{Q}}(t)+C_{\mathcal{P}}(t)$, with

$\displaystyle C_{\mathcal{Q}}(t)=\mathbb{E}[\Delta_{x(0)}]\delta(t),$ (24)
$\displaystyle
C_{\mathcal{P}}(t)=\mathbb{E}\left[\braket{\hat{I}}_{x(t)}\braket{\hat{I}}_{x(0)}\right]-\mathbb{E}\left[\braket{\hat{I}}_{x(0)}\right]^{2}.$
(25)

In these expressions, all averages $\mathbb{E}[\bullet]$ are taken over
oscillator trajectories $x(t)$ sampled from the Langevin equation (14), with
initial condition $x(0)$ distributed according to the steady-state position
distribution, $P(x)$. Equation (24) derives from the $\mathcal{Q}$ component
of the correlation function, which is effectively delta-correlated in
comparison to the slow timescales of the oscillator dynamics. This term
represents the quantum shot noise of the current, $\Delta_{x}$, averaged over
the steady-state position distribution. Equation (25) is the $\mathcal{P}$
component of $C(t)$: it describes slow fluctuations of the current induced by
stochastic variations of the oscillator position, $x(t)$.

For our model, both the conditional average current $\braket{I}_{x}$ and
current noise $\Delta_{x}$ can be found using Landauer-Büttiker theory [58],
as

$\displaystyle\langle\hat{I}\rangle_{x}=\frac{1}{\pi}\int\differential
E\,\tau_{x}(E)\left[f_{L}(E)-f_{R}(E)\right],$ (26)
$\displaystyle\Delta_{x}=\frac{2}{\pi}\int\differential
E\tau_{x}(E)\left[f_{L}(E)(1-f_{L}(E))+f_{R}(E)(1-f_{R}(E)\right]$
$\displaystyle\qquad+\frac{2}{\pi}\int\differential
E\tau_{x}(E)\left[1-\tau_{x}(E)\right]\left(f_{L}(E)-f_{R}(E)\right)^{2}.$
(27)

Here, $f_{\alpha}(E)=\left[\mathrm{e}^{\beta(E-\mu_{\alpha})}+1\right]^{-1}$
is the Fermi-Dirac distribution for lead $\alpha$, while $\tau_{x}(E)$ is the
transmission function of the resonant-level model with the position-dependent
energy $\epsilon-Fx$ [58]:

$\tau_{x}(E)=\frac{\kappa_{L}(E)\kappa_{R}(E)}{\left|E-\epsilon+Fx-\tilde{\chi}_{L}(E)-\tilde{\chi}_{R}(E)\right|^{2},}$
(28)

where $\kappa_{\alpha}(E)$ is the spectral density [Eq. (5)] of lead
$\alpha=L,R$ and the corrresponding self-energy is

$\tilde{\chi}_{\alpha}(E)=\fint\frac{\text{d}E^{\prime}}{2\pi}\frac{\kappa_{\alpha}(E^{\prime})}{E-E^{\prime}}-\mathrm{i}\frac{\kappa_{\alpha}(E)}{2},$
(29)

with $\fint$ a principal-value integral.

##  III Results

###  III.1 Analysis of current fluctuations

####  III.1.1 Current trajectories and timescale separation

As described in Sec. II.3, the mechanical motion induces changes in the
measured current. Therefore, mechanical self-oscillations will lead to a
periodic current signal, which can be used to measure the passage of time. The
accuracy of such timekeeping is then limited by stochastic fluctuations in the
current signal, as characterised by the autocorrelation function $C(t)$. The
most relevant fluctuations for timekeeping are encoded in the slow part of the
correlation function $C_{\mathcal{P}}(t)$ [Eq. (25)] because this varies on
the same timescale as the clockwork dynamics. The fast fluctuations described
by $C_{\mathcal{Q}}(t)$ merely add a constant to the current noise spectrum,
as discussed at the end of this section.

Figure 2: Portion of a single current trajectory above the self-oscillation
voltage threshold, with $V=100\omega_{0}$. The inset shows the oscillations on
a shorter time scale where the oscillations appear more sinusoidal. The red
vertical lines indicate times when a tick is recorded. Parameters:
$\lambda=0.5\omega_{0}$, $\Gamma=10\omega_{0}$, $\beta\omega_{0}=0.1$,
$\omega_{L}=-\omega_{R}=2.5\omega_{0}$, and $\delta=5\omega_{0}$.

To characterise the properties of the measured current signal, we solve the
Langevin equation (14) using a stochastic Euler method and convert
trajectories $x(t)$ into a corresponding current signal
$\braket{\hat{I}}_{x(t)}$ via Eq. (26). Fig. 2 shows a sample current
trajectory within the parameter regime where self-oscillations occur. We use
the same parameters as in Ref. [49], where the transition to self-sustained
oscillations was found to occur at $V\approx 30\omega_{0}$. We observe near-
periodic oscillations of the current on shorter timescales (Fig. 2 inset),
with a peak-to-peak amplitude that varies significantly over longer timescales
(Fig. 2 main panel). We also see that the oscillation frequency is
approximately twice the natural frequency of the mechanical oscillator, and
that the maximum oscillation amplitude has a fixed value that does not
fluctuate. Both of these observations can be understood by the fact that the
maximum current is obtained whenever the oscillator’s position crosses $x=0$.
A natural way to measure time is therefore to register a tick of the clock
whenever the current reaches its maximum value, as illustrated in the inset of
Fig. 2. Note that while the precise position $x$ at which the current is
maximised depends on the specific choice of parameters, other parameters would
yield qualitatively similar results.

Figure 3: Auto-correlation function of the current for the same parameters as
Fig. 2 and three different voltages: $V=5\omega_{0}$ (top), $V=50\omega_{0}$
(middle), $V=100\omega_{0}$ (bottom). Two timescales emerge under high voltage
bias: a fast-decaying component due to phase fluctuations and and a longer-
lived component associated with amplitude noise.

Figure 3 shows the correlation function $C(t)$ for three different applied
voltages, spanning the transition from thermal noise at low bias to self-
sustained oscillations at high bias. At low voltage we observe decaying
oscillations, as would be expected for the thermal fluctuations of an
underdamped harmonic oscillator [59]. Meanwhile, at voltages above the
threshold for self-sustained oscillations, two distinct timescales can be
discerned. The first, fast timescale governs the decay of oscillations, which
results from phase diffusion on the limit cycle. The second, slow timescale
dictates the asymptotic decay of the correlation function: for example, in the
lower panel of Fig. 3, oscillations decay over a few hundred oscillation
periods, while the correlation function itself does not decay completely until
approximately $10^{4}$ periods have elapsed. This latter decay is monotonic
rather than oscillatory and can, therefore, be understood as a consequence of
long-lived amplitude correlations, as we now explain.

####  III.1.2 Modelling amplitude and phase fluctuations

To elucidate the different roles of amplitude and phase fluctuations, in
Appendix B we derive an effective equation of motion describing small
fluctuations of the oscillator trajectory away from the limit cycle. We
express the results in polar coordinates, i.e. $x=A\cos(\phi)$ and
$\dot{x}/\omega_{0}=-A\sin(\phi)$. An ideal limit cycle would have a constant
amplitude $A=A_{0}$ and a linearly increasing phase $\phi=\omega_{0}t$. Taking
into account fluctuations away from the limit cycle up to first order, we
obtain an Ornstein-Uhlenbeck process describing amplitude variations, and a
biased diffusion process describing the phase evolution, viz.

$\displaystyle\text{d}A=-\gamma_{A}(A-A_{0})\text{d}t+\sqrt{D_{A}}\text{d}w_{A},$
(30)
$\displaystyle\text{d}\phi=\omega_{0}\text{d}t+\sqrt{D_{\phi}}\text{d}w_{\phi}.$
(31)

Here, $\gamma_{A}$ is the damping rate governing the regression of amplitude
fluctuations back to the limit cycle $A_{0}$, $\differential w_{A}$ and
$\differential w_{\phi}$ are independent Wiener increments describing
amplitude and phase fluctuations, respectively, and $D_{A}$ and $D_{\phi}$ are
the corresponding diffusion coefficients. Our simplified model assumes that
the limit-cycle amplitude is large and that the timescale of phase and
amplitude diffusion is small compared to the oscillation frequency, which
justifies approximating the amplitude and phase noise as independent even
though they derive from the same Wiener process. The explicit expressions for
the above parameters are cumbersome (see Appendix B for details), but a
straightforward numerical evaluation shows that $\gamma_{A}\gg D_{\phi}$ in
the high-bias regime where Eqs. (30) and (31) are valid. This indicates that
small amplitude variations evolve significantly faster than phase decoherence,
and therefore this model alone cannot explain the timescale separation evident
in Fig. 3.

We thus turn to an alternative toy model that captures the large-scale
amplitude variations seen in Fig. 2. Specifically, the trajectory shown in
Fig. 2 exhibits rare periods where the oscillation amplitude reduces
dramatically, so that the current remains close to its maximum value for a
short time before returning to its previous oscillatory state. We model this
behaviour as a bistable process, where the stochastic current switches between
two states labelled by $i=1,2$, with corresponding currents $I_{1}(t)$ and
$I_{2}(t)$. This kind of bistable behaviour was predicted for
electromechanical self-oscillations in Ref. [44], and electromechanical
bistabilities have recently been observed experimentally [60]. We assume that
these two states have different steady-state averages
$\overline{I_{i}}=\mathbb{E}[I_{i}(t)]$, and that they are statistically
independent from each other and from the switching process itself. Under these
assumptions, we show in Appendix C that the steady-state autocorrelation
function takes the form

$\displaystyle
C_{\mathcal{P}}(t)=\sum_{i=1}^{2}p_{i}p_{i|i}(t)C_{i}(t)+C_{1\leftrightarrow
2}(t),$ (32) $\displaystyle C_{1\leftrightarrow
2}(t)=\sum_{i,j=1}^{2}p_{i}\left[p_{j|i}(t)-p_{j}\right]\overline{I_{i}}\,\overline{I_{j}},$
(33)

where $p_{i}$ is the steady-state probability to be in state $i$, $p_{j|i}(t)$
is the probability to transition $i\to j$ over a time $t$, and
$C_{i}(t)=\mathbb{E}[I_{i}(t)I_{i}(0)]-\overline{I_{i}}^{2}$ is the auto-
correlation function of $I_{i}(t)$.

The first term in Eq. (32) clearly quantifies correlations within each state.
The second term, $C_{1\leftrightarrow 2}(t)$, is the contribution from the
bistability itself: its time dependence arises only from the statistics of the
switching process, while it depends on the current only through the
conditional averages, $\overline{I_{i}}$. Therefore, the autocorrelation
function will exhibit a separation of timescales whenever the switching
process is much slower than phase and amplitude fluctuations of the current
within each state. In Appendix C we demonstrate this explicitly, showing that
this simple bistable model recovers the features of the auto-correlation
function above threshold in Fig. 3. Crucially, using the conservation of
probability, $\sum_{i}p_{i}=\sum_{j}p_{j|i}(t)=1$, one can show that
$C_{1\leftrightarrow 2}(t)\propto\overline{I_{1}}-\overline{I_{2}}$.
Therefore, it is important that the cycle-averaged current depends on the
oscillation amplitude, which is indeed the case in our setup (cf. Fig. 2).

####  III.1.3 Current noise power spectrum

Finally, it is instructive to analyse the current noise in the Fourier domain
via the power spectrum

$\mathfrak{S}(\omega)=\int\differential t\,\mathrm{e}^{i\omega t}C(t).$ (34)

In addition to a flat white-noise contribution from the quantum shot-noise
term $C_{\mathcal{Q}}(t)$, we observe two distinct peaks in the power spectrum
at $\omega=0$ and $\omega\approx 2\omega_{0}$, reflecting the separation of
timescales inherent to $C_{\mathcal{P}}(t)$. The height of both peaks grows
with voltage bias, but for different reasons. The zero-frequency peak is
associated with amplitude noise and its height scales with the lifetime of
amplitude correlations (see Appendix C), which becomes very large well above
threshold. This growth of the zero-frequency noise is well-known in other non-
equilibrium systems with long-lived metastable states [44, 61]. Meanwhile, the
peak in $\mathfrak{S}(\omega)$ at $\omega\approx 2\omega_{0}$ is related to
electromechanical oscillations, and is plotted in Fig. 4. The height of this
peak scales with the oscillation amplitude, which is of course larger in the
self-oscillating regime.

Figure 4: Power spectrum of the output current near the oscillation frequency
$\omega\approx 2\omega_{0}$, for four different voltage biases; all other
parameters are the same as Fig. 2. The frequency-independent contribution from
$C_{\mathcal{Q}}(t)$ is small and only makes a significant contribution for
smaller biases as shown in the inset.

As the voltage increases, the position of the finite-frequency peak shifts
slightly downwards in frequency, while its width becomes narrower. This
implies that the ticks of the clock become more regular in time, while their
frequency slightly decreases. As discussed in Sec. III.2, these results
indicate an increase in the clock’s accuracy as the rate of dissipation
increases, accompanied by a small decrease in the clock’s resolution. This is
consistent with previous studies of autonomous clocks [11, 14, 50, 16].

###  III.2 Accuracy and resolution from the tick distribution

We now turn to a more detailed analysis of the individual ticks of the clock.
The tick distribution or waiting-time distribution (WTD) is defined as the
probability, $W(\tau)$, of a time delay $\tau$ between two consecutive ticks.
In simple models of autonomous clocks, each tick is assumed to be generated by
an identical process, leading to the waiting times being independent and
identically distributed (i.i.d.) random variables. Then, the WTD fully
characterises the statistical properties of the time signal.

In particular, the mean and variance of the WTD are

$\displaystyle\mu=\int_{0}^{\infty}\differential\tau\,\tau W(\tau),$ (35)
$\displaystyle\sigma^{2}=\int_{0}^{\infty}\differential\tau\,(\tau-\mu)^{2}W(\tau).$
(36)

The clock’s resolution $\nu$ and accuracy $N$ are then defined as [11]

$\displaystyle\nu=\mu^{-1},$ (37) $\displaystyle
N=\frac{\mu^{2}}{\sigma^{2}}.$ (38)

The resolution is simply the rate at which the ticks are generated, while the
accuracy is the expected number of ticks before the clock’s reading is off by
one tick. To see this, suppose that $n$ ticks spaced by waiting times
$\tau_{1},\tau_{2},\ldots\tau_{n}$ have been generated. The clock’s time
reading is $t_{n}=\sum_{i=1}^{n}\tau_{i}$, with average
$\mathbb{E}[t_{n}]=n\mu$ and root-mean-square error $\Delta
t_{n}=\sqrt{n}\sigma$ due to the assumption of i.i.d. ticks. The accuracy is
then the effective number of ticks $N$ such that $\Delta t_{N}=\mu$, which
leads directly to Eq. (38).

Figure 5: Waiting time distributions of the ticks for different voltage
biases, with other parameters the same as Fig. 2. Dots are the data from the
numeric calculations, the solid curves are the lines of best fit to the
inverse Gaussian distribution (39), with fitting parameters $\mu$ and
$\sigma$.

We construct the WTD for our model numerically as a normalised histogram of
the waiting times between all pairs of consecutive ticks in the time series.
The results are plotted in Fig. 5 for three different voltages. Above
threshold, we find that the WTD is well approximated by an inverse Gaussian
(Wald) distribution

$W(\tau)=\sqrt{\frac{\mu^{3}}{2\pi\sigma^{2}\tau^{3}}}\text{exp}\left(-\frac{\mu\left(\mu-\tau\right)^{2}}{2\sigma^{2}\tau}\right).$
(39)

As shown in Ref. [62], Eq. (39) can be derived by considering the distribution
of the first-passage time $\tau$ for the oscillation phase on the limit cycle
to complete one full $2\pi$ rotation. This result holds when the phase is
described by a continuous-time biased random walk, as in Eq. (31). The solid
curves in Fig. 5 show the best fits of the numerical WTDs to Eq. (39), with
$\mu$ and $\sigma$ as fitting parameters.

Figure 6: Accuracy (red line), and resolution (blue line) for increasing bias
voltage, plotted against the entropy produced per tick of the clock.
Parameters are the same as Fig. 2.

Notably, as the bias across the system is increased, the mean of the
distribution increases and the variance of the distribution is decreased. This
indicates that the clock’s accuracy increases as the bias across the system is
increased, while its resolution decreases slightly. Increasing the bias also
increases the average rate of entropy production, which is proportional to the
power dissipated in the circuit:

$\dot{\Sigma}=\beta\langle\hat{I}\rangle V.$ (40)

Figure 6 shows the accuracy and resolution of the electromechanical clock
against the entropy produced per tick, $\Delta S_{\rm tick}=\dot{\Sigma}/\nu$.
We observe an approximately linear relationship between accuracy and entropy
production during the onset of self-sustained oscillations, showing that
accurate timekeeping comes at the cost of increased dissipation. However, as
the voltage is increased further, the growth of accuracy slows and begins to
saturate: at this point, additional dissipation gives little further advantage
to timekeeping.

Meanwhile, the clock’s resolution decreases and eventually saturates with
increasing dissipation. However, the change in resolution is minimal for this
model, which is unsurprising given that the tick period is fundamentally fixed
by the physical oscillator frequency. The small decrease in resolution is
probably due to effective non-linearities associated with the position
dependence mean-field force $\propto F\braket{\hat{n}}_{x}$ and damping rate
$\gamma_{x}$, whose effect will change as the oscillator explores a larger
region of phase space.

###  III.3  Allan variance

The measures of accuracy and resolution considered in Sec. III.2 are only
strictly applicable when the waiting times are i.i.d. random variables. The
latter condition is equivalent to the statement that the total number of ticks
is a renewal process [63]. However, this assumption does not exactly hold
because the limit cycle is mesoscopic in scale and subject to strong amplitude
fluctuations. While the amplitude and phase noise can be approximated as
independent for very small deviations from the limit cycle [see Sec. III.1.2],
this conclusion does not hold in general. Larger amplitude fluctuations will
become correlated with phase noise because they both derive from the same
underlying stochastic force induced by electron tunnelling. Therefore,
differences in the state of the system after each tick may cause the next tick
to be correlated with the previous ones.

Given that the ticks of our mesoscopic electromechanical clock are correlated,
we need more a sophisticated measure of performance than the accuracy and
resolution defined in Eqs. (37) and (38). We focus instead on the Allan
variance [53], which measures the frequency stability of a stochastic process
and is widely used to test the stability of sensors and clocks. The Allan
variance quantifies the deviation of a measured clock signal $\theta$, from a
second reference clock. For simplicity, we take this reference to be an ideal
clock whose reading is the coordinate time, $t$, so that the measured
deviation is

$X(t)=\theta(t)-t.$ (41)

Then, $y(t)=\differential X/\differential t$ is the instantaneous frequency
difference between the measured clock and the reference clock at a time $t$.
This frequency deviation is estimated by observing the clock’s reading over a
finite integration time $T$, leading to the average fractional frequency
difference

$\mathrm{bar}{y}(t,T)=\frac{1}{T}\int^{T}_{0}\text{d}t^{\prime}\,y(t+t^{\prime})=\frac{X(t+T)-X(t)}{T}.$
(42)

The Allan variance is then proportional to the mean-square deviation of
successive average fractional frequency differences

$\sigma^{2}_{y}(T)=\frac{1}{2}\mathbb{E}\left\\{\left[\mathrm{bar}{y}(t+T,T)-\mathrm{bar}{y}(t,T)\right]^{2}\right\\},$
(43)

thus quantifying how instability in the average frequency scales with
integration time. In the case of deterministic signals, where the frequency
remains constant, the anticipated frequency deviation is zero. Consequently, a
clock that closely adheres to deterministic behaviour — and is, therefore,
more accurate — would exhibit a smaller Allan variance.

We compute the Allan variance for our electromechanical clock by partitioning
the stochastic current trajectories into non-overlapping regions of equal
duration $T$. The measured deviation pertaining to the $n^{\rm th}$ region,
$X_{n}=X(n\tau)$, is evaluated by counting the number of registered ticks up
to time $t=n\tau$. Finally, the Allan variance is computed as

$\sigma_{y}^{2}(T)=\frac{1}{2T^{2}}\mathbb{E}\left[\left(X_{n-1}-2X_{n}+X_{n+1}\right)^{2}\right].$
(44)

The behaviour of the Allan variance as a function of $T$ may be sensitive to
the precise choice of $T$ values — we have repeated our calculations for
different values of $T$ to ensure that the results are robust.

Figure 7: Allan variance as a function of integration time $T$ for various
different voltage biases. Parameters are the same as Fig. 2.

The results for the Allan variance are plotted in Fig. 7, for a range of
different voltages. In general, the Allan variance decreases systematically
with $T$: evidently, integrating the clock signal over a longer time yields
greater precision. This dynamic highlights the trade-off between the clock’s
timekeeping accuracy and resolution [16]. Note that here $\sigma_{y}^{2}$
decreases indefinitely with $T$, whereas real experiments usually find that
the Allan variance begins to grow beyond a certain integration time. This is
typically due to the influence of low-frequency fluctuations, e.g. due to
temperature variations or the gradual degradation of the apparatus over time,
which are absent in our theoretical model.

At short and intermediate timescales, we observe from Fig. 7 that higher
voltage leads to a smaller values of $\sigma_{y}^{2}$, demonstrating a
quantitative link between accuracy and the rate of energy dissipation. It
follows that the relationship between timekeeping precision and entropy
production extends over multiple ticks of the clock, beyond the simple
accuracy measure of Eq. (38) that only accounts for the single-tick WTD. At
long times, all curves in Fig. 7 tend to the approximate scaling
$\sigma_{y}^{2}\sim T^{-1}$. This is consistent with known results for the
long-time Allan variance for a renewal process (see Appendix D or Ref. [30]):

$\sigma_{y}^{2}\to\frac{\mu}{NT},$ (45)

which is fully determined by the clock’s accuracy $N$ defined in Eq. (38). A
similar $1/T$ scaling is observed in Fig. 7 across essentially all times for
voltages below the threshold for self-oscillations.

Interestingly, however, the Allan variance initially decreases slightly faster
than $T^{-1}$ for voltages above threshold. This trend continues up to an
intermediate timescale on the order of a few thousand ticks, after which the
slope of $\sigma_{y}^{2}(T)$ changes dramatically and the advantage in
timekeeping precision due to self-oscillations is eventually lost altogether.
The non-trivial behaviour of the Allan variance as a function of $T$ indicates
that temporal correlations play a significant role in determining the clock’s
perfomance over multiple ticks.

###  III.4 Quantifying correlations between ticks

In this section we turn to explicitly quantifying correlations between the
ticks of the clock, which are responsible for the faster decrease of
$\sigma_{y}^{2}(T)$ above threshold. We consider two different measures of
statistical correlations. First, we consider the clock reading after $n$
ticks, $t_{n}=\sum_{i=1}^{n}\tau_{i}$, which is the sum of $n$ consecutive
waiting times $\tau_{1},\tau_{2},\ldots,\tau_{n}$. Let $P_{n}(t_{n})$ be the
corresponding probability distribution. If the ticks were generated by a
renewal process, so that the individual waiting times were i.i.d., this
distribution would be given by the $n$-fold convolution of single-tick WTDs,
i.e.

$\displaystyle\left[\mathop{\scalebox{1.5}{\raisebox{-0.86108pt}{$\circledast$}}}^{n}W\right](t_{n})=\int_{0}^{t_{n}}$
$\displaystyle\text{d}t_{1}\int_{t_{1}}^{t_{n}}\differential
t_{2}\cdots\int_{t_{n-2}}^{t_{n}}\text{d}t_{n-1}W(t_{1})\times$ $\displaystyle
W(t_{2}-t_{1})\times\cdots\times W(t_{n}-t_{n-1}).$ (46)

We use the relative entropy (Kullback-Leibler divergence) to measure the
distance between $P_{n}(t_{n})$ and the i.i.d. result in Eq. (III.4):

$D_{\rm
KL}(P_{n}||\mathop{\scalebox{1.5}{\raisebox{-0.86108pt}{$\circledast$}}}^{n}W)=\int_{0}^{\infty}\text{d}t\,P_{n}(t)\ln\left(\frac{P_{n}(t)}{\left[\displaystyle\mathop{\scalebox{1.5}{\raisebox{-0.86108pt}{$\circledast$}}}^{n}W\right](t)}\right).$
(47)

Since we are considering a stationary process, if the waiting times are not
i.i.d. then they are not independent. Equation (47) can thus be understood as
a measure of the total statistical correlations among all $n$ ticks.

Figure 8: The main panel shows the mutual information [Eq. (48)] between a
tick and the $m^{\rm th}$ subsequent tick, i.e. the pairwise correlations
between the waiting times $\tau_{i}$ and $\tau_{i+m}$. The inset shows the
relative entropy [Eq. (47)] between the distribution of the $n$-tick clock
reading $P(t_{n})$ and the $n$-fold convolution of single-tick WTD, which
measures the total correlations among all $n$ ticks. Parameters are the same
as Fig. 2.

The inset of Fig. 8 demonstrates that the ticks are indeed correlated: the
relative entropy grows as a function of $n$. Moreover, this growth is faster
for higher voltages, indicating that correlations between ticks become more
significant in the self-oscillation regime. Since these correlations should
have a finite range in time, we expect the total correlations to eventually
saturate as $n$ increases. Unfortunately, we are unable to reliably extend our
calculation of $D_{\rm
KL}(P_{n}||\mathop{\scalebox{1.5}{\raisebox{-0.86108pt}{$\circledast$}}}^{n}W)$
to higher values of $n$: the results become very noisy due to the limited
quantity of available numerical data and the relative complexity of the
$n$-fold convolution.

To estimate the timescale over which correlations decay, we therefore focus on
a simpler measure of correlations: the mutual information between pairs of
ticks that are $m$ ticks apart. This is defined as

$\mathcal{I}(\tau_{i}:\tau_{i+m})=2H[W(\tau_{i})]-H[W(\tau_{i},\tau_{i+m})],$
(48)

where $W(\tau_{i},\tau_{i+m})$ is the joint distribution of the waiting times
$\tau_{i}$ and $\tau_{i+m}$, $W(\tau_{i})$ is the WTD (identical for
$\tau_{i}$ and $\tau_{i+m}$) and $H[W(\tau)]=-\int\differential\tau W(\tau)\ln
W(\tau)$ denotes the Shannon entropy. Unlike the relative entropy, which
measures the total correlations among all ticks, Eq. (48) quantifies only
correlations between specific pairs of ticks and can be computed reliably over
long timescales using our numerical data. The mutual information is plotted in
the main panel of Fig. 8. This plot demonstrates that indeed the correlations
are short-ranged in time: $\mathcal{I}(\tau_{i}:\tau_{i+m})$ is largest for
small $m$ and at lower voltages. However, especially at voltages above
threshold, these correlations persist over long timescales on the order of
several thousand ticks. This is the same timescale at which the non-monotonic
behaviour of the Allan variance in Fig. 7 appears.

We thus conclude that $\sigma_{y}^{2}$ decreases more rapidly above threshold
because of temporal correlations between the ticks. These correlations provide
an advantage in timekeeping precision, so long as the integration time $T$ is
not much greater than the correlation time. By contrast, when $T$ far exceeds
the range of temporal correlations, the integrated signals in two neighbouring
time periods become effectively independent and the behaviour of the Allan
variance reduces to that of a renewal process, i.e. $\sigma_{y}^{2}\sim
T^{-1}$.

These correlations are likely dictated by several factors acting in tandem.
The timescale over which they persist appears to be far longer than the
timescale of pure phase decoherence apparent in Fig. 3, indicating that
amplitude fluctuations also play a role. This can be understood from the fact
that amplitude and phase fluctuations are generally not independent, meaning
that a variation in amplitude after one cycle can influence the time taken for
the next cycle to be completed. The timescale over which tick correlations
persist will then be determined by a combination of phase decoherence,
amplitude fluctuations near the limit cycle, and large-scale amplitude
variations due to bistability in the current signal.

##  IV Conclusions and outlook

In this work, we have introduced a self-contained model of an autonomous
clock, where the clockwork is an electromechanical self-oscillator and the
ticks of the clock are read out from the electrical current passing through
the system. A peculiar and novel feature of our setup is that the clockwork is
driven by the same irreversible process that is used for tick readout — the
tunnelling of electrons through the device. We have confirmed that the clock’s
accuracy increases with the rate of entropy production, as the current drives
the oscillator into a limit cycle with long-lived phase coherence. We have
identified that both phase and amplitude fluctuations play an important role
in determining the statistical properties of the measured current.

A key finding of our work is that the ticks of the clock are generally
correlated with each other and that these correlations can improve the clock’s
accuracy, as quantified by the Allan variance. It remains an open problem to
understand under what conditions such correlations can improve clock accuracy
in general. We note that temporal correlations have recently been considered
in the context of (classical) parameter estimation, where they may either help
or hinder depending on their nature [64]. Yet far more work is needed to unify
the theory of quantum ticking clocks with the standard framework of quantum
parameter estimation, especially beyond the asymptotic regime of infinitely
many samples (or in this case, ticks).

Our model is strongly motivated by modern electromechanical platforms, meaning
that our predictions could be probed directly in state-of-the-art experiments.
In the present work, we have limited ourselves to the quasi-adiabatic regime,
which is often the most experimentally relevant one and is also the simplest
regime to treat theoretically. However, the quasi-adiabatic regime is very far
from thermodynamically optimal, since many electrons tunnel through the device
per tick of the clock. This regime is also sub-optimal from the perspective of
readout, since the mechanical oscillations are “oversampled” by the tunnelling
events [16]. It is thus important to explore the complementary regime where
each oscillation is sampled by only a handful of electrons, approaching the
lower sensitivity limit of electrical readout and entailing thus a far smaller
rate of dissipation. At the level of theory and modelling, going beyond the
quasi-adiabatic limit requires a non-perturbative methodology capable of
tackling situations where the electromechanical coupling and the tunnelling
rates are comparable. In future work, we plan to apply the mesoscopic-
leads/extended-reservoir approach [65, 66] to address this problem.

Finally, it would be interesting to explore the thermodynamics of timekeeping
for other kinds of self-oscillators in the quantum regime. Conditions for
self-oscillations to emerge from open quantum dynamics have only recently been
established [67, 68, 69, 70], within the broader context of research on broken
time-translation symmetry [71, 72, 73]. However, a general recipe to exploit
such symmetry breaking for accurate timekeeping remains lacking.

###### Acknowledgements.

We are grateful to N. Ares, J. Dexter, M. Huber and F. Meier for informative
discussions. This work was funded by the Irish Research Council, the European
Research Council Starting Grant ODYSSEY (Grant Agreement No. 758403), the
EPSRC-SFI joint project QuamNESS, and the European Commission via the Horizon
Europe project ASPECTS (Grant Agreement No. 101080167). J.G. is supported by a
SFI-Royal Society University Research Fellowship. M.J.K. acknowledges the
financial support from a Marie Skłodwoska-Curie Fellowship (Grant Agreement
No. 101065974). A.S. acknowledges financial support from PNRR MUR project
PE0000023-NQSTI and Quantera project SuperLink. M.T.M. is supported by a Royal
Society-Science Foundation Ireland University Research Fellowship
(URF\R1\221571). We acknowledge the Irish Centre for High End Computing for
the provision of computational facilities. Some calculations were performed on
the Kelvin cluster maintained by the Trinity Centre for High Performance
Computing. This cluster was funded through grants from the Higher Education
Authority, through its PRTLI program.

##  Appendix A Derivation of the electronic correlation function

In this appendix, we derive the expression (25) for the electronic correlation
function. The starting point is the auxiliary density operator defined in Eq.
(20), i.e.

$\hat{K}(t)=\text{e}^{-i\hat{H}t}\hat{I}\hat{\rho}\text{e}^{i\hat{H}t}.$ (49)

Taking the Wigner transform as in Eq. (8), we have the equation of motion

$\frac{\partial}{\partial
t}\hat{K}(x,p,t)=(\mathcal{L}_{x}+\mathcal{H}+\mathcal{F})\hat{K}(x,p,t),$
(50)

with the superoperators given by Eqs. (10)–(12). The projection operator
$\mathcal{P}$ defined in Eq. (13), and its orthogonal complement
$\mathcal{Q}=1-\mathcal{P}$, obey the following properties [49]

$\displaystyle\mathcal{P}^{2}$ $\displaystyle=\mathcal{P},$ (51a)
$\displaystyle\mathcal{P}\mathcal{L}_{x}$
$\displaystyle=0=\mathcal{L}_{x}\mathcal{P},$ (51b)
$\displaystyle\mathcal{P}\mathcal{F}\mathcal{P}$ $\displaystyle=0,$ (51c)
$\displaystyle[\mathcal{P},\mathcal{H}]$
$\displaystyle=\frac{p}{m}\frac{\partial\hat{\rho}_{x}}{\partial
x}\tr[\bullet],$ (51d) $\displaystyle\mathcal{P}[\mathcal{P},\mathcal{H}]$
$\displaystyle=0.$ (51e)

As in Ref. [49], we assume that $\mathcal{L}_{x}\gg\mathcal{H},\mathcal{F}$,
which follows from the quasi-adiabatic conditions (7). Then the electronic
degrees of freedom can be assumed to relax quickly to the stationary state of
$\mathcal{L}_{x}$, while $\mathcal{H}$ and $\mathcal{F}$ can be treated
perturbatively.

To proceed, we define the compact notation
$\hat{r}(t)=\mathcal{P}\hat{K}(x,p,t)$ and
$\hat{q}(t)=\mathcal{Q}\hat{K}(x,p,t)$, suppressing phase-space arguments for
concision. Inserting appropriate factors of the identity
$1=\mathcal{Q}+\mathcal{P}$ on both sides of Eq. (50), and using the
properties (51), one has

$\displaystyle\frac{\partial\hat{r}}{\partial t}$
$\displaystyle=\mathcal{P}\mathcal{H}\hat{r}+\mathcal{P}\left(\mathcal{H}+\mathcal{F}\right)\hat{q},$
(52) $\displaystyle\frac{\partial\hat{q}}{\partial t}$
$\displaystyle=\mathcal{L}_{x}\hat{q}+\mathcal{Q}\left(\mathcal{H}+\mathcal{F}\right)\hat{q}+\mathcal{Q}(\mathcal{H}+\mathcal{F})\hat{r}.$
(53)

Formally integrating Eq. (53) yields

$\displaystyle\hat{q}(t)=\mathcal{G}(t)\hat{q}(0)+\int_{0}^{t}\differential
t^{\prime}\mathcal{G}(t-t^{\prime})\mathcal{Q}\left(\mathcal{H+\mathcal{F}}\right)\hat{r}(t^{\prime}),$
(54)
$\displaystyle\mathcal{G}(t)=\exp\left\\{\left[\mathcal{L}_{x}+\mathcal{Q}(\mathcal{H}+\mathcal{F})\right]t\right\\}.$
(55)

Substituting this back into Eq. (52), we obtain the closed equation of motion

$\displaystyle\frac{\partial\hat{r}}{\partial t}$
$\displaystyle=\mathcal{P}\mathcal{H}\hat{r}(t)+\int_{0}^{t}\differential
t^{\prime}\mathcal{P}\left(\mathcal{H}+\mathcal{F}\right)\mathcal{G}(t-t^{\prime})\mathcal{Q}\left(\mathcal{H}+\mathcal{F}\right)\hat{r}(t^{\prime})$
$\displaystyle\quad+\mathcal{P}\left(\mathcal{H}+\mathcal{F}\right)\mathcal{G}(t)\hat{q}(0).$
(56)

The first two terms on the right-hand side (RHS) of Eq. (A) are identical to
Eq. (A14) in the Supplemental Material of Ref. [49], and can be treated by
identical manipulations. In brief, we make the following approximations: (i)
since $\mathcal{H},\mathcal{F}\ll\mathcal{L}_{x}$ we write
$\mathcal{G}(t)\approx\mathrm{e}^{\mathcal{L}_{x}t}$ to lowest non-trivial
order in small quantities (Born approximation), (ii) assuming that the
integrand decays rapidly as a function of $t^{\prime}$, we approximate
$\hat{r}(t-t^{\prime})\approx\hat{r}(t)$ and extend the upper integration
limit to infinity (Markov approximation). After tracing over the electronic
degrees of freedom, therefore, the first two terms on the RHS of Eq. (A)
reduce to the Fokker-Planck generator

$\displaystyle\mathcal{L}_{\rm FP}K=$ $\displaystyle\frac{\partial}{\partial
p}\left(m\omega_{0}^{2}x-F\braket{\hat{n}}_{x}+\gamma_{x}p\right)K-\frac{p}{m}\frac{\partial
K}{\partial x}+\frac{D_{x}}{2}\frac{\partial^{2}K}{\partial p^{2}},$ (57)

where $K(x,p,t)=\tr\hat{K}(x,p,t)$. The same Fokker-Planck operator dictates
the evolution of the Wigner function, $\partial W/\partial t=\mathcal{L}_{\rm
FP}W$, which can therefore be reconstructed by sampling the equivalent
Langevin equation (14). However, unlike the Wigner function, the equation for
$K(x,p,t)$ has an additional inhomogenity deriving from the third term on the
RHS of Eq. (A), i.e.

$\frac{\partial K}{\partial t}=\mathcal{L}_{\rm
FP}K+\tr\left[\mathcal{P}\left(\mathcal{H}+\mathcal{F}\right)\mathcal{G}(t)\hat{q}(0)\right].$
(58)

This equation has the solution

$K(t)=\mathrm{e}^{\mathcal{L}_{\rm FP}t}K(0)+\int_{0}^{t}\differential
t^{\prime}\,\mathrm{e}^{\mathcal{L}_{\rm
FP}(t-t^{\prime})}\tr\left[\mathcal{P}\left(\mathcal{H}+\mathcal{F}\right)\mathcal{G}(t^{\prime})\hat{q}(0)\right].$
(59)

Since $\mathcal{P}\hat{K}(x,p,t)=K(x,p,t)\hat{\rho}_{x}$ and
$\mathcal{P}\hat{K}(x,p,t)$ is given by Eq. (54), we have a formal solution
for both components of $\hat{K}(x,p,t)$.

Now, since we are considering the steady-state current correlation function,
the initial condition of Eq. (50) is

$\hat{K}(x,p,0)=\hat{I}\hat{\rho}(x,p)\approx\hat{I}\hat{\rho}_{x}W(x,p),$
(60)

to lowest order in small quantities, where $\hat{\rho}(x,p)$ and $W(x,p)$
correspond to the stationary state. Applying the projection operators, we get
at lowest order

$\displaystyle\mathcal{P}\hat{K}(x,p,0)\approx\hat{\rho}_{x}\braket{\hat{I}}_{x}W(x,p),$
(61)
$\displaystyle\mathcal{Q}\hat{K}(x,p,0)\approx\left(\hat{I}-\braket{\hat{I}}_{x}\right)\hat{\rho}_{x}W(x,p).$
(62)

To the same order, the second term on the right-hand side of both Eqs. (54)
and (59) can be neglected and we can write
$\mathcal{G}(t)\approx\mathrm{e}^{\mathcal{L}_{x}t}$. Therefore, we finally
obtain the approximate solution

$\displaystyle\hat{K}(x,p,t)$
$\displaystyle=\mathcal{P}\hat{K}(x,p,t)+\mathcal{Q}\hat{K}(x,p,t)$
$\displaystyle\approx\hat{\rho}_{x}\mathrm{e}^{\mathcal{L}_{\rm
FP}t}\\!\left[\braket{W(x,p)}\right]+\mathrm{e}^{\mathcal{L}_{x}t}\\!\left[\left(\hat{I}-\braket{\hat{I}}_{x}\right)\hat{\rho}_{x}\right]\\!W(x,p).$
(63)

Plugging this into Eq. (21) yields

$\displaystyle\braket{\hat{I}(t)\hat{I}(0)}$ $\displaystyle=\int\differential
x\int\differential p\,\braket{\hat{I}}_{x}\mathrm{e}^{\mathcal{L}_{\rm
FP}t}\\!\left[\braket{\hat{I}}_{x}W(x,p)\right]$
$\displaystyle\quad+\int\differential x\int\differential
p\,W(x,p)\left[\braket{\mathrm{e}^{\mathcal{L}^{\dagger}_{x}t}(\hat{I})\hat{I}}_{x}-\braket{\hat{I}}^{2}_{x}\right].$
(64)

The first term in Eq. (A) can be interpreted as an autocorrelation function of
the stochastic variable $\braket{\hat{I}}_{x(t)}$ with respect to the Langevin
evolution of the oscillator. To see this, note that
$\mathrm{e}^{\mathcal{L}_{\rm FP}t}\\!\left[\braket{\hat{I}}_{x}W(x,p)\right]$
is the solution of the Fokker-Planck equation with initial condition
$\braket{\hat{I}}_{x}W(x,p)$. We define the Fokker-Planck propagator
$\Pi(x^{\prime},p^{\prime},t^{\prime}|x,p,t)$ as the conditional probability
to find the oscillator at the phase space point $(x^{\prime},p^{\prime})$ at
time $t^{\prime}$ given that it was at $(x,p)$ at time $t$. We may therefore
write the first term in Eq. (A) as

$\displaystyle\int\differential x\int\differential
p\,\braket{\hat{I}}_{x}\mathrm{e}^{\mathcal{L}_{\rm
FP}t}\\!\left[\braket{\hat{I}}_{x}W(x,p)\right]$
$\displaystyle=\int\differential x^{\prime}\int\differential
p^{\prime}\int\differential x\int\differential
p\,\braket{\hat{I}}_{x}\Pi(x,p,t|x^{\prime},p^{\prime},0)\braket{\hat{I}}_{x^{\prime}}W(x^{\prime},p^{\prime})$
$\displaystyle=\int\differential x^{\prime}\int\differential
p^{\prime}\int\differential x\int\differential
p\,\braket{\hat{I}}_{x}\braket{\hat{I}}_{x^{\prime}}P(x,p,t;x^{\prime},p^{\prime},0)$
$\displaystyle\equiv\mathbb{E}\left[\braket{\hat{I}}_{x(t)}\braket{\hat{I}}_{x(0)}\right],$
(65)

where
$P(x,p,t;x^{\prime},p^{\prime},0)=\Pi(x,p,t|x^{\prime},p^{\prime},0)W(x^{\prime},p^{\prime})$
is the joint probability to find the system initially (time $t=0$) at
$(x^{\prime},p^{\prime})$ and finally (time $t$) at $(x,p)$, while
$\mathbb{E}[\bullet]$ denotes an average over the joint distribution which can
be computed by sampling Langevin trajectories.

The second term in Eq. (A) represents the instantaneous current fluctuations
conditioned on the oscillator position $x$. Within the quasi-adiabatic limit,
this correlation decays very quickly compared to the slow timescales of the
oscillator. Thus, on timescales relevant for the oscillator, the real part of
this correlation function can by approximated by a Dirac delta function:

$\real\left(\braket{\mathrm{e}^{\mathcal{L}^{\dagger}_{x}t}(\hat{I})\hat{I}}_{x}-\braket{\hat{I}}^{2}_{x}\right)\approx\Delta_{x}\delta(t),$
(66)

where $\Delta_{x}$ is the steady-state current noise conditioned on $x$, i.e.

$\Delta_{x}=\int_{-\infty}^{\infty}\differential
t\,\real\left(\braket{\mathrm{e}^{\mathcal{L}^{\dagger}_{x}t}(\hat{I})\hat{I}}_{x}-\braket{\hat{I}}^{2}_{x}\right).$
(67)

Putting everything together, we obtain Eq. (25). Finally, we note that the
above derivation holds for an arbitrary operator $\hat{I}$ whose support is
restricted to the electronic system in the Schrödinger picture (otherwise Eq.
(60) does not hold).

##  Appendix B  Modelling small amplitude and phase fluctuations near the
limit cycle

In this appendix, we derive a simple description of phase and amplitude
fluctuations on the oscillator’s limit cycle. Our starting point is the
Langevin equation (14) governing the evolution of the oscillator. We first
transform the equation to polar coordinates. Using the substitution
$z\omega_{0}=v$ and converting the system to polar coordinates,
$x=A\text{cos}\phi$, $z=A\text{sin}\phi$ the equations of motion become

$\displaystyle\text{d}A$
$\displaystyle=\text{cos}\phi\text{d}x+\text{sin}\phi\text{d}z+\frac{\text{cos}^{2}\phi}{2A}\text{d}z^{2},$
(68) $\displaystyle\text{d}\phi=$
$\displaystyle-\frac{\text{sin}\phi}{A}\text{d}x+\frac{\text{cos}\phi}{A}\text{d}z-\frac{\text{cos}\phi\text{sin}\phi}{2A^{2}}\text{d}z^{2}.$
(69)

Here we have used Itô’s lemma, where
$\text{d}z^{2}=\frac{D_{x}}{\omega_{0}^{2}}\text{d}t$, $\text{d}x^{2}=0$ and
$\text{d}x\text{d}z=0$. Focusing first on the amplitude equation, the equation
of motion is given as

$\displaystyle\text{d}A=-m\gamma_{(A,\phi)}$ $\displaystyle
A\text{sin}^{2}(\phi)\text{d}t+\frac{F}{\omega_{0}}\langle
n\rangle_{(A,\phi)}\text{sin}(\phi)\text{d}t$
$\displaystyle+\frac{\sqrt{D_{(A,\phi)}}}{w_{0}}\text{sin}(\phi)\text{d}w+\frac{\text{cos}^{2}(\phi)}{2A}\frac{D_{(A,\phi)}}{\omega_{0}^{2}}\text{d}t.$
(70)

Here $dw$ is a Wiener increment, and the subscript $(A,\phi)$ is used to
highlight dependence on both amplitude and the phase. Assuming the amplitude
change over a single cycle is small, we can replace the deterministic terms
with their mean over a single oscillation and the stochastic term with its RMS
over a single oscillation;

$\displaystyle\text{d}A$
$\displaystyle=-\frac{A}{2\pi}\oint\Big{(}\gamma_{(A,\phi)}\text{sin}^{2}(\phi)\Big{)}\text{d}\phi\text{d}t+\frac{F}{2\pi\omega_{0}}\oint\Big{(}\langle
n\rangle_{(A,\phi)}\text{sin}(\phi)\Big{)}\text{d}\phi\text{d}t+\frac{1}{\sqrt{2\pi}\omega_{0}}\sqrt{\oint\Big{(}\sqrt{D_{(A,\phi)}}\text{sin}(\phi)\Big{)}^{2}\text{d}\phi}\text{d}w_{A}$
$\displaystyle\quad+\frac{1}{4\pi
A\omega_{0}^{2}}\oint\Big{(}D_{(A,\phi)}\text{cos}^{2}(\phi)\Big{)}\text{d}\phi\text{d}t.$
(71)

Here $\differential w_{A}$ is an effective coarse-grained Wiener increment.
Due to symmetry in the $\langle\hat{n}\rangle_{x}$ term, its average over a
single cycle will be zero. We can find the limit cycle of the system by taking
the expectation value of both sides of the equation and setting
$\frac{\Delta\mathbb{E}\left[A\right]}{\Delta t}=0$. The limit cycle is given
by the solution of the transcendental equation

$1=\frac{1}{2\mathbb{E}[A]^{2}\omega_{0}^{2}}\frac{\mathbb{E}\left[\oint\text{cos}^{2}(\phi)D_{(A,\phi)}\text{d}\phi\right]}{\mathbb{E}\left[\oint\text{sin}^{2}(\phi)\gamma_{(A,\phi)}\text{d}\phi\right]}.$
(72)

We can also truncate Eq. 71 to lowest order for the dissipation and the
fluctuation terms around the limit cycle. Using the substitutions

$\displaystyle\gamma_{A}=\frac{\differential}{\differential
A}\Bigg{(}-\frac{A}{2\pi}\oint\Big{(}\gamma_{(A,\phi)}\text{sin}^{2}(\phi)\Big{)}\text{d}\phi$
$\displaystyle\qquad\qquad\qquad+\frac{1}{4\pi
A\omega_{0}^{2}}\oint\Big{(}D_{(A,\phi)}\text{cos}^{2}(\phi)\Big{)}\text{d}\phi\Bigg{)}\Bigg{\rvert}_{A=A_{0}},$
(73)
$\displaystyle\sqrt{D_{A}}=\frac{1}{\sqrt{2\pi}\omega_{0}}\sqrt{\oint\left(\sqrt{D_{(A,\phi)}}\text{sin}(\phi)\right)^{2}\differential\phi}\Bigg{\rvert}_{A=A_{0}},$
(74)

Substituting these into the Eq. 71, one arrives at the much simpler equation

$\text{d}A=-\gamma_{A}\left(A-A_{0}\right)\text{d}t+\sqrt{D_{A}}\text{d}w_{A}.$
(75)

This equation describes the amplitude of the system as an Ornstein–Uhlenbeck
process [55].

We can perform a similar analysis on the equation of motion for the phase:

$\displaystyle\text{d}\phi=$
$\displaystyle-\omega_{0}\text{d}t-\gamma_{x}\text{cos}(\phi)\text{sin}(\phi)\text{d}t+\frac{\text{cos}\phi
F}{\omega_{0}A}\text{d}t$
$\displaystyle+\frac{\text{cos}(\phi)\sqrt{D_{x}}}{A\omega_{0}}\text{d}w-\frac{\text{cos}(\phi)\text{sin}(\phi)D_{x}}{2A^{2}\omega_{0}^{2}}\text{d}t.$
(76)

We will assume that the oscillator’s natural frequency is the dominate term in
the deterministic part of the equation. The phase evolution is given by

$\text{d}\phi=-\omega_{0}\text{d}t+\frac{\text{cos}(\phi)\sqrt{D_{x}}}{A\omega_{0}}\text{d}w.$
(77)

Assuming sufficiently large $\omega_{0}^{2}A/\sqrt{D_{x}}$ we can linearise
Eq. (77) by replacing the stochastic term prefactor with its RMS over a single
cycle so that

$\text{d}\phi=-\omega_{0}\text{d}t+\sqrt{D_{\phi}}\text{d}w_{\phi},$ (78)

here $\differential w_{\phi}$ is the effective Wiener increment originating
from the linearisation. The diffusion term is given by

$\sqrt{D_{\phi}}=\frac{1}{A\omega_{0}\sqrt{2\pi}}\sqrt{\oint\left[\text{cos}^{2}(\phi)D_{x}\right]\text{d}\phi}\Bigg{\rvert}_{A=A_{0}}.$
(79)

Note that in the main text near Eqs. (30) and (31) we have replaced
$\phi\to-\phi$ merely for clarity of presentation.

Figure 9: Ratio between the phase diffusion coefficient $D_{\phi}$ and the
amplitude damping rate $\gamma_{A}$ computed for small fluctuations away from
the ideal limit cycle. The relatively small values indicate that phase
coherence on the limit cycle is long-lived compared to the rate of amplitude
variation. Parameters are the same as Fig. 2.

Since the amplitude and phase noise are independent within our approximations,
it is straightforward to compute the steady-state position autocorrelation
function [55]

$\displaystyle\mathbb{E}[x(t)x(0)]-\mathbb{E}[x(0)]^{2}$
$\displaystyle=\mathbb{E}[A(t)A(0)]\mathbb{E}[\cos\phi(t)\cos\phi(0)]$
$\displaystyle=\left(A_{0}^{2}+\frac{D_{A}\mathrm{e}^{-2\gamma_{A}t}}{2\gamma_{A}}\right)\frac{\cos(\omega_{0}t)\mathrm{e}^{-D_{\phi}t/2}}{2}.$
(80)

This function exhibits oscillatory decay at a rate that is ultimately dictated
by the phase decoherence time, $D_{\phi}^{-1}$. This conclusion holds
irrespectively of the scale, $D_{A}$, and regression rate, $\gamma_{A}$, of
amplitude fluctuations. This point is discussed in more detail at the end of
Appendix A.

To compare the characteristic timescales of amplitude and phase fluctuations
that result from this description, we numerically compute the quantities
$\gamma_{A}$ and $D_{\phi}$ from Eqs. (B) and (79), respectively. Their ratio
is plotted in Fig. 9 as a function of the voltage bias, for values of $V$
above and beyond the threshold for self-oscillations. We see that for all
voltages considered, we have $D_{\phi}/4\gamma_{A}\ll 1$, implying that phase
coherence lasts longer than the characteristic timescale of amplitude
variations.

##  Appendix C  Bistable model of large amplitude fluctuations

In this appendix, we derive the autocorrelation function (32) for a bistable
current that switches between two states, finding that it reproduces the
features seen in Fig. 3 when the switching process is sufficiently slow.
Concretely, consider a stochastic process $I(t)$ that switches randomly
between two different oscillatory states labelled by $i=1,2$, giving rise to
stochastic sub-signals $I_{1}(t)$ and $I_{2}(t)$. We are interested in the
steady-state correlation function

$C(t)=\mathbb{E}[I(t)I(0)]-\mathbb{E}[I(0)]^{2}.$ (81)

The average signal is

$\mathbb{E}[I(0)]=p_{1}\overline{I_{1}}+p_{2}\overline{I_{2}},$ (82)

where $p_{i}$ is the steady-state probability to be in state $i$ and
$\overline{I_{i}}=\mathbb{E}[I_{i}(0)]$ is the steady-state average of the
corresponding sub-signal. Assuming that switching between the two states
occurs independently of the value of $I_{i}(t)$, we can write

$\mathbb{E}\left[I(t)I(0)\right]=\sum_{i,j=1}^{2}p_{j|i}(t)p_{i}\mathbb{E}\left[I_{j}(t)I_{i}(0)\right],$
(83)

where $p_{j|i}(t)$ is the conditional probability to be in state $j$ at time
$t$, given state $i$ at $t=0$. Here we take $t>0$ without loss of generality,
since the stationary correlation function obeys $C(t)=C(-t)$. We further
assume that the two sub-signals are statistically independent from each other,
i.e.

$\mathbb{E}[I_{1}(t)I_{2}(t^{\prime})]=\mathbb{E}[I_{1}(t)]\mathbb{E}[I_{2}(t^{\prime})].$
(84)

Now, splitting the sum in Eq. (83) into diagonal ($i=j$) and off-diagonal
($i\neq j$) terms, and subtracting the square of Eq. (82), we arrive at Eqs.
(32) and (33), i.e.

$\displaystyle C(t)=\sum_{i=1}^{2}p_{i}p_{i|i}(t)C_{i}(t)+C_{1\leftrightarrow
2}(t),$ (85) $\displaystyle C_{1\leftrightarrow
2}(t)=\sum_{i,j=1}^{2}p_{i}\left[p_{j|i}(t)-p_{j}\right]\overline{I_{i}}\,\overline{I_{j}},$
(86)

where $C_{i}(t)=\mathbb{E}[I_{i}(t)I_{i}(0)]-\overline{I_{i}}^{2}$ is the
auto-correlation function of $I_{i}(t)$.

The switching term $C_{1\leftrightarrow 2}(t)$ has the properties expected of
an autocorrelation function: in particular, it vanishes as $t\to\infty$ since
$p_{j|i}(t)\to p_{j}$ at long times. Using the conservation of probability,
$\sum_{i}p_{i}=\sum_{j}p_{j|i}(t)=1$, one can also show that

$C_{1\leftrightarrow
2}(t)=\left(\overline{I_{2}}-\overline{I_{1}}\right)\left(p_{1}\overline{I}_{1}\left[p_{2|1}(t)-p_{2}\right]-p_{2}\overline{I_{2}}\left[p_{1|2}(t)-p_{1}\right]\right).$
(87)

Therefore, the switching term is non-zero only if
$\overline{I_{2}}\neq\overline{I_{1}}$, i.e. the two sub-signals have a
different mean value. If this is the case, and if the switching probabilities
$p_{j|i}(t)$ evolve on a slow timescale compared to $C_{i}(t)$, then Eq. (85)
will exhibit two well-separated timescales as in Fig. 3.

To demonstrate this explicitly, we assume that the switching occurs via a
simple telegraph process described by the master equation

$\dot{p}_{1}(t)=\lambda_{2}p_{2}(t)-\lambda_{1}p_{1}(t)=-\dot{p}_{2}(t),$ (88)

where $\lambda_{i}$ are transition rates between the sub-signals. The steady-
state probabilities are given by

$\begin{pmatrix}p_{1}\\\
p_{2}\end{pmatrix}=\frac{1}{\lambda_{1}+\lambda_{2}}\begin{pmatrix}\lambda_{2}\\\
\lambda_{1}\end{pmatrix},$ (89)

The transition probability of switching from state $i$ to state $j\neq i$
after a time $t>0$ is given by the solution of Eq. (88) with $p_{i}(0)=1$,
i.e.

$p_{j|i}(t)=\frac{\lambda_{i}}{\lambda_{1}+\lambda_{2}}\left(1-\mathrm{e}^{-(\lambda_{1}+\lambda_{2})t}\right).$
(90)

Substituting these results into Eq. (87) gives the two-point correlation
function of the telegraph process

$C_{1\leftrightarrow
2}(t)=\frac{\left(\overline{I_{2}}-\overline{I_{1}}\right)^{2}\lambda_{1}\lambda_{2}\mathrm{e}^{-(\lambda_{1}+\lambda_{2})|t|}}{\left(\lambda_{1}+\lambda_{2}\right)^{2}}.$
(91)

Therefore, the total correlation function [Eq. (85)] is the sum of two terms:
an oscillatory function determined by the $C_{i}(t)$, which decays due to
amplitude and phase fluctuations within each state, and a pure exponential
decay given by Eq. (91). So long as the total switching rate
$\lambda_{1}+\lambda_{2}$ is slow compared to the decay of the $C_{i}(t)$, the
latter contribution will survive after the oscillations have died out. This is
precisely the qualitative behaviour above threshold seen in Fig. 3.

This model also predicts a large zero-frequency component to the current
noise, consistent with our observations in Sec. III.1.3. From Eq. (34) one has
[74]

$\mathfrak{S}(0)=2\int_{0}^{\infty}\differential t\,C(t).$ (92)

Therefore the switching term contributes a term that scales as

$\mathfrak{S}(0)\sim\frac{\lambda_{1}\lambda_{2}}{(\lambda_{1}+\lambda_{2})^{3}},$
(93)

which diverges as $\lambda_{i}\to 0$. For example, if
$\lambda_{1}=\lambda_{2}=\lambda$ then we simply have
$\mathfrak{S}(0)\sim\lambda^{-1}$, where $\lambda^{-1}$ is the lifetime of
each state.

Note that other models of slow amplitude fluctuations give rise to similar
behaviour. However, it is crucial that the cycle-averaged signal depends on
the amplitude, which is reflected here by the condition that
$\overline{I_{1}}\neq\overline{I_{2}}$. For example, the model of independent
amplitude and phase fluctuations developed in Appendix B yields a position
autocorrelation function given by Eq. (B). This does not show a separation of
timescales, even when the amplitude fluctuations are much slower than phase
fluctuations, i.e. $D_{A},\gamma_{A}\ll D_{\phi}$. This is because the average
of $x(t)=A(t)\cos\phi(t)$ around one cycle is zero, independently of the value
of $A(t)$.

To illustrate the importance of an amplitude-dependent cycle average, let us
instead consider a model of stochastic oscillations with an offset:

$I(t)=A(t)\left\\{\cos[\phi(t)]-b\right\\}+c,$ (94)

where $b$ and $c$ are constants, and $A(t)$ and $\phi(t)$ are described by an
Ornstein-Uhlenbeck process [Eq. (30)] and a biased diffusion process [Eq.
(31)], respectively. The resulting stochastic dynamics is qualitatively
similar to the behaviour seen in Fig. 2 in regions where the peak-to-peak
amplitude is large, with $c$ representing the maximum current and the scale
factor $b$ determining how much oscillations reduce the current on average. In
particular, the cycle average of this signal is now amplitude dependent and
given by $\mathbb{E}_{\phi}[I(t)]=c-bA(t)$. We find the autocorrelation
function for this model is given by

$C(t)=\left(A_{0}^{2}+\frac{D_{A}\mathrm{e}^{-2\gamma_{A}t}}{2\gamma_{A}}\right)\frac{\cos(\omega_{0}t)\mathrm{e}^{-D_{\phi}t/2}}{2}+\frac{b^{2}D_{A}\mathrm{e}^{-2\gamma_{A}t}}{2\gamma_{A}}.$
(95)

The first term above is identical to Eq. (B), but now another term appears
that depends only on the amplitude fluctuations. This latter term can be
understood as a continuous analogue of Eq. (86), where now $b$ instead of
$\overline{I_{1}}-\overline{I_{2}}$ determines how amplitude fluctuations
affect the cycle-averaged signal.

##  Appendix D Allan variance for a renewal clock

In this appendix we derive the asymptotic Allan variance for a clock whose
ticks are generated by a renewal process. Our derivation takes a different
approach from that of Ref. [30], although the results are equivalent up to a
rescaling. Let $n(t)$ denote the number of ticks generated after a coordinate
time $t$ has elapsed. Let us assume that $n(t)$ is a renewal process for which
the mean and variance are well known [63]

$\displaystyle\mathbb{E}\left[n(t)\right]=\frac{t}{\mu}+\mathcal{O}(t^{0}),$
(96)
$\displaystyle\mathbb{E}\left[n(t)^{2}\right]-\mathbb{E}[n(t)]^{2}=\frac{\sigma^{2}t}{\mu^{3}}+\mathcal{O}(t^{0}),$
(97)

where $\mu$ and $\sigma$ are the mean and variance of the waiting-time
distribution.

The clock reading at time $t$ is taken to be

$\theta(t)=\mu n(t),$ (98)

which is the “asymptotically unbiased” choice since then
$\mathbb{E}[\theta(t)]=t$ at long times. The fractional frequency deviation
averaged over a time $T$ is then, from Eq. (42),

$\mathrm{bar}{y}(t,T)=\frac{\mu}{T}\left[n(t+T)-n(t)\right].$ (99)

Now, the Allan variance is defined as

$\displaystyle\sigma^{2}_{y}(T)$
$\displaystyle=\frac{1}{2}\mathbb{E}\left\\{\left[\mathrm{bar}{y}(t+T,T)-\mathrm{bar}{y}(t,T)\right]^{2}\right\\}$
$\displaystyle=\frac{1}{2}\mathbb{E}\left[\mathrm{bar}{y}(t+T,T)^{2}\right]+\frac{1}{2}\mathbb{E}\left[\mathrm{bar}{y}(t,T)^{2}\right]$
$\displaystyle\quad-\mathbb{E}\left[\mathrm{bar}{y}(t+T,T)\mathrm{bar}{y}(t,T)\right]$
$\displaystyle\underset{T\to\infty}{\longrightarrow}\mathbb{E}\left[\mathrm{bar}{y}(t,T)^{2}\right]-\mathbb{E}\left[\mathrm{bar}{y}(t,T)\right]^{2}$
(100)

On the second line, the first two terms are equal because the process is
stationary. The third line then follows in the limit of large $T$ because
$\mathrm{bar}{y}(t+T,T)$ and $\mathrm{bar}{y}(t,T)$ become statistically
independent, i.e. the number of renewals in time interval $[t,t+T]$ becomes
independent of the number of renewals in the following time interval
$[t+T,t+2T]$. Therefore, the Allan variance for large $T$ is (up to scale
factors) nothing but the variance of, $n(t+T)-n(t)$, i.e. the number of
renewals within a time interval of duration $T$. Analogously to Eq. (97), this
result is known asymptotically to be ${\rm
Var}\left[n(t+T)-n(t)\right]=\sigma^{2}T/\mu^{3}+\mathcal{O}(T^{0})$ [63].
Finally, this yields the asymptotic Allan variance for a renewal process

$\sigma^{2}_{y}(T)\longrightarrow\frac{\sigma^{2}}{\mu T},$ (101)

which is equivalent to Eq. (45).

## References

  * Stray  _et al._ [2022] B. Stray, A. Lamb, A. Kaushik, J. Vovrosh, A. Rodgers, J. Winch, F. Hayati, D. Boddice, A. Stabrawa, A. Niggebaum, M. Langlois, Y.-H. Lien, S. Lellouch, S. Roshanmanesh, K. Ridley, G. de Villiers, G. Brown, T. Cross, G. Tuckwell, A. Faramarzi, N. Metje, K. Bongs, and M. Holynski, Nature 602, 590 (2022), number: 7898 Publisher: Nature Publishing Group.
  * Bothwell  _et al._ [2022] T. Bothwell, C. J. Kennedy, A. Aeppli, D. Kedar, J. M. Robinson, E. Oelker, A. Staron, and J. Ye, Nature 602, 420 (2022), number: 7897 Publisher: Nature Publishing Group.
  * Salecker and Wigner [1958] H. Salecker and E. P. Wigner, Physical Review 109, 571 (1958), publisher: American Physical Society.
  * Peres [1980] A. Peres, American Journal of Physics 48, 552 (1980).
  * Page and Wootters [1983] D. N. Page and W. K. Wootters, Physical Review D 27, 2885 (1983), publisher: American Physical Society.
  * Janzing and Beth [2003] D. Janzing and T. Beth, IEEE Transactions on Information Theory 49, 230 (2003), conference Name: IEEE Transactions on Information Theory.
  * Ranković  _et al._ [2015] S. Ranković, Y.-C. Liang, and R. Renner, Quantum clocks and their synchronisation - the Alternate Ticks Game (2015), arXiv:1506.01373 [quant-ph].
  * Woods  _et al._ [2019] M. P. Woods, R. Silva, and J. Oppenheim, Annales Henri Poincaré 20, 125 (2019).
  * Woods [2021] M. P. Woods, Quantum 5, 381 (2021), publisher: Verein zur Förderung des Open Access Publizierens in den Quantenwissenschaften.
  * Woods  _et al._ [2022] M. P. Woods, R. Silva, G. Pütz, S. Stupar, and R. Renner, PRX Quantum 3, 010319 (2022), publisher: American Physical Society.
  * Erker  _et al._ [2017] P. Erker, M. T. Mitchison, R. Silva, M. P. Woods, N. Brunner, and M. Huber, Phys. Rev. X 7, 031022 (2017).
  * Barato and Seifert [2016] A. C. Barato and U. Seifert, Physical Review X 6, 10.1103/physrevx.6.041053 (2016).
  * Milburn [2020] G. J. Milburn, Contemp. Phys. 61, 69 (2020), arXiv:2007.02217 .
  * Schwarzhans  _et al._ [2021] E. Schwarzhans, M. P. E. Lock, P. Erker, N. Friis, and M. Huber, Phys. Rev. X 11, 011046 (2021).
  * Dost and Woods [2023] A. P. T. Dost and M. P. Woods, Quantum advantages in timekeeping: dimensional advantage, entropic advantage and how to realise them via Berry phases and ultra-regular spontaneous emission (2023), arXiv:2303.10029 [quant-ph].
  * Meier  _et al._ [2023] F. Meier, E. Schwarzhans, P. Erker, and M. Huber, Fundamental accuracy-resolution trade-off for timekeeping devices (2023), arXiv:2301.05173 [quant-ph].
  * Barato and Seifert [2015] A. C. Barato and U. Seifert, Physical Review Letters 114, 158101 (2015), publisher: American Physical Society.
  * Gingrich  _et al._ [2016] T. R. Gingrich, J. M. Horowitz, N. Perunov, and J. L. England, Physical Review Letters 116, 120601 (2016), publisher: American Physical Society.
  * Garrahan [2017] J. P. Garrahan, Physical Review E 95, 032134 (2017), publisher: American Physical Society.
  * Terlizzi and Baiesi [2018] I. D. Terlizzi and M. Baiesi, Journal of Physics A: Mathematical and Theoretical 52, 02LT03 (2018), publisher: IOP Publishing.
  * Ptaszyński [2018] K. Ptaszyński, Physical Review B 98, 085425 (2018), publisher: American Physical Society.
  * Agarwalla and Segal [2018] B. K. Agarwalla and D. Segal, Physical Review B 98, 155438 (2018), publisher: American Physical Society.
  * Guarnieri  _et al._ [2019] G. Guarnieri, G. T. Landi, S. R. Clark, and J. Goold, Physical Review Research 1, 033021 (2019), publisher: American Physical Society.
  * Malabarba  _et al._ [2015] A. S. L. Malabarba, A. J. Short, and P. Kammerlander, New Journal of Physics 17, 045027 (2015), publisher: IOP Publishing.
  * Frenzel  _et al._ [2016] M. F. Frenzel, D. Jennings, and T. Rudolph, New Journal of Physics 18, 023037 (2016), publisher: IOP Publishing.
  * Woods and Horodecki [2023] M. P. Woods and M. Horodecki, Physical Review X 13, 011016 (2023), publisher: American Physical Society.
  * Ball  _et al._ [2016] H. Ball, W. D. Oliver, and M. J. Biercuk, npj Quantum Information 2, 10.1038/npjqi.2016.33 (2016), publisher: Nature Publishing Group.
  * Xuereb  _et al._ [2023] J. Xuereb, P. Erker, F. Meier, M. T. Mitchison, and M. Huber, The Impact of Imperfect Timekeeping on Quantum Control (2023), arXiv:2301.10767 [quant-ph].
  * Hauge and Støvneng [1989] E. H. Hauge and J. A. Støvneng, Reviews of Modern Physics 61, 917 (1989), publisher: American Physical Society.
  * Silva  _et al._ [2023] R. Silva, N. Nurgalieva, and H. Wilming, Ticking clocks in quantum theory (2023), arXiv:2306.01829 [quant-ph].
  * Pietzonka [2022] P. Pietzonka, Physical Review Letters 128, 130606 (2022), publisher: American Physical Society.
  * Jenkins [2013] A. Jenkins, Physics Reports Self-oscillation, 525, 167 (2013).
  * Bachtold  _et al._ [2022] A. Bachtold, J. Moser, and M. I. Dykman, Rev. Mod. Phys. 94, 045005 (2022).
  * Gorelik  _et al._ [1998] L. Y. Gorelik, A. Isacsson, M. V. Voinova, B. Kasemo, R. I. Shekhter, and M. Jonson, Physical Review Letters 80, 4526–4529 (1998).
  * Park  _et al._ [2000] H. Park, J. Park, A. K. L. Lim, E. H. Anderson, A. P. Alivisatos, and P. L. McEuen, Nature 407, 57 (2000).
  * Erbe  _et al._ [2001] A. Erbe, C. Weiss, W. Zwerger, and R. H. Blick, Physical Review Letters 87, 096106 (2001), publisher: American Physical Society.
  * Pistolesi and Fazio [2006] F. Pistolesi and R. Fazio, New Journal of Physics 8, 113 (2006).
  * Novotný  _et al._ [2003] T. Novotný, A. Donarini, and A.-P. Jauho, Physical Review Letters 90, 10.1103/physrevlett.90.256801 (2003).
  * Novotný  _et al._ [2004] T. Novotný, A. Donarini, C. Flindt, and A.-P. Jauho, Physical Review Letters 92, 10.1103/physrevlett.92.248302 (2004).
  * Blanter  _et al._ [2004] Y. M. Blanter, O. Usmani, and Y. V. Nazarov, Phys. Rev. Lett. 93, 136802 (2004).
  * Blanter  _et al._ [2005] Y. M. Blanter, O. Usmani, and Y. V. Nazarov, Phys. Rev. Lett. 94, 049904 (2005).
  * Clerk and Bennett [2005] A. A. Clerk and S. Bennett, New Journal of Physics 7, 238 (2005).
  * Bennett and Clerk [2006] S. D. Bennett and A. A. Clerk, Phys. Rev. B 74, 201301 (2006).
  * Usmani  _et al._ [2007] O. Usmani, Y. M. Blanter, and Y. V. Nazarov, Phys. Rev. B 75, 195312 (2007).
  * Wen  _et al._ [2020] Y. Wen, N. Ares, F. J. Schupp, T. Pei, G. A. D. Briggs, and E. A. Laird, Nature Physics 16, 75 (2020).
  * Urgell  _et al._ [2020] C. Urgell, W. Yang, S. L. De Bonis, C. Samanta, M. J. Esplandiu, Q. Dong, Y. Jin, and A. Bachtold, Nat. Phys. 16, 32 (2020), arXiv:1903.04892 .
  * Wächtler  _et al._ [2019] C. W. Wächtler, P. Strasberg, S. H. L. Klapp, G. Schaller, and C. Jarzynski, New Journal of Physics 21, 073009 (2019).
  * Strasberg  _et al._ [2021] P. Strasberg, C. W. Wächtler, and G. Schaller, Phys. Rev. Lett. 126, 180605 (2021), arXiv:2101.05027 .
  * Culhane  _et al._ [2022] O. Culhane, M. T. Mitchison, and J. Goold, Physical Review E 106, 10.1103/physreve.106.l032104 (2022).
  * Pearson  _et al._ [2021] A. N. Pearson, Y. Guryanova, P. Erker, E. A. Laird, G. A. D. Briggs, M. Huber, and N. Ares, Phys. Rev. X 11, 021029 (2021).
  * He  _et al._ [2023] X. He, P. Pakkiam, A. A. Gangat, M. J. Kewming, G. J. Milburn, and A. Fedorov, Quantum clock precision studied with a superconducting circuit (2023), arXiv:2207.11043 [quant-ph].
  * Manikandan [2023] S. K. Manikandan, Autonomous quantum clocks using athermal resources (2023), arXiv:2207.07909 [cond-mat, physics:gr-qc, physics:quant-ph].
  * Allan [1966] D. Allan, Proceedings of the IEEE 54, 221 (1966).
  * [54] C. W. Gardiner and P. Zoller,  _Quantum Noise_ , 3rd ed. (Springer, Berlin).
  * Gardiner [2009] C. Gardiner,  _Stochastic Methods_ (Springer Berlin Heidelberg, 2009).
  * Vigneau  _et al._ [2023] F. Vigneau, F. Fedele, A. Chatterjee, D. Reilly, F. Kuemmeth, M. F. Gonzalez-Zalba, E. Laird, and N. Ares, Applied Physics Reviews 10, 021305 (2023).
  * Clerk  _et al._ [2010] A. A. Clerk, M. H. Devoret, S. M. Girvin, F. Marquardt, and R. J. Schoelkopf, Reviews of Modern Physics 82, 1155 (2010), publisher: American Physical Society.
  * Ryndyk [2016] D. A. Ryndyk, Green functions, in  _Theory of Quantum transport at nanoscale an introduction_ (Springer International Publishing, 2016).
  * Risken [1996] H. Risken,  _The Fokker-Planck Equation: Methods of Solution and Applications_, edited by H. Haken, Springer Series in Synergetics, Vol. 18 (Springer, Berlin, Heidelberg, 1996).
  * Tabanera-Bravo  _et al._ [2022] J. Tabanera-Bravo, F. Vigneau, J. Monsel, K. Aggarwal, L. Bresque, F. Fedele, F. Cerisola, G. A. D. Briggs, J. Anders, A. Aufèves, J. M. R. Parrondo, and N. Ares, Stability of long-sustained oscillations induced by electron tunneling (2022), arXiv:2211.04074 [cond-mat, physics:quant-ph].
  * Kewming  _et al._ [2022] M. J. Kewming, M. T. Mitchison, and G. T. Landi, Physical Review A 106, 033707 (2022), publisher: American Physical Society.
  * Aminzare  _et al._ [2019] Z. Aminzare, P. Holmes, and V. Srivastava, in  _2019 IEEE 58th Conference on Decision and Control (CDC)_ (IEEE, 2019).
  * Cox [1967] D. R. Cox,  _Renewal theory_ (Methuen, 1967).
  * Radaelli  _et al._ [2023] M. Radaelli, G. T. Landi, K. Modi, and F. C. Binder, New Journal of Physics 25, 053037 (2023), arXiv:2206.00463 [cond-mat, physics:quant-ph, stat].
  * Brenes  _et al._ [2020] M. Brenes, J. J. Mendoza-Arenas, A. Purkayastha, M. T. Mitchison, S. R. Clark, and J. Goold, Phys. Rev. X 10, 031040 (2020).
  * Lacerda  _et al._ [2023] A. M. Lacerda, A. Purkayastha, M. Kewming, G. T. Landi, and J. Goold, Physical Review B 107, 195117 (2023), publisher: American Physical Society.
  * Chan  _et al._ [2015] C.-K. Chan, T. E. Lee, and S. Gopalakrishnan, Physical Review A 91, 051601 (2015), publisher: American Physical Society.
  * Iemini  _et al._ [2018] F. Iemini, A. Russomanno, J. Keeling, M. Schirò, M. Dalmonte, and R. Fazio, Physical Review Letters 121, 035301 (2018), publisher: American Physical Society.
  * Buča  _et al._ [2019] B. Buča, J. Tindall, and D. Jaksch, Nature Communications 10, 1730 (2019), number: 1 Publisher: Nature Publishing Group.
  * Guarnieri  _et al._ [2022] G. Guarnieri, M. T. Mitchison, A. Purkayastha, D. Jaksch, B. Buča, and J. Goold, Physical Review A 106, 022209 (2022), publisher: American Physical Society.
  * Wilczek [2012] F. Wilczek, Physical Review Letters 109, 160401 (2012), publisher: American Physical Society.
  * Sacha and Zakrzewski [2018] K. Sacha and J. Zakrzewski, Reports on Progress in Physics 81, 016401 (2018), arXiv:1704.03735 [cond-mat, physics:quant-ph].
  * Zaletel  _et al._ [2023] M. P. Zaletel, M. Lukin, C. Monroe, C. Nayak, F. Wilczek, and N. Y. Yao, Colloquium: Quantum and Classical Discrete Time Crystals (2023), arXiv:2305.08904 [cond-mat, physics:nlin, physics:quant-ph].
  * Landi  _et al._ [2023] G. T. Landi, M. J. Kewming, M. T. Mitchison, and P. P. Potts, Current fluctuations in open quantum systems: Bridging the gap between quantum continuous measurements and full counting statistics (2023), arXiv:2303.04270 [cond-mat, physics:quant-ph].