\section{INTRODUCTION}
A \emph{flywheel} is a device that stores kinetic energy in the rotational motion of the wheel and supplies it on demand.
In many devices the flywheel is an essential component for extracting work from an engine.
The main tasks of a flywheel are twofold: transducing discrete energy into continuous power and storing useful \emph{work}.
This energy reserve can be rapidly drained on demand, ultimately gaining more power than the power supply charging it can deliver.
Miniaturizing heat engines and refrigerators received much attention in the past decade.
Experimental set-ups of such devices were constructed in the micrometers domain \cite{clemens12,steeneken11}, and recently the operation of a single-atom heat engine was reported \cite{rossnagel15}.
Many theoretical studies of these devices were extended to the quantum domain, concentrated on the study of efficiency, power extraction and the study of thermodynamic laws (see reviews  \cite{levy14a,goold15,janet15,millen16,gelbwaser15b} and references therein).
Work extraction from quantum systems was also studied extensively \citep{alahverdyan04,alicki13,skrzypczyk14}.
\par
Any realistic engine is regulated by monitoring and a feedback loop.
The purpose is to control its timing, adjust its frequency and amplitude to match the other parts of the device, and  compensate for unpredictable disturbances.
Recent theoretical studies demonstrated that quantum properties such as coherence and correlations enhance the work  extracted from the system \cite{uzdin15,marti15,korzekwa15,scully03}.
Future quantum technologies aiming to exploit these quantum features will encounter the issue of regulating the device.
Standard ideal quantum measurements will demolish these features.
Therefore, to overcome this problem a conceivable approach to regulate the quantum device is by continuous weak measurements (\emph{monitoring}) and \emph{feedback control}.
Another fundamental problem which is demonstrated in this study and that is resolved by monitoring and feedback control is the unlimited entropy increase of the work repository, i.e. proliferating fluctuations catastrophically heat up the flywheel.
\par
In this paper we introduce the concept of a quantum flywheel as an integral part of a quantum \emph{heat engine}.
The flywheel is composed of a quantum harmonic oscillator (HO) interacting with a two qubit quantum heat engine.
If the flywheel was simply replaced by an external classical field then the engine would operate in steady state and power can continuously be extracted from the engine(see Appendix E for details).
However, if the task is to \emph{store} work in a quantum HO by coupling it to the engine,  the HO will be subject to a fatal growth of fluctuations and establishment of \emph{steady state} is impossible.
The HO is unstable even when an external \emph{driving} field is utilized to extract power and stabilize it.
By applying monitoring and feedback control we obtain a steady state operation for the flywheel, continuously gaining power, and storing useful work in the flywheel that later can be extracted.
\par
Monitored and controlled quantum heat engines are still to be realized experimentally, however, the individual components already exist.
Quantum monitoring and feedback control experiments exists for various HO's such as electromagnetic cavity, nanomechanical oscillators, trapped particles and superconducting circuits, see review \cite{nori14} and reference therein.
Single microscopic quantum heat engines realisations are still under development with only few examples as of today \cite{koski14,rossnagel15}.
\section{Heat Engine Operation}
 The basic concept of a quantum heat engine (similar to the classical one) is consists of two thermal heat baths at different temperature, a working medium, and a work repository.
In the quantum counterpart the working medium is quantized and the work repository can be either an external classical field \cite{levy14a} or it can be quantized as well \cite{gelbwaser15}.
Here we consider the operation of a continuous quantum engine for which the heat baths and the work repository are coupled simultaneously and continuously to the working medium \cite{levy14a}.
The working medium is comprised of two qubits, with the Hamiltonians $\hat{H}_a=\omega_h\hat{a}^{\dagger}\hat{a}$ and $\hat{H}_b=\omega_c\hat{b}^{\dagger}\hat{b}$ respectively.
Each qubit is weakly coupled to a different heat bath with the inverse temperature $\beta_h$ and $\beta_c$, where the indexes $h$ and $c$ stand for $hot$ and $cold$.
The dynamics of the qubits follow the standart thermalizing master equation of Lindblad-Gorini- Kossakowski-Sudarshan (LGKS) \cite{lindblad76,gorini276}.
The asymptotic two-qubit state $\hat\rho_h^\infty\otimes\hat\rho_c^\infty$ is the product of the thermal equilibrium Gibbs states of the two qubits, respectively at hot and cold temperatures $1/\beta_h$ and $1/\beta_c$.
Satisfying the heat engine conditions, $\beta_h/\beta_c < \omega_c/\omega_h<1$,  population inversion is obtained between the 3rd level $\ket{10}$ and the 2nd level $\ket{01}$ (see Fig.\ref{fig:schem}a).
The populations of these states are given by $p_{10}=n_h(1-n_c)$ and  $p_{01}=n_c(1-n_h)$.
Here, $n_{h(c)}=\left[\exp(\beta_{h(c)}\omega_{h(c)})+1\right]^{-1} $  are the thermal occupation numbers in $\hat\rho_{h(c)}^\infty$.
The 2nd and the 3rd levels are treated as an effective two level system (TLS) with the energy gap $\omega_o = \omega_h-\omega_c $ (we take $\hbar=k_B=1$).
The state of this TLS is a Gibbs state with a negative effective temperature
\begin{equation}
\frac{1}{\beta_e^-}=\frac{\omega_h-\omega_c}{\beta_h\omega_h-\beta_c\omega_c}<0.
\end{equation}
\par
We exploit the TLS population inversion to "charge" a quantum harmonic oscillator (HO) with useful work.
The Hamiltonian of the HO and the TLS-HO interaction Hamiltonian are given by $\hat H_o=\omega_o \hat{c}^{\dagger}\hat{c}$, and $\hat{K}=ig(\hat{a}^{\dagger}\hat{b}\hat{c}- \hat{a}\hat{b}^{\dagger}\hat{c}^{\dagger})$, respectively.
Given that the termalization time of the qubits is much shorter then the internal time scale, then
\begin{equation}
\label{eq:weak_coupling}
g\sqrt{\ave{\hat{c}^{\dagger}\hat{c}}+1}\ll\Gamma_{h(c)} \left(1+\exp(-\beta_{h(c)}\omega_{h(c)})\right),
\end{equation}
the TLS can be considered heuristically as a heat bath with negative temperature weakly coupled to the HO. We prove that indeed the state $\hat\rho$ of the HO  satisfies the standard thermalizing master equation extended to negative temperature $1/\beta_e^-$, which in the interaction picture of $\hat H_o$ takes the form
\begin{equation}
\label{eq:ME_e}
\frac{d\hat\rho}{dt}=\mathcal{L}_e\hat\rho\equiv
\Gamma_e(\hat{c}\hat\rho\hat{c}^\dagger-\mathcal{H}\hat{c}^\dagger\hat{c}\hat\rho)
      +\Gamma_e e^{-\beta_e^-\omega_o}(\hat{c}^\dagger\hat\rho\hat{c}-\mathcal{H}\hat{c}\cop^\dagger\hat\rho).
\end{equation}
The damping rate $\Gamma_e$ is proportional to the squared coupling $g^2$, and depends on the parameters of the engine, such as the occupations $n_{h(c)}$ and the rates $\Gamma_{h(c)}$ (see Appendix A).
The notation $\mathcal{H}$ stands for the Hermitian part of everything coming after it (different from the convention in  ref \cite{wisemanbook}).
A rigorous derivation of equation (\ref{eq:ME_e}) can be found in Appendix A.
\par
Since $\beta_e^- < 0$ the master equation \eqref{eq:ME_e} has no steady state solution, energy will constantly flow into the flywheel.
The parameters containing the superscript $^-$ are negative.
The standard equations remain valid for the
mean amplitude $\ave{\hat{c}}_t$ and the occupation $\ave{\hat{c}^{\dagger}\hat{c}}_t$:
\begin{eqnarray}
\label{eq:c_t}
\frac{d\ave{\hat{c}}_t}{dt}&=&-\left(\kappa_e^-+i\omega_o\right)\ave{\hat{c}}_t\\
\label{eq:n_t}
\frac{d\ave{\hat{c}^\dagger\hat{c}}_t}{dt}&=&-2\kappa_e^-\ave{\hat{c}^\dagger\hat{c}}_t+\Gamma_e e^{-\beta_e^-\omega_o}
\end{eqnarray}
where the amplitude damping rate
\begin{equation}
\label{eq:kappa_neg}
\kappa_e^- = \frac{1}{2}\Gamma_e\left(1-e^{-\beta_e^-\omega_o}\right)
\end{equation}
takes negative values since $\beta_e^-<0$ (see Appendix A).
Therefore both $\ave{\hat{c}}_t$ and $\ave{\hat{c}^\dagger\hat{c}}_t$ (and all higher moments) diverge exponentially with time (see Fig.\ref{fig:schem}b)
resulting in the instability of the dynamics against small perturbations.
In particular, an initial Gibbs state maintains its form but with an exponentially growing temperature
$1/\beta_t=\omega_o/\log\left(1+ \ave{\hat{c}^{\dagger}\hat{c}}_t^{-1} \right)$.
Thus, $\hat\rho(t)\propto\exp\left(-\beta_t\omega_o\hat{c}^\dagger\hat{c}\right)$
is an unstable solution of the master eq. \eqref{eq:ME_e}. Any small
perturbation will divert it from the class of Gibbs states. A more general
class of solutions, displaced Gibbs states
\begin{equation}
\hat\rho(t)\propto\exp\left[-\beta_t\omega_o
(\hat{c}-\ave{\hat{c}}_t)^\dagger(\hat{c}-\ave{\hat{c}}_t) \right],\nonumber
\end{equation}
with effective temperature
 $1/\beta_t=\omega_o/\log\left[1+ \left( \ave{\hat{c}^{\dagger}\hat{c}}_t-|\ave{\hat{c}}_t|^2\right)^{-1} \right]$ will, in principle, be suitable for work extraction.
But this option is misleading since the instability of the above solutions is not yet resolved.
\par
A reasonable approach to stabilize the flywheel while extracting additional power is achieved by driving the HO via a resonant oscillating external field.
The field is expressed by the time dependent Hamiltonian, $\hat H_d(t)=-i\epsilon_d\hat{c}^{\dagger}e^{-i\omega_o t}+\mathrm{H.c.}$.
The master equation in the interaction picture \eqref{eq:ME_e} becomes modified by a static Hamiltonian term (see Appendix B):
\begin{equation}
\label{eq:ME_e_d}
\frac{d\hat\rho}{dt}=\mathcal{L}_e\hat\rho-\epsilon_d\left[\hat{c}^{\dagger} -\hat{c},\hat\rho_o\right].
\end{equation}
Indeed, equation (\ref{eq:ME_e_d}) leads to a stationary amplitude with a rotating phase:
$\ave{\hat{c}}_t=-(\epsilon_d/\kappa_e^-)e^{-i\omega_o t}$.
Nevertheless, the stationary state remains unstable,
the occupation number and higher moments  diverge invariably.
Driving in itself can not solve the instability issue.
Unlimited growth of quantum and thermal fluctuations must be
suppressed by active control of the flywheel.

\begin{figure}
\centering
\begin{tabular}{cc}
\subfloat[]{\includegraphics[height=3cm]{Picture33.png}}&
\subfloat[]{\includegraphics[height=3cm]{Picture44.png}}
\\
\subfloat[]{\includegraphics[height=3cm]{Picture22.png}}&
\subfloat[]{\includegraphics[height=3cm]{Picture11.png}}
\end{tabular}
\caption{Scheme of the different components and operations of the heat engine with a flywheel.
(a) The state of two qubits of the heat engine, coupled to heat baths at temperatures $T_h$ and $T_c$, is represented as a two-qubit state with population inversion between the 2nd and 3rd energy levels. The size of the sphere represents the population in each level.
(b) The population inversion in the engine corresponds to a heat bath with the inverse negative temperature $\beta_e^-$. This bath is coupled to the harmonic oscillator (flywheel), increasing exponentially its energy and the width of phase-space probability distribution.
(c) Measurement of the quadratures of the harmonic oscillator, resulting in the signal $\bar{c}$. The signal is then fed back to the oscillator to ensure a steady-state.
(d) Energy flow chart of the different components in steady-state of the flywheel.}
\label{fig:schem}
\end{figure}

\section{Measurement and Feedback Control}
Continuous measurement, i.e., monitoring, is  the first task towards implementing feedback control \cite{wisemanbook}.
By applying monitoring and feedback control we can stabilize the  flywheel and charge it with useful work.
Consider a time-continuous measurement of  both quadratures $\hat{x}=\frac{1}{\sqrt{2}}(\hat{c}^\dagger +\hat{c})$ and $\hat{y}=\frac{i}{\sqrt{2}}(\hat{c}^\dagger -\hat{c})$ of the HO.
Generalizing the result of \cite{diosi88}, we simultaneously monitor $\hat{x}$ and $\hat{y}$ (see Fig.\ref{fig:schem}c).
The dynamics is described by a stochastic master equation (SME) for the density operator $\hat{\sigma}$ conditioned on both measurement signals $\bar{x},\bar{y}$ (see Appendix C).
The stochastic mean ($\mathbf{M}$) of the conditional state yields the unconditional state, i.e., $\mathbf{M}\hat{\sigma}=\hat\rho$ satisfying a corresponding master equation of the usual LGKS structure.
It differs from the master equation eq.(\ref{eq:ME_e_d})
by the additional monitoring term
\begin{equation}
\label{eq:L_m}
\mathcal{L}_m\hat\rho= \frac{\gamma_m}{4}(\hat{c}\hat\rho\hat{c}^\dagger-\mathcal{H}\hat{c}^\dagger\hat{c}\hat\rho + \hat{c}^\dagger\hat\rho\hat{c}-\mathcal{H}\hat{c}\cop^\dagger\hat\rho),
\end{equation}
 where $\gamma_m$ is the measurement strength. This generator corresponds to an infinite temperature heat bath.
Hence, the act of monitoring additionally heats the flywheel, and
contributes to the undesirable proliferating fluctuations of the HO.
\par
Stabilization is accomplished by a feedback loop conditioned on the measured signals $\bar{x},\bar{y}$.
As a result, the HO is kept in the vicinity of the constant rotating amplitude set by the external driving.
The feedback Hamiltonian in the  Schr\"{o}dinger picture is given by $\hat H_f(t)=-i\kappa_f\bar{c}(t)\hat{c}^\dagger +\mathrm{H.c.}$, where $\bar{c}=\frac{1}{\sqrt{2}}(\bar{x}+i\bar{y})$ is the complex representation of the two real signals $\bar{x}$ and $\bar{y}$, and $\kappa_f$ is the feedback strength.
By setting the value of $\kappa_f$ the steady-state of the flywheel is guaranteed.
The feedback is applied on top of the monitored evolution,
$\hat{\sigma}+d\hat{\sigma}\rightarrow e^{-i\hat H_f dt}(\hat{\sigma}+d\hat{\sigma}) e^{i\hat H_f dt}$,
yielding a SME for the conditional state (see Appendix D for details).
The master equation of the unconditional state reads
\begin{equation}
\label{eq:ME_e_d_m_f}
\frac{d\hat\rho}{dt}=(\mathcal{L}_e+\mathcal{L}_m+\mathcal{L}_f)\hat\rho-\epsilon_d\left[\hat{c}^{\dagger}-\hat{c},\hat\rho\right].
\end{equation}
The dissipative contribution of the feedback is:
\begin{eqnarray}
\label{eq:L_f}
\mathcal{L}_f\hat{\sigma}=
 \left(\frac{\kappa_f^2}{\gamma_m}+\kappa_f\right)\left(\hat{c}\hat{\sigma}\hat{c}^\dagger-\mathcal{H}\hat{c}^\dagger\hat{c}\hat{\sigma}\right)
 \\  \nonumber
 + \left(\frac{\kappa_f^2}{\gamma_m}-\kappa_f\right)\left(\hat{c}^\dagger\hat{\sigma}\hat{c}-\mathcal{H}\hat{c}\cop^\dagger\hat{\sigma}\right).
\end{eqnarray}

For $\kappa_f>\gamma_m$ this corresponds to a thermal bath of positive temperature.  Entering the regime $0<\kappa_f<\gamma_m$,
the cooling effect of $\mathcal{L}_f$ within the sum  $\mathcal{L}_e+\mathcal{L}_m+\mathcal{L}_f$ becomes enhanced although
$\mathcal{L}_f$ ceases to be a mathematically correct dissipator in itself.
Eq.(\ref{eq:ME_e_d_m_f}) can be written in a compact form,
\begin{eqnarray}
\label{eq:ME_final}
\begin{array}{lr}
\frac{d\hat\rho}{dt}=
 \Gamma\left(\hat{c}\hat\rho\hat{c}^\dagger-\mathcal{H}\hat{c}^\dagger\hat{c}\hat\rho\right)
+\Gamma e^{-\beta\omega_o}\left(\hat{c}^\dagger\hat\rho\hat{c}-\mathcal{H}\hat{c}\cop^\dagger\hat\rho\right)
\\
 ~~~ - ~\epsilon_d[\hat{c}^\dagger-\hat{c},\hat\rho]
\end{array}
\end{eqnarray}
where $\Gamma$ and $\beta$ are determined by
\begin{equation}
\label{eq:rate_damp}
\Gamma=\Gamma_e+\frac{\gamma_m}{4}+\frac{\kappa_f^2}{\gamma_m}+\kappa_f
\end{equation}
\begin{equation}
\label{eq:rate_pump}
e^{-\beta\omega_o}\Gamma=\Gamma_e e^{-\beta_e^-\omega_o}+\frac{\gamma_m}{4}+\frac{\kappa_f^2}{\gamma_m}-\kappa_f.
\end{equation}
\par
The effective temperature $1/\beta$ becomes positive by setting the feedback strength above the following threshold:
\begin{equation}
\label{eq:kappa_thresh}
\kappa_f > -\kappa_e^-=\frac{1}{2}\Gamma_e\left(e^{-\beta_e^-\omega_o}-1\right).
\end{equation}
To summarize, as a result of the feedback the negative temperature
heat bath and the negative amplitude damping rate $\kappa_e^-$ for HO, becomes an effective positive temperature heat bath with amplitude damping rate $\kappa_f+\kappa_e^->0$.
\section{Steady State and Work Extraction}
For sufficient strong feedback $\kappa_f$, satisfying $\kappa_f+\kappa_e^->0$, the master equation \eqref{eq:ME_final} is a standard thermalizing master equation with resonant external driving.
It has a unique stationary state which in the Schr\"{o}dinger picture is a thermal state with rotating displacement (see Appendix D),
\begin{equation}
\label{eq:steady-Sch}
\hat\rho^{\infty}\propto\exp\left[-\beta\omega_o\left(\hat{c}-c_\infty e^{-i\omega_o t}\right)^{\dagger}\left(\hat{c}-c_\infty e^{-i\omega_o t}\right)\right].
\end{equation}
where
\begin{equation}
\label{eq:c_infty}
c_\infty=-\frac{\epsilon_d}{\kappa_f+\kappa_e^-}<0.
\end{equation}
Hence, the mean amplitude rotates, $\ave{\hat{c}}_\infty=c_\infty e^{-i\omega_o t}$, its phase is shifted by $-\pi/2$ with respect to the external driving .
The average population is given by the sum of the Bose statistic  $n_o$ and the yield of displacement
\begin{equation}
\label{eq:n_infty}
\ave{\hat{c}^{\dagger}\hat{c}}_\infty =  \frac{1}{e^{\beta\omega_o}-1} +\vert c_\infty\vert^2 \equiv n_o+\vert c_\infty\vert^2.
\end{equation}
\par
We distinguish two opposing regimes of the steady-state operation of the flywheel.
The first is the deep quantum regime, $n_o,\vert c_\infty\vert^2\ll1$, where the flywheel is operating in the vicinity of its ground state.
The second is the classical regime in which both the thermal occupation and the displacement are large  numbers, $n_o,\vert c_\infty\vert^2\gg1$.
The two crossed regimes also present peculiar quantum features.
Recall that weak coupling condition \eqref{eq:weak_coupling} sets an
asymptotic upper limit on the total occupation (\ref{eq:n_infty}).
This implies asymptotic upper limits on the temperature $1/\beta$, excluding too high thermal occupations $n_o$, as well as on the driving strength $\epsilon_d$, confining the displacements $\ave{\hat{c}}_\infty$.
Thus, accessibility to the classical regime depends on the physical properties of the two-qubit heat engine and its coupling $g$ to the flywheel.
\begin{figure}[htbp]
\center{\includegraphics[width=8.6cm]{finalfig.png}}
\caption{ Charging efficiency as function of measurement strength $\gamma_m$ and feedback strength $\kappa_f$.   The percentage of useful work out of the entire energy stored in the flywheel has a maximum for the ratio $\gamma_m / \kappa_f = 2$, and it is further maximised for $\kappa_f$ approaching its threshold $\vert\kappa_e^-\vert$.
Here: $\omega_o=1$, $\beta_e^-=-10^{-1}$, $\Gamma_e=10^{-6}$ and $\epsilon_d=9\cdot10^{-2}$.}
\label{fig:Charging efficiency}
\end{figure}
The steady  state \eqref{eq:steady-Sch} becomes a displaced Gibbs state and as such, it is suitable for work extraction.
The internal energy of the steady-state is given by
$\mathcal{E} =  \omega_o(n_o+\vert c_\infty\vert^2)$.
Applying a unitary displacement transformation can bring the state eq.(\ref{eq:steady-Sch}) into a Gibbs state (passive state) with the temperature $1/\beta$. Thus, the part of the internal energy that is due to $c_\infty$ can all be extracted by the unitary operation as the maximum useful work
\begin{equation}
\label{eq:Work}
\mathcal{W}=\omega_o\vert c_\infty\vert^2=\frac{\omega_o\epsilon_d^2}{(\kappa_f+\kappa_e^-)^2},
\end{equation}
which is independent of the strength $\gamma_m$ of the monitoring.
The \emph{charging efficiency} of the flywheel can be defined as the ratio between useful work and the internal energy stored in the HO (see Fig.\ref{fig:Charging efficiency}),
\begin{equation}
\label{eq:Wefficiency}
\eta=\frac{\mathcal{W}}{\mathcal{E}}=\frac{1}{1+n_o/\vert c_\infty\vert^2}.
\end{equation}
The efficiency is improved for small thermal occupation $n_o$ and large displacement $c_\infty$.
The occupation $n_o$ becomes small when the effective temperature $1/\beta$ is reduced.
Interestingly, this singles out the optimum measurement strength $\gamma_m$ which has so far remained unconstrained.
From eqs. (\ref{eq:rate_damp},\ref{eq:rate_pump}) we find that $1/\beta$ takes its minimum value with the choice $\gamma_m=2\kappa_f$ obtaining the expression
\begin{equation}
n_o\vert_{min}=\frac{1}{(1+2\kappa_f/\Gamma_e)e^{\beta_e^-\omega_o}-1}
\end{equation}
of the minimum thermal occupation. The charging efficiency yields:
\begin{equation}
\eta\vert_{\gamma_m=2\kappa_f}
=\frac{1}{1+\frac{\Gamma_e}{2\epsilon_d^2}e^{-\beta_e^-\omega_o}(\kappa_f+\kappa_e^-)^2}.
\end{equation}
The efficiency $\eta_{\gamma_m=2\kappa_f}$ together with the extractable work $\mathcal{W}$ reach higher values if we increase the displacement $\vert c_\infty\vert$.
In particular, the efficiency approaches its maximal value $1$
when the feedback $\kappa_f$ approaches its lower threshold, $\kappa_f \rightarrow -\kappa_e^-$.
A different technique to maximize both the efficiency and the work is by increasing $\epsilon_d$, i.e., applying a stronger driving field.
Nevertheless, as was already mentioned, these two approaches are limited by the weak-coupling condition (\ref{eq:weak_coupling}).
\section{Energy Flows in Steady-State}
A macroscopic flywheel at rest requires an input work (initial push) to reach the vicinity of steady-state, at which point the output power is larger than the input power.
Regulating the flywheel also has energetic costs that should be accounted for.
These energetic considerations, in principle, also apply to the quantum flywheel.
However, the related calculations would need a novel approach to  heat flow and power in quantum systems under stochastic control.
\par
The standard definition of thermodynamic heat flow $\mathcal{J}$ and power
$\mathcal{P}$ in open quantum systems is given \cite{alicki79} by the time derivative of the internal energy $\mathcal{E}=\mathrm{tr}[\hat\rho \hat H]$ in the following manner:
\begin{equation}
\label{eq:dE}
d\mathcal{E}=\mathrm{tr}[d\hat\rho \hat H]+\mathrm{tr}[\hat\rho d\hat H]\equiv\mathcal{J} dt+\mathcal{P} dt.
\end{equation}
The Hamiltonian and the state of the system are typically stochastic in the theory of monitoring and feedback control.
Since stochastic fluctuations are microscopic, the thermodynamic definition of the internal energy is given by the stochastic mean  of the microscopic energy, $\mathbf{M}\mathrm{tr}(\hat{\sigma}\hat H)$.
This leads to the following generalization of the standard thermodynamic relation:
 \begin{equation}
 \label{eq:dE_stoch}
d\mathcal{E}=\mathbf{M}\mathrm{tr}[d\hat{\sigma} \hat H]+\mathbf{M}\mathrm{tr}[\hat{\sigma} d\hat H]\equiv\mathcal{J} dt+\mathcal{P} dt.
\end{equation}
The differentials in (\ref{eq:dE_stoch}) must be Stratonovich ones instead of those of Ito.
For the Ito differentials the rhs should contain the so-called Ito correction  $\mathbf{M}\mathrm{tr}[d\hat{\sigma} d\hat H]$ which would jeopardize the split of $d\mathcal{E}$ between heat flow and power.
\par
Any systematic calculation of heat flow and power requires to transform the final SME from Ito into Stratonovich form.
We postpone this very novel task to future research.
Rather, we focus on the minimal calculations and considerations confirming that our model represents a genuine heat engine.
\par
Next, we show that there is a consumable output power in the steady state operation of the flywheel.
The total Hamiltonian has two time-dependent contributions $H_d(t)$ and $H_f(t)$.
Accordingly, the power $\mathcal{P}$ consists of two contributions corresponding to the power  invested by the driving  and the power  gained from the feedback.
The first, in the steady state $\mathbf{M}\hat{\sigma}=\hat\rho^\infty$ reads
\begin{eqnarray}
\begin{array}{lr}
\label{eq:P_d}
\mathcal{P}_d=\mathbf{M}\mathrm{tr}\left[\hat{\sigma}\frac{d\hat H_d}{dt}\right]
=\\
\mathrm{tr}\left[\hat\rho^\infty\frac{d}{dt}(-i\epsilon_d\hat{c}^{\dagger}e^{-i\omega_o t}
+\mathrm{H.c.})\right]
=-2\epsilon_d \omega_o c_{\infty}>0,

\end{array}
\end{eqnarray}
where the positivity  indicates power going into (consumed by) the flywheel.
We restrict our calculations for the deterministic part of feedback, i.e., we replace $\hat H_f(t)$ by its deterministic part
$\hat H_{f,det}=-i\kappa_f\ave{\hat{c}}_\sigma\hat{c}^\dagger+\mathrm{H.c.}$.
As was mentioned before, considering the stochastic part
$\hat H_{f,sto}=-i\kappa_f/\sqrt{\gamma_m})\hat{c}^\dagger d\xi+\mathrm{H.c.}$ requires the Stratonovich calculus.
The power reads
\begin{widetext}
\begin{equation}
\label{eq:p_fdet}
\mathcal{P}_{f,det}=\mathbf{M}\mathrm{tr}\left[\hat{\sigma}\frac{d\hat H_{f,det}}{dt}\right]
 =
  \mathbf{M}\mathrm{tr}\left[\hat{\sigma}\frac{d}{dt}(-\kappa_f\ave{\hat{c}}_\sigma\hat{c}^\dagger+\mathrm{H.c.})\right]
=
-i\kappa_f\mathbf{M}\mathrm{tr}\left[\frac{d\hat{\sigma}}{dt}\hat{c}\right]\ave{\hat{c}^\dagger}_\sigma+\mathrm{c.c}.
\end{equation}
\end{widetext}

The power in eq.(\ref{eq:p_fdet}) is proportional to the (weighted) mean of the phase drift of the amplitude $\ave{\hat{c}}_\sigma$.
To calculate $d\hat{\sigma}$ we apply the final SME given in  Appendix D.
The only relevant yield is the unitary rotation $-i\omega_o\ave{\hat{c}}_\sigma dt$ since the dissipative part does not alter the phase of   $\ave{\hat{c}}_\sigma$
and the  Ito stochastic part will cancel out by the mean operation $\mathbf{M}$. Therefore we get
\begin{equation}
\label{eq:P_f}
\mathcal{P}_{f,det}=-2\kappa_f \omega_o \mathbf{M}|\ave{\hat{c}}_{\sigma}|^2<0.
\end{equation}
Negativity means that power is gained (supplied) by feedback.
Although analytical solutions for similar SMEs such as ours exist \cite{gisin91}, we restrict ourselves to a simple guess.
Using  the Cauchy-Schwartz relation $\mathbf{M}\vert\ave{\hat{c}}_\sigma
\vert^2 \geq \vert\mathbf{M}\ave{\hat{c}}_\sigma\vert^2$, we obtain the lower bound
$-\mathcal{P}_{f,det}\geq 2\kappa_f\omega_o |c_{\infty}|^2$
for the stationary power gained by feedback in the steady-state.
Hence the overall stationary power satisfies the inequality
\begin{equation}
\label{eq:P_lowerbound}
-\mathcal{P}_{det}=-\mathcal{P}_d-\mathcal{P}_{f,det} \geq 2\omega_o\epsilon_d^2\frac{-\kappa_e^-}{(\kappa_f+\kappa_e^-)^2}.
\end{equation}
The sign is negative and thus the consumable power of the flywheel is positive and bounded from below.
We conjecture that the contribution of the stochastic part $\hat H_{f,sto}(t)$ of driving can not invalidate the positivity of the consumable power.
\par
We summarize the plausible structure of energy currents (see Fig.\ref{fig:schem}d).
The steady state energy balance contains five different currents:
\begin{equation}
\dot\mathcal{E} = \mathcal{J}_e + \mathcal{J}_m + \mathcal{J}_f + \mathcal{P}_d + \mathcal{P}_f=0.
\end{equation}

The heat flowing into the flywheel has two contributions,  the first is from the engine, $\mathcal{J}_e$, the second is from the monitoring device,  $\mathcal{J}_m$.
Power from the driving field, $\mathcal{P}_d$, is also consumed by the flywheel, and serves as an input power activating the flywheel.
This power is overcompensated by the output power $\mathcal{P}_f$ realized by the feedback.
In addition, the outflow $\mathcal{J}_f$ cools the flywheel, thereby stabilizing it and lowering the entropy produced in the flywheel as a result of the engine and the monitoring operations.
In the case $\beta_e^- \rightarrow0^-$ of no population inversion in the engine, the heat flow $\mathcal{J}_e$ and the consumable power must vanish.
The work \eqref{eq:Work} stored in the flywheel reaches its minimal, yet positive, value $\mathcal{W}=\omega_o \epsilon_d^2/\kappa_f^2$.
\section{Summary}
Population inversion, corresponding to negative temperature $1/\beta_e^-$ in a few-level quantum heat engine has been established a long time ago \cite{ramsey56}.
In this paper we have shown that the heat engine operation is equivalent to a negative temperature heat bath in the standard dynamical sense.
Thus, its influence on the work repository is the typical thermalizing master equation extended to negative temperature $1/\beta_e^-$.
\par
Work extraction  is still an outstanding issue because of the spread of thermal and quantum noise over the work repository which in our case is a quantum HO.
If the HO is replaced by an idealized classical field, all the energy flowing out of the engine can in principle be extracted as power (see Appendix E).
This results from the fact that classical field generates a unitary operation on the system which does not change its entropy.
However, when a work repository is quantized the work exchange is accompanied with heat exchange which may degrade the charging efficiency.
In this paper we introduced a generic novel approach that can be applied to resolve such problems.
Specifically, we demonstrated the difficulties of storing useful work in a quantum harmonic oscillator.
Overcoming the unlimited growth of fluctuations, regulating and stabilizing the flywheel is achieved by applying monitoring and feedback control to the system.
\par
The steady-state, the power, and the stored extractable energy of the flywheel are determined analytically.
While the amount of work stored in the flywheel is  independent of the accuracy of the monitoring, the charging efficiency is optimized for a special ratio between the monitoring and the feedback strength.
In other words there is a balance between information gained by monitoring the flywheel and the information fed back to the flywheel.
The balance coincides with minimum temperature of the flywheel.
Breaking this balance implies that the phase space distribution is no longer optimal for work extraction from the flywheel.
\par
This model is a prototype of an analytically tractable model of a quantum heat engine coupled to a single degree-of-freedom work repository, operating continuously in steady state under quantum control.
Experiments which employ quantum monitoring and feedback strategies are becoming common in recent years \cite{,clerk08,mabuchi04,vijay12,nori14}.
Future advances in quantum technologies depend on our ability to control and manipulate quantum system.
A firm theoretical foundation relating systems that are subject to quantum monitoring and feedback control with basic concepts of thermodynamics is still missing.
\section{Acknowledgments}
We thank Raam Uzdin, Saar Rahav, David Gelbwaser-Klimovsky and Walter Singaram  for fruitful discussions. This work was supported by the Israeli Science Foundation and the Hungarian Scientific Research Fund under Grant No. 103917.
Part of this work was supported by the COST Action MP1209 ``Thermodynamics in the quantum
regime''.
\par
\section*{Appendix A: Tripartite Heat Engine}
\label{sec:a}
We use an interaction picture  for its convenience especially for our master equations.
The stochastic master equations of monitoring and feedback are presented in the Schr\"{o}dinger picture
for transparency. Heat flow and power are, as a rule, defined in the Schr\"{o}dinger picture.
We derive the master equation for the harmonic oscillator (HO) subject to the operation of the engine.
The quantum heat engine is comprised of two two-level systems (TLS's), with the Hamiltonians $\hat H_h = \omega_h \hat{a}^{\dagger}\hat{a}$ and $\hat H_c= \omega_c \hat{b}^{\dagger} \hat{b}$.
 The two TLS's are coupled to a hot and a cold heat bath, respectively, at temperatures $T_h>T_c$.
The dynamics follow the Lindblad-Gorini-Kossakowski-Sudarshan dynamics \cite{lindblad76,gorini276},
in the interaction picture of $\hat H_{h(c)}$ the corresponding master equations read
\begin{widetext}
\begin{eqnarray}
\label{eq:LGKS}
\begin{array}{ll}
 \frac{d\hat\rho_h}{dt}=
\Gamma_h\left(\hat{a}\hat\rho_h\hat{a}^\dagger-\mathcal{H} \hat{a}^\dagger\hat{a}\hat\rho_h\ +e^{-\beta_h\omega_h}(\hat{a}^\dagger\hat\rho_h\hat{a}- \mathcal{H}  \hat{a}\aop^\dagger\hat\rho_h)\right)\equiv\mathcal{L}_h\hat\rho_h\\
 \frac{d\hat\rho_c}{dt}=
\Gamma_c\left(\hat{b}\hat\rho_c\hat{b}^\dagger-\mathcal{H}\hat{b}^\dagger\hat{b}\hat\rho_c\ +e^{-\beta_c\omega_c}(\hat{b}^\dagger\hat\rho_c\hat{b}- \mathcal{H}\hat{b} \hat{b}^\dagger\hat\rho_c)\right)\equiv\mathcal{L}_c\hat\rho_c\\
\end{array}
\end{eqnarray}
\end{widetext}
where $\Gamma_{h(c)}$ are the damping rates.
(In our convention, different from  ref.'s \cite{wisemanbook}, $\mathcal{H}$ denotes the Hermitian part of all what stands after it.)
The heat baths bring the TLS's to thermal equilibrium states $\hat\rho_{h(c)}^\infty$
with the occupation numbers $n_{h(c)}=1/(e^{\beta_{h(c)}\omega_{h(c)}}+1)$, and with the inverse temperatures $\beta_{h(c)}=1/T_{h(c)}$, respectively.

The two TLS's are then weekly coupled to a quantum HO of self-Hamiltonian $\hat H_o=\omega_o\hat{c}^{\dagger}\hat{c}$, via the tripartite Hamiltonian
\begin{equation}
\label{eq:Hint}
 \hat K =-ig\hat{a}\hat{b}^{\dagger}\hat{c}^{\dagger} + H.c.
\end{equation}
We work in resonance, $\omega_o=\omega_h-\omega_c$, and in the weak coupling regime for which a local master equation holds \cite{levy14b}.
The master equation in the interaction picture for the tripartite state $\hat\rho_3$ of the TLS's coupled to the HO is written as
\begin{equation}
\label{eq:master1}
 \frac{d\hat\rho_3}{dt} = (\mathcal{L} + \mathcal{K})\hat\rho_3,
\end{equation}
with  $\mathcal{L}=\mathcal{L}_h + \mathcal{L}_c$ and
\begin{equation}
\mathcal{K} \hat\rho_3 = -i[\hat K ,\hat\rho_3].
\end{equation}
We will to derive the effective master equation for the HO's state $\hat\rho$ assuming that the TLS's are initially in their equilibrium states
$\hat\rho_{hc}^\infty=\hat\rho_h^\infty\otimes\hat\rho_c^\infty$ and the initial state of the tripartite system is the product state $\hat\rho_3(0)=\hat\rho_{hc}^\infty\otimes\hat\rho(0)$.
The solution of the master equation  \eqref{eq:master1} can be written
in the implicit form
\begin{equation}
 \hat\rho_3(t)=\hat\rho_3(0) + \int^{t}_0 ds  e^{\mathcal{L}(t-s)}\mathcal{K} \hat\rho_3(s),
\end{equation}
which we can confirm by taking the time derivative of both sides of the equation, and using the relation $\mathcal{L}\hat\rho_3(0)=0$.
Inserting the above solution into the rhs  of eq.\eqref{eq:master1}, we obtain
\begin{equation}
 \frac{d \hat\rho_3(t)}{dt}=\mathcal{K}\hat\rho_3(0) + (\mathcal{L}+\mathcal{K})\int^{t}_0 ds  e^{\mathcal{L}(t-s)}\mathcal{K} \hat\rho_3(s).
\end{equation}
We assume that $\hat\rho_3(s)\approx \hat\rho_{hc}^\infty\otimes\hat\rho(s)$.
This assumption is justified when the thermalization time of the TLS's is faster than the time scale in which the system is changed significantly due to coupling (\ref{eq:Hint}).
Taking the partial trace over the TLS's:
\begin{equation}
 \frac{d \hat\rho(t)}{dt}=\mathrm{tr}_{hc}[\mathcal{K}\int^{t}_0 ds  e^{\mathcal{L}(t-s)}\mathcal{K} \hat\rho_{hc}^\infty\otimes\hat\rho(s)].
\end{equation}
Here we have used the relations
$\mathrm{tr}_{hc}[\mathcal{K}\hat\rho_{hc}^\infty ]=0$ and $\mathrm{tr}_{hc}[\mathcal{L}\int_0^t e^{\mathcal{L}(t-s)} \mathcal{K} \hat\rho_{hc}^\infty]=0$.
Performing the standard Markovian approximations \cite{breuer} we obtain
\begin{equation}
 \frac{d \hat\rho(t)}{dt}=\mathrm{tr}_{hc}[\mathcal{K}\int^{\infty}_0 ds  e^{\mathcal{L}s}\mathcal{K} \hat\rho_{hc}^\infty\otimes\hat\rho(t)],
\end{equation}
which can be written explicitly as
\begin{eqnarray}
 \frac{d \hat\rho}{dt}=-\mathrm{tr}_{hc}\left[\hat K,\int^{\infty}_0 ds  e^{\mathcal{L}s}\left[\hat K , \hat\rho_{hc}^\infty \otimes \hat\rho\right]\right]
\\ \nonumber
=-\mathrm{tr}_{hc}\int^{\infty}_0 ds \left[\left(e^{\mathcal{L^{\dagger}}s}\hat K \right), \left[ \hat K, \hat\rho_{hc}^\infty\otimes\hat\rho\right]\right].
\end{eqnarray}
Making use of the relation
\begin{equation}
 e^{\mathcal{L^{\dagger}}s}\hat K  = \hat K\exp[-\frac{1}{2} \sum_{l=h,c}\Gamma_l(1+e^{-\beta_l \omega_l})s],
\end{equation}
we have
 %
\begin{widetext}
\begin{equation}
\frac{d\hat\rho}{dt}= \frac{(2g)^2}{\sum_{l=h,c} \Gamma_l(1+ e^{-\beta_l \omega_l})}
\left[\langle\hat{a}\aop^\dagger\rangle\!_\infty\langle\hat{b}^\dagger\hat{b}\rangle\!_\infty(\hat{c}\hat\rho \hat{c}^\dagger-\mathcal{H}\hat{c}^\dagger\hat{c}\hat\rho)
 +\langle\hat{a}^\dagger\hat{a}\rangle\!_\infty\langle\hat{b}\bop^\dagger\rangle\!_\infty(\hat{c}^\dagger\hat\rho \hat{c}-\mathcal{H}\hat{c}\cop^\dagger\hat\rho) \right]
\end{equation}
\end{widetext}
where $\ave{\cdot}_\infty$ stands for the expectation value with respect to the TLS's thermal equilibrium states $\hat\rho_{h(c)}^\infty$.
Finally, the master equation for the HO subject to the engine operation takes the form
\begin{equation}
\label{eq:L_e}
\frac{d\hat\rho}{dt}\equiv\mathcal{L}_e\hat\rho=
\Gamma_e(\hat{c}\hat\rho\hat{c}^\dagger-\mathcal{H}\hat{c}^\dagger\hat{c}\hat\rho)
      +\Gamma_e e^{-\beta_e^-\omega_o}(\hat{c}^\dagger\hat\rho\hat{c}-\mathcal{H}\hat{c}\cop^\dagger\hat\rho),
\end{equation}
where
\begin{equation}
\Gamma_e=(2g)^2\frac{(1-n_h)^2(1-n_c)n_c}{\Gamma_h(1-n_c)+\Gamma_c(1-n_h)},
\end{equation}
and the output temperature of the heat-engine
\begin{equation}
\beta_e^-=\frac{\beta_h\omega_h - \beta_c\omega_c}{\omega_h-\omega_c}
\end{equation}
which is a function of the TLS's excitation energies and temperatures only.
We operate the system as a heat engine, i.e.  $T_h/T_c > \omega_h/\omega_c >1$,
the effective temperature is negative, i.e. $1/\beta_e^- < 0$, and the HO will not reach a stable asymptotic state, as we show below.
The master equation \eqref{eq:L_e} together with the Hamiltonian $\hat H_o$ yield closed evolution equations for the mean
amplitude $\ave{\hat{c}}_t$ as well as for the occupation $\ave{\hat{c}^\dagger\hat{c}}_t$:
\begin{eqnarray}
\label{eq:c_t}
\frac{d\ave{\hat{c}}_t}{dt}&=&-\left(\kappa_e^-+i\omega_o\right)\ave{\hat{c}}_t\\
\label{eq:n_t}
\frac{d\ave{\hat{c}^\dagger\hat{c}}_t}{dt}&=&-2\kappa_e^-\ave{\hat{c}^\dagger\hat{c}}_t+\Gamma_e e^{-\beta_e^-\omega_o}
\end{eqnarray}
where
\begin{equation}
\kappa_e^- = \frac{1}{2}\Gamma_e\left(1-e^{-\beta_e^-\omega_o}\right)<0
\end{equation}
is the standard amplitude damping constant. This time it is negative
since $\beta_e^-<0$ therefore both $\ave{\hat{c}}_t$ and $\ave{\hat{c}^\dagger\hat{c}}_t$ diverge exponentially with time.
In particular, a thermal state remains thermal, the temperature is increasing
exponentially as can be shown by the simple solution of \eqref{eq:n_t} for
the occupation. Note, however, that our model is only valid in the weak coupling regime where the thermalization time is shorter than the internal time scale. This implies that the occupation must be limited by
 %
\begin{equation}
\label{eq:weak}
g\sqrt{\ave{\hat{c}^\dagger\hat{c}}+1}\ll\Gamma_{h(c)} \left(1+e^{-\beta_{h(c)}\omega_{h(c)}}\right).
\end{equation}

\section*{Appendix B: External Driving}
Coupling the HO to a resonant oscillating external field. Via such \emph{driving} one would expect to extract power.
Consider the time-dependent Hamiltonian in the Schr\"{o}dinger picture,
\begin{equation}
\label{eq:Hext}
\hat H_d(t)=-i\epsilon_d\left(\hat{c}^{\dagger}e^{-i\omega_o t} -\hat{c} e^{i\omega_o t}\right)
\end{equation}
where $\epsilon_d>0$.
In the interaction picture, the master equation \eqref{eq:L_e} is modified by an additional static Hamiltonian:
\begin{equation}
\label{eq:L_e_ext}
\frac{d\hat\rho}{dt}=\mathcal{L}_e\hat\rho-\epsilon_d[\hat{c}^{\dagger} -\hat{c},\hat\rho].
\end{equation}
Now the rhs of eq. \eqref{eq:c_t} of the mean amplitude acquires an additional term $-\epsilon_d e^{-i\omega_o t}$.
This allows for an exceptional  stationary solution of constant amplitude with the rotating phase:
\begin{equation}
\label{eq:c-const}
\ave{\hat{c}}_t=-\frac{\epsilon_d}{\kappa_e^-}e^{-i\omega_o t}=\mathrm{const}\times e^{-i\omega_o t}.
\end{equation}
This solution is unstable since all neighbouring solutions exponentially diverge with $t$.
As to the occupation $\ave{\hat{c}^\dagger\hat{c}}_t$,
 the rhs of eq.\eqref{eq:n_t} acquires the additional linear term
$-\epsilon_d(\ave{\hat{c}^\dagger}_t-\ave{\hat{c}}_t)$ hence the occupation remains exponentially divergent,
there is no steady state solution under external driving.
The stability issue of the HO is still not resolved.
\section*{Appendix C: Monitoring}
Continuous measurement, i.e.: \emph{monitoring}, is  the first task towards feedback control on the system \cite{wisemanbook}.
Here we consider the time-continuous measurement of both quadratures
$\hat{x}=\frac{1}{\sqrt{2}}(\hat{c}^\dagger +\hat{c})$ and $\hat{y}=\frac{i}{\sqrt{2}}(\hat{c}^\dagger -\hat{c})$ of the HO.
Generalizing the result of \cite{diosi88} for monitoring simultaneously $\hat{x}$ and $\hat{y}$,
we can write the following stochastic master equation (SME) in the Schr\"{o}dinger picture for the density matrix $\hat{\sigma}$ conditioned on both measurement signals $\bar{x},\bar{y}$:
\begin{widetext}
\begin{eqnarray}
\label{eq:SME}
d\hat{\sigma}&=&-i\left[\hat H_o,\hat{\sigma}\right]dt-\frac{\gamma_m}{8}\left[\hat{x},\left[\hat{x},\hat{\sigma}\right]\right]dt-\frac{\gamma_m}{8}\left[\hat{y},\left[\hat{y},\hat{\sigma}\right]\right]dt
+ \mathcal{H}\sqrt{\gamma_m}\left(\hat{x}-\left\langle\hat{x}\right\rangle_\sigma\right)\hat{\sigma} d\xi_{x}
     + \mathcal{H}\sqrt{\gamma_m}\left(\hat{y}-\left\langle\hat{y}\right\rangle_\sigma\right)\hat{\sigma} d\xi_{y}.
\end{eqnarray}
\end{widetext}
All expectation values $\ave{\cdot}_\sigma$ are understood in the stochastic conditional state $\hat{\sigma}$.
The measurement signals satisfy
\begin{equation}
\label{eq:signals_xy}
\bar{x}dt=\ave{\hat{x}}_\sigma dt+\frac{d\xi_x}{\sqrt{\gamma_m}},~~~~
\bar{y}dt=\ave{\hat{y}}_\sigma dt+\frac{d\xi_y}{\sqrt{\gamma_m}}.
\end{equation}
Here $d\xi_x,d\xi_y$ are Ito increments of independent standard Wiener processes, satisfying
\begin{equation}
(d\xi_x)^2 =(d\xi_y)^2 = dt,~~d\xi_x d\xi_y=0,~~
\mathbf{M} d\xi_x=\mathbf{M} d\xi_y=0,
\end{equation}
with the symbol $\mathbf{M}$ for stochastic mean, and $\gamma_m$ for the measurement  strength.
(Note that we changed $\gamma_m$ in ref. \cite{diosi88} for $\gamma_m/2)$.
We can return to complex notation, i.e.,  we rewrite the above equations in terms of $\hat{c},\hat{c}^\dagger$ and the corresponding complex signal $\bar{c}=(\bar{x}+i\bar{y})/\sqrt{2}$.
We define the complex Wiener increment
\begin{equation}
d\xi=\frac{d\xi_x+i d\xi_y}{\sqrt{2}}
\end{equation}
that satisfies
\begin{equation}
(d\xi)^2=(d\xi^\ast)^2=0,~~~~~d\xi^\ast d\xi=dt,~~~~~\mathbf{M} d\xi=\mathbf{M} d\xi^\ast=0.
\end{equation}
The SME \eqref{eq:SME} of the conditional state becomes
\begin{widetext}
\begin{eqnarray}
d\hat{\sigma}&=&-i\left[\hat H_o,\hat{\sigma}\right]dt
+\frac{\gamma_m}{4}\left(\hat{c}\hat\rho\hat{c}^\dagger-\mathcal{H}\hat{c}^\dagger\hat{c}\hat{\sigma}
                                            +\hat{c}^\dagger\hat\rho\hat{c}-\mathcal{H}\hat{c}\cop^\dagger\hat{\sigma}\right)
+ \sqrt{\gamma_m}\mathcal{H}\left(\left(\hat{c}-\ave{\hat{c}}_\sigma\right)d\xi^{*}+\mathrm{H.c.}\right)\hat{\sigma}\nonumber\\
&\equiv&-i\left[\hat H_o,\hat{\sigma}\right]dt +\mathcal{L}_m\hat{\sigma} dt
+ \sqrt{\gamma_m}\mathcal{H}\left(\left(\hat{c}-\ave{\hat{c}}_\sigma\right)d\xi^{*}+\mathrm{H.c.}\right)\hat{\sigma}
\end{eqnarray}
\end{widetext}
Eqs.\eqref{eq:signals_xy} of the real signals take the following form for the complex signal:
\begin{equation}
\label{eq:signal_c}
\bar{c}dt=\ave{\hat{c}}_\sigma dt+\frac{d\xi}{\sqrt{\gamma_m}}.
\end{equation}
Applying this time-continuous measurement to the HO which is coupled to the heat engine
and driven by the external field,  cf. eq. \eqref{eq:L_e_ext}, we get the following SME:
\begin{widetext}
\begin{equation}
\label{eq:smemonitor}
d\hat{\sigma}=-i[\hat H_o,\hat{\sigma}]dt+(\mathcal{L}_e+\mathcal{L}_m)\hat{\sigma} dt-\epsilon_d[\hat{c}^\dagger e^{i\omega_o t}-\hat{c},\hat{\sigma}]dt
+\sqrt{\gamma_m} \mathcal{H}[(\hat{c}-\ave{\hat{c}}_\sigma)d\xi^\ast+\mathrm{H.c.}]\hat{\sigma}.
\end{equation}
\end{widetext}
The state $\hat{\sigma}$ of the HO is the conditioned state on the measured signal \eqref{eq:signal_c},
its stochastic mean is the unconditional density matrix: $\mathbf{M}\hat{\sigma}=\hat\rho$.
Taking the stochastic mean $\mathbf{M}$ of both sides of the SME,
we are left with the master equation of the unconditional state:
\begin{equation}
\frac{d\hat\rho}{dt}= (\mathcal{L}_e+\mathcal{L}_m)\hat\rho-\epsilon_d[\hat{c}^\dagger-\hat{c},\hat\rho].
\end{equation}
As a result of the measurement, additional heat flows into the oscillator, the damping
rate becomes  $\Gamma_e+\gamma_m$ and the inverse "temperature" $\beta_e^-$ is
modifies but remains negative.
The exceptional steady amplitude \eqref{eq:c-const} exists with the
modified parameters, but it is unstable like all other solutions.
\section*{Appendix D: Feedback Control}
Using the measured signal \eqref{eq:signal_c}, we control the state of the HO in the vicinity of the constant rotating
amplitude set by the external driving, in such a way that we get a true stable
steady state. Consider the following \emph{feedback} Hamiltonian in the Schr\"{o}dinger picture:
\begin{equation}
\hat H_f(t)=-i\kappa_f\bar{c}(t)\hat{c}^\dagger +\mathrm{H.c.}.
\end{equation}
Here $\kappa_f$ is the feedback strength.
We apply the feedback \cite{diosi94} on top of the monitored evolution described by eq.(\ref{eq:smemonitor}):
\begin{equation}
\label{eq:drho}
\hat{\sigma}+d\hat{\sigma}\rightarrow e^{-i\hat H_f dt}(\hat{\sigma}+d\hat{\sigma}) e^{i\hat H_f dt}.
\end{equation}
Expanding the rhs into a series, keeping first order terms in $dt $ and keeping in mind that  $\vert d\xi\vert^2=dt$,
the terms that are left for evaluation are: $-i\left[H_fdt,\hat{\sigma}\right]$, $-i\left[H_fdt,g\hat{\sigma}\right]$ and $-\frac{1}{2} \left[H_fdt,\left[H_fdt,\hat{\sigma}\right]\right]$.
The final SME including feedback reads:
\begin{widetext}
\begin{equation}
\label{eq:smecontrol}
d\hat{\sigma}  =  -i[\hat H_o,\hat{\sigma}]dt+(\mathcal{L}_e+\mathcal{L}_m+\mathcal{L}_f)\hat{\sigma} dt-\epsilon_d\left[\hat{c}^\dagger e^{i\omega_o t}-\hat{c} e^{-i\omega_o t},\hat{\sigma}\right]dt
  -\frac{\kappa_f}{\sqrt{\gamma_m}}\left[\hat{c}^\dagger d\xi-\hat{c} d\xi^{*},\hat{\sigma}\right]
         +\sqrt{\gamma_m}\mathcal{H}\left[\hat{c}-\ave{\hat{c}}_\sigma)d\xi^{*}+\mathrm{H.c.}\right]\hat{\sigma}.
\end{equation}
\end{widetext}
The dissipative contribution of the feedback reads:
\begin{eqnarray}
\label{eq:L_f}
\mathcal{L}_f\hat{\sigma}=
 \left(\frac{\kappa_f^2}{\gamma_m}+\kappa_f\right)\left(\hat{c}\hat{\sigma}\hat{c}^\dagger-\mathcal{H}\hat{c}^\dagger\hat{c}\hat{\sigma}\right)
\\ \nonumber
+\left(\frac{\kappa_f^2}{\gamma_m}-\kappa_f\right)\left(\hat{c}^\dagger\hat{\sigma}\hat{c}-\mathcal{H}\hat{c}\cop^\dagger\hat{\sigma}\right).
\end{eqnarray}
For $\kappa_f>\gamma_m$ this corresponds to a thermal bath of positive temperature.  Entering the regime $0<\kappa_f<\gamma_m$,
the cooling effect of $\mathcal{L}_f$ within the sum  $\mathcal{L}_e+\mathcal{L}_m+\mathcal{L}_f$ becomes enhanced although
$\mathcal{L}_f$ ceases to be a mathematically correct dissipator in itself.
Taking the stochastic mean over eq. (\ref{eq:smecontrol}) we obtain the master equation of the unconditional state
which in the interaction picture takes this form :
\begin{equation}
\label{eq:me}
\frac{d\hat\rho}{dt}=(\mathcal{L}_e+\mathcal{L}_m+\mathcal{L}_f)\hat\rho-\epsilon_d\left[\hat{c}^{\dagger}-\hat{c},\hat\rho\right].
\end{equation}

What we have got for the HO dynamics is the following:
The HO is excited by the negative temperature ($1/\beta_e^-$) bath $\mathcal{L}_e$ due to population inversion, heated by the infinite temperature bath  $\mathcal{L}_m$ due
to noise of monitoring, and cooled by the feedback $\mathcal{L}_f$.
On top of this, the external driving shifts the Hamiltonian $\hat H_o$.
We write the full master equation \eqref{eq:me} in a compact form:
\begin{eqnarray}
\begin{array}{lr}
\label{eq:me_final}
\frac{d\hat\rho}{dt}=
 \Gamma\left(\hat{c}\hat\rho\hat{c}^\dagger-\mathcal{H}\hat{c}^\dagger\hat{c}\hat\rho\right)
+\Gamma e^{-\beta\omega_o}\left(\hat{c}^\dagger\hat\rho\hat{c}-\mathcal{H}\hat{c}\cop^\dagger\hat\rho\right)
\\
~~~ - ~\epsilon_d[\hat{c}^\dagger-\hat{c},\hat\rho]
\end{array}
\end{eqnarray}
where $\Gamma$ and $\beta$ are determined by
\begin{equation}
\label{eq:rate_d}
\Gamma=\Gamma_e+\frac{\gamma_m}{4}+\frac{\kappa_f^2}{\gamma_m}+\kappa_f,
\end{equation}
\begin{equation}
\label{eq:rate_u}
e^{-\beta\omega_o}\Gamma=\Gamma_e e^{-\beta_e^-\omega_o}+\frac{\gamma_m}{4}+\frac{\kappa_f^2}{\gamma_m}-\kappa_f.
\end{equation}
\par
We turn the effective temperature $\beta$ positive by choosing the feedback strength above the following threshold:
\begin{equation}
\kappa_f > -\kappa_e^-=\frac{1}{2}\Gamma_e\left(e^{-\beta_e^-\omega_o}-1\right).
\end{equation}
Note that the driving on the rhs of the master equation \eqref{eq:me_final} can be absorbed into the standard
thermal dissipator at (inverse) temperature $\beta$ if we displace $\hat{c},\hat{c}^\dagger$ by a suitable real number.
Accordingly, the master equation  \eqref{eq:me_final} must have a unique stationary state which is the following displaced
thermal state of the HO:
\begin{equation}
\label{eq:final_state}
\hat\rho^{\infty} =\mathcal{N}\exp\left[-\beta\omega_o\left(\hat{c}-c_\infty\right)^{\dagger}\left(\hat{c}-c_\infty\right)\right]
\end{equation}
with the static real displacement in interaction picture:
\begin{equation}
\label{eq:c-const-displacement}
c_\infty=-\frac{\epsilon_d}{\kappa_f+\kappa_e^-}<0.
\end{equation}
In the Schr\"{o}dinger picture  the stationary state is a thermal state with the rotating displacement:
\begin{equation}
\label{eq:final_state-Sch}
\hat\rho^{\infty} \Rightarrow \mathcal{N}\exp\left[-\beta\omega_o\left(\hat{c}-c_\infty e^{-i\omega_o t}\right)^{\dagger}\left(\hat{c}-c_\infty e^{-i\omega_o t}\right)\right].
\end{equation}
Hence the mean amplitude rotates, its phase is shifted by $-\pi/2$ with respect to the external driving:
\begin{equation}
\label{eq:mean_ampl}
\ave{\hat{c}}_\infty=c_\infty e^{-i\omega_o t}.
\end{equation}
The average population is the Planckian thermal value plus the yield of displacement:
\begin{equation}
\ave{\hat{c}^{\dagger}\hat{c}}_\infty =  \frac{1}{e^{\beta\omega_o}-1} +\vert c_\infty\vert^2 \equiv  n_o +\vert c_\infty\vert^2 ,
\end{equation}
we use the redundant expression   $\vert c_\infty\vert^2$ for $c_\infty^2$ to capture an occasionally
different phase convention of driving. Both terms on the rhs diverge at the edge of the regime of operation
$\kappa_f+\kappa_e^-\rightarrow+0$ where the model breaks down because it violates the weak coupling condition \eqref{eq:weak}.
\section*{Appendix E: Power Extraction Via Classical periodic Field }
\label{sec:e}
Steady state power extraction \emph{without} storing work is possible by just by driving the engine directly without the flywheel.
Power is gained by  amplification of a classical rotating field in resonance with the two TLS's.
The interaction Hamiltonian is given by  $\hat K(t)=-i\epsilon \left(\hat{a} \hat{b}^{\dagger}e^{i\omega_o t} -\hat{a}^{\dagger} \hat{b} e^{-i\omega_o t}\right)$.
For weak driving the master equation for the two TLS's, $\hat\rho_{hc}$, in the interaction picture of $\hat H_h$ and $\hat H_c$ is,
\begin{equation}
\label{eq:mse_classical_field }
\frac{d\hat\rho_{hc}}{dt}=-\left[\epsilon\left(\hat{a}\hat{b}^{\dagger}-\hat{a}^{\dagger}\hat{b} \right),\hat\rho_{hc}  \right]+\mathcal{L}_h\hat\rho_{hc}+\mathcal{L}_c\hat\rho_{hc},
\end{equation}
where $\mathcal{L}_{h(c)}$ are defined in (\ref{eq:LGKS}).
The master equation (\ref{eq:mse_classical_field }) posses a unique stationary state.
The stationary output power,
\begin{widetext}
\begin{equation}
-\mathcal{P}^{\infty}=\frac{4\epsilon^2 \omega_o (n_h-n_c)}{4\epsilon^2\left( \Gamma_h^{-1}(1-n_h)+\Gamma_c^{-1}(1-n_c)\right)+\Gamma_h\left( 1+e^{-\beta_h\omega_h}\right)+\Gamma_c\left( 1+e^{-\beta_c\omega_c}\right)  }>0,
\end{equation}
\end{widetext}
is positive.
This implies that steady state power extraction can be obtained from a periodically driven field.
Note that for strong driving there is also a steady state power extraction from the engine.
Nevertheless, the master equation (\ref{eq:mse_classical_field }) must be modified.
Derivation of a master equation  driven by a strong periodic field for such a model can be found in \cite{levy212}.