\section{Introduction}
Consider two systems, C for `cold' and H `hot', in thermal states at temperatures $T_{\rm C} < T_{\rm H}$. If the overall system is isolated, energy is conserved, and there are no initial correlations between C and H, heat will flow on average  from the hot to the cold body:
\begin{equation}
\label{eq:uncorrelatedbound}
Q \equiv Q({\rm H} \rightarrow {\rm C}) \leq 0,
\end{equation}
as mandated by the second law of thermodynamics. Eq.~\eqref{eq:uncorrelatedbound} can be derived for both classical and quantum systems  and understood as an average over heat fluctuations that satisfy Jarzynski's exchange fluctuation theorem (XFT)~\cite{jarzynski2004classical}.

The situation is less straightforward if the initial state is locally thermal but correlated, since correlations allow for temporary  `backflows', i.e. $Q > 0$. \cite{lloyd1989use, jennings2010entanglement}.  However, as long as the system is simply \emph{classically} correlated, we can still understand these flows as emerging from averages of underlying classical fluctuations~\cite{jennings2012exchange}.

This picture breaks down in the presence of quantum correlations. In quantum theory, access to the energy of the system, necessary for defining fluctuations classically, can be acquired only via a measurement process that destroys the very quantum properties we wish to investigate. A striking example is given by \emph{strong heat backflows}, i.e. the observation of an amount of heat flowing from C to H larger than $(\Delta \beta)^{-1} \log d$, where \mbox{$\Delta \beta = \beta_{\rm C} - \beta_{\rm H}$, $\beta_{\rm X} = 1/(kT_{\rm X})$}, with $k$ Boltzmann's constant and $d$ the local dimension. Observation of such flows implies that the underlying state is entangled~\cite{jennings2010entanglement}.
If we attempt to understand this flow as originating from  underlying (classical) heat fluctuations, the energy measurements required to access them destroy the entanglement and so the effect disappears. In fact, any stochastic interpretation of heat flows that reproduces the observed $Q$ faces severe no-go theorems~\cite{perarnau2017no, lostaglio2018quantum}. This is the reason why the leading proposal for defining and measuring heat and work fluctuations, the `two-projective-measurement' (TPM), clearly cannot reproduce quantum effects such as \textit{strong} heat  backflows.

 We propose that \emph{quantum} heat fluctuations should be associated with a \emph{quasiprobability} distribution. This arises naturally from a correction to the TPM distribution and can be accessed by weak measurements. We show that the associated negative `probabilities' have a clear thermodynamic interpretation as contributions to the heat flows which cannot be explained within the framework of stochastic thermodynamics. These quasiprobabilities satisfy an XFT that incorporates quantum correlations in the initial state and includes in the relevant classical limits the original XFT of Ref.~\cite{jarzynski2004classical} as well as its extension to classically correlated systems of Ref.~\cite{jennings2012exchange}.  More importantly, our results provide heat flow inequalities that can only be violated in the presence of negativities. These violations are a proxy for a strong form of nonclassicality known as contextuality and remarkably can be evaluated by means of TPM alone.

The challenges posed by the construction of a theory of heat fluctuations reproducing the observed  flows are a strong indication that certain flow configurations are inherently quantum mechanical. We propose that quasiprobabilities, a tool already successfully applied in fields such as quantum optics and quantum computing, can be important diagnostic tools for exploring quantum effects in thermodynamics as our ability to measure heat flows improves.

In Sec.~\ref{sec:problem} we present the problem of capturing quantum phenomena in fluctuation theorems and introduce the concept of quasiprobabilities as the solution.
In Sec.~\ref{sec:quasi}
we introduce a quasiprobability distribution for heat fluctuation and discuss its thermodynamic interpretation, its quantum nature, and the classical limit.
Sec.~\ref{sec:qubits} focuses on the two qubit system, and presents heat flow inequalities that witness negativity. We further demonstrate the result on recent experimental data.
In Sec.~\ref{sec:qxft} we present the full quantum heat exchange fluctuation theorem (QXFT), and extend the results on heat flow inequalities to arbitrary dimension systems. In Sec.~\ref{sec:highD} we show that witnessing negativity in higher dimensional systems is possible using  projective energy measurement alone, and we demonstrate this on a two qutrit system. Finally in Sec.~\ref{sec:outlook} we summarize and discuss future possible inquires.

\subsection*{Comparison with previous works}

 In recent years several proposals have attempted to go beyond some of the limitations faced by TPM schemes. In \cite{deffner2016quantum} a modified quantum Jarzynski equality is derived, which accounts for the thermodynamic cost of the measurement. This in turn led to approaches based on an (initial) one-time measurement scheme for closed \cite{beyer2020work} and open \cite{sone2020quantum} systems, which solve experimental and conceptual difficulties related to performing the second projective measurement in the TPM scheme. Here, on the other hand, we want to	confront the loss of information about quantum effects caused by the first projective measurement of the TPM scheme.

	The definition of fluctuations in the quantum regime is not unique. There is in fact a zoo of such proposals that reduce to the TPM scheme in the case of initial states diagonal in the energy basis, and whose first moments coincide with the average energy change (of the undisturbed system), see review~\cite{baumer2018}.
These proposals necessarily fall into one of two categories~\cite{perarnau2017no, lostaglio2018quantum}:
	\begin{enumerate}
		\item \emph{Lack of positivity:} approaches that give up the idea of assigning probabilities to fluctuations, by allowing `negative probabilities' (Margenau-Hill quasiprobability~ \cite{allahverdyan2014nonequilibrium}, full counting statistics \cite{solinas2015full, solinas2016probing, xu2018effects}, consistent histories \cite{miller2016time, goldstein1995linearly}). These approaches call for an interpretation of such negativities.
		\item \emph{Lack of $\sigma$-additivity:} approaches based on probabilities, which give up the linearity in the density operator of the initial state (Hamilton-Jacobi framework \cite{sampaio2017impossible}, quantum
		Bayesian network \cite{micadei2020quantum}, one-time final energy measurement \cite{gherardini2020role}). For example, if we consider two preparations $\rho_{\rm Heads}$, $\rho_{\rm Tails}$, in these approaches the work probability for  $\frac{1}{2}\rho_{\rm Heads} + \frac{1}{2}\rho_{\rm Tails}$ (tossing a fair coin and preparing $\rho_{\rm Heads}$ or  $\rho_{\rm Tails}$ accordingly) is \emph{not} the equal mixture of the work probability for  $\rho_{\rm Heads}$ and $\rho_{\rm Tails}$. These approaches require an interpretation for the lack of this property, which is normally a consequence of the standard probabilistic interpretation.
	\end{enumerate}
	Independently of the preferred alternative, given the proliferation of  definitions one needs to set  clear objectives for these investigations. It is not complex to come up with novel definitions of quantum fluctuations, leading to  suitably complicated fluctuation theorems which differ from the standard ones due to the presence of quantum coherence. To make a proposal relevant, however, one should attempt to solve three main problems:
	\begin{enumerate}
		\item[(a)] Give an answer to the above-mentioned conceptual questions, providing a thermodynamic interpretation of the failure of properties expected in the statistical interpretation and showing how these are recovered in the classical limit.
		\item[(b)] Find a definition of fluctuations leading to a notion of `quantum signature' that can be identified with a clear-cut notion of nonclassicality. In other words, one should define a precise notion of classicality and show that the identified quantum signature is an \emph{operational} feature that cannot be reproduced by \emph{any} model satisfying that classicality property.
		\item[(c)] Show that such quantum signatures are \emph{experimentally accessible} without requiring a very large amount of control (like knowledge of the initial state or the dynamics). Ideally witnessing quantum signatures should necessitate no more control than that required by the original TPM scheme.
	\end{enumerate}
What we set out to do in this work is to show that all these problems can be solved simultaneously in the context of heat fluctuations.

\section{The problem of quantum heat fluctuations}
\label{sec:problem}
 The impasse in the definition of heat and work fluctuations can be traced back to the information-disturbance tradeoff of quantum mechanics. Imagine the cold and hot systems, C and H, initially described by a density operator $\rho_{{\rm C} {\rm H} }$, evolve under a closed dynamics described by a unitary $U$, $\rho_{{\rm C} {\rm H} }(\tau) = U \rho_{{\rm C} {\rm H}} U^\dag$. The average heat flow reads
\begin{equation}
\label{eq:weak}
Q := \tr{}{\rho_{{\rm C} {\rm H}} H_C} - \tr{}{\rho_{{\rm C} {\rm H} }(\tau) H_C}.
\end{equation}
In order to study fluctuations, classically we want to see this quantity as emerging from a microscopic process where energy $q$ is exchanged with probability $p(q)$ such that $Q=~\int dq p(q) q$.  In this case, the system  C starts at some phase space point $z_0$ with probability $p(z_0)$ and evolves deterministically to $z_\tau$, in which case $q = E(z_0) - E(z_\tau)$, with $E(z)$ the energy at point $z$. Hence
\begin{equation}
p(q) = \int dz_0 p(z_0) \delta(q - (E(z_0) - E(z_\tau)))
\end{equation}
 Quantumly, however,  the initial energy  is generally not well-defined. Non-commutativity forbids the construction of a joint probability (linear in $\rho_{{\rm C} {\rm H} }$) that correctly reproduces the energy distributions of the initial and final states \cite{perarnau2017no}.

This long-standing problem dates back many decades and the standard way to overcome it has been to define fluctuations by a two-point-measurement (TPM) scheme \cite{campisi2011colloquium, esposito2009nonequilibrium}. The TPM scheme involves (1) projectively measuring the energy of C at the start, (2) letting the system dynamics unfold, and (3) projectively measuring the energy of C again at the end.  $p^{\rm TPM}(q)$ is then defined as the probability of recording an energy difference between the outcomes of the two projective measurements equal to $q$.
 However, in this scheme $Q^{
\rm TPM} := \int dq  p^{\rm TPM}(q) q \neq Q$, since effects due to non-commutativity are expunged by the initial invasive measurement. The widely used TPM method has the advantage that, by construction, classical fluctuation theorems hold for $p^{\rm TPM}(q)$. The drawback is that we are effectively describing the heat fluctuations of a state disturbed (decohered) by the initial energy measurement, rather than those of the original state $\rho_{{\rm C} {\rm H} }$. Quantum effects such as strong heat backflows provably never appear. Hence this solution cannot be satisfactory in the fully quantum regime, when $\rho_{{\rm C} {\rm H}}$ includes energetic coherence and quantum correlations.

\subsection*{A different solution: quasiprobabilities} We argue that the root cause of the present difficulties is the attempt to give a classical stochastic account of quantum heat fluctuations.
Since quantum mechanics does not allow for the construction of joint probabilities  of noncommuting quantities, a natural framework is to describe them by means of quasiprobabilities~\cite{hofer2017quasi}, i.e. real-valued functions $p(q)$ that sum to $1$ but are allowed to take values outside $[0,1]$. While relatively uncommon in the thermodynamic setting (but see e.g. \cite{allahverdyan2014nonequilibrium, solinas2015full, solinas2016, xu2018effects, lostaglio2018quantum}), the use of quasiprobabilities to probe quantum effects is widespread in other fields of quantum sciences, including quantum transport \cite{belzig2001full,bednorz2010quasiprobabilistic,hovhannisyan2019quantum},  quantum optics~\cite{bucca, wigner1932quantum}, quantum computing \cite{pashayan2015estimating,seddon2020quantifying} and condensed matter \cite{halpern2017jarzynski, alonso2019out}.
 Quasiprobabilities originated with the attempt to give a phase space description to quantum mechanics despite non-commutativity of position and momentum. The most famous example is the Wigner function \cite{wigner1932quantum}, which associates a quantum state with a quasiprobability distribution over positions and momenta whose marginals correctly reproduce the statistics collected in the corresponding quantum measurements \footnote{Interestingly, Wigner originally introduced his function to compute quantum corrections to thermodynamic quantities.}.

The analogy to our scenario is that, in order to recover Eq.~\eqref{eq:weak}, we wish to reproduce the correct quantum statistics of energy measurements performed on the initial $\rho_{{\rm C} {\rm H} }$ and final $\rho_{{\rm C} {\rm H} }(\tau)$ states.
Formally, if $p^{\rm W}_{i_\Co i_\Ho \veryshortarrow  f_\Co f_\Ho}$ denotes the quasiprobability of starting with energies $(E^{\rm C}_{i_{\rm C}}, E^{\rm H}_{i_{\rm H}})$ before the process $U$ and ending up in the final energies  $(E^{\rm C}_{f_{\rm C}}, E^{\rm H}_{f_{\rm H}})$ afterwards, then we wish
\begin{equation}
\nonumber
\sum_{f_{\rm C}, f_{\rm H}} p^{\rm W}_{i_\Co i_\Ho \veryshortarrow  f_\Co f_\Ho} = \tr{}{\rho_{{\rm C} {\rm H}} \Pi^{i_{{\rm C}} i_{{\rm H}}}}:=  p^{i_{{\rm C}} i_{{\rm H}}},
\end{equation}
 \begin{equation}
\label{eq:marginalproperty}
\sum_{i_{\rm C}, i_{\rm H}} p^{\rm W}_{i_\Co i_\Ho \veryshortarrow  f_\Co f_\Ho} = \tr{}{\rho_{{\rm C} {\rm H}}(\tau) \Pi^{f_{{\rm C}} f_{{\rm H}}}}:=  p^{f_{{\rm C}} f_{{\rm H}}}(\tau),
\end{equation}
where \mbox{$\Pi^{i_{{\rm C}} i_{{\rm H}}} \equiv \Pi^{i_{{\rm C}}} \otimes \Pi^{i_{{\rm H}}}$} and $\Pi^{i_{\rm X}}$ are energy projection operators such that $H_{\rm X} = \sum_{i_{\rm X}} E^{\rm X}_{i_{\rm X}} \Pi^{i_{\rm X}}$ (X = C, H). Also note that $p^{i_{{\rm C}} i_{{\rm H}}}$ ($p^{i_{{\rm C}} i_{{\rm H}}}(\tau)$) denotes the statistics of an energy measurements carried out at the beginning (end) of the experiment.
The condition in Eq.~\eqref{eq:marginalproperty} ensures that the fluctuations  $p^{\rm W}_{i_\Co i_\Ho \veryshortarrow  f_\Co f_\Ho}$ reproduce the average heat flows of Eq.~\eqref{eq:weak}:
\begin{equation}
\label{QMargenau} \sum_{i_{\rm C}, i_{\rm H}, f_{\rm C}, f_{\rm H}} p^{\rm W}_{i_\Co i_\Ho \veryshortarrow  f_\Co f_\Ho} (E_{i_{\rm C}} - E_{f_{{\rm C}}}) = Q.
\end{equation}
 \section{A quasiprobability for quantum heat fluctuations}
\label{sec:quasi}
 How do we choose a quasiprobability distribution among the many possible candidates? The analogy with quantum optics suggests different choices will be relevant in different scenarios. Nevertheless, one candidate stands out.  Recall that the TPM distribution reads
 	\begin{equation}
 	p^{\rm TPM}_{i_\Co i_\Ho \veryshortarrow  f_\Co f_\Ho} := \tr{}{ \Pi^{f_{\rm C} f_{\rm H}}(\tau)   \Pi^{i_{\rm C} i_{\rm H}} \mathcal{D} \rho_{{\rm C} {\rm H}}},
 	\end{equation}
 	where $ \Pi^{f_{\rm C} f_{\rm H}}(\tau) :=  U^\dag  \Pi^{f_{\rm C} f_{\rm H}} U$ and $\mathcal{D}$ denotes dephasing, the operation that removes all off-diagonal elements in the energy basis. Moreover, the TPM characteristic function is associated with a two-time correlation function \emph{computed on a decohered version} $\mathcal{D}\rho_{{\rm C} {\rm H}}$ \emph{of the initial state} $\rho_{{\rm C} {\rm H}}$ \cite{campisi2011colloquium}. This naturally suggests reintroducing coherence and quantum correlations by removing such dephasing:
 \begin{equation}
\label{eq:weakquasiprobability}
p^{\rm W}_{i_\Co i_\Ho \veryshortarrow  f_\Co f_\Ho} := {\rm Re} \, \tr{}{ \Pi^{f_{\rm C} f_{\rm H}}(\tau)   \Pi^{i_{\rm C} i_{\rm H}} \rho_{{\rm C} {\rm H}}},
\end{equation}
where we focus here on the real part (the imaginary part may be of independent interest \cite{kunjwal2018anomalous}). This is our proposal for a quasiprobability distribution for heat fluctuations in the quantum regime.

Since every $\rho_{{\rm C} {\rm H}}$ can be additively decomposed as  a term $\mathcal{D} \rho_{{\rm C} {\rm H}}$ (block-)diagonal in the energy basis plus energetic coherences $\chi_{{\rm C} {\rm H}}$,  i.e. $\rho_{{\rm C} {\rm H}} = \mathcal{D} \rho_{{\rm C} {\rm H}} + \chi_{{\rm C} {\rm H}}$,  corrections to the standard scheme emerge only if \mbox{$\chi_{{\rm C} {\rm H}} \neq 0$}:
\begin{equation*}
p^{\rm W}_{i_\Co i_\Ho \veryshortarrow  f_\Co f_\Ho} = p^{\rm TPM}_{i_\Co i_\Ho \veryshortarrow  f_\Co f_\Ho} + {\rm Re} \, \tr{}{ \Pi^{f_{\rm C} f_{\rm H}}(\tau)   \Pi^{i_{\rm C} i_{\rm H}} \chi_{{\rm C} {\rm H}}}.
\end{equation*}
The quantum correction to $p^{\rm TPM}_{i_\Co i_\Ho \veryshortarrow  f_\Co f_\Ho}$ is bounded by $\| \chi_{{\rm C} {\rm H}}\|$, which is a coherence measure \cite{aberg2006superposition, marvian2016how}.
 $p^{\rm W}$ satisfies the marginal properties of Eq.~\eqref{eq:marginalproperty} and hence it reproduces the quantum heat flow of Eq.~\eqref{QMargenau}. As expected,  when $\chi_{{\rm C} {\rm H}} \neq 0$ in general $p^{\rm W}$ does not describe a stochastic process, since it can be negative: in fact $p^{\rm W}_{i_\Co i_\Ho \veryshortarrow  f_\Co f_\Ho}~\in~[-1/8,1]$\footnote{Ref.~\cite{allahverdyan2014nonequilibrium} (Eq. 41) or  \cite{zhang2011matrix} (Theor.~7.5).}.  It is worth noting that $p^{\rm W}$ recovers the definition of a quasiprobability distribution for two noncommuting observables
given by Margenau and Hill in 1961~\cite{margenau1961correlation} and known since the late 1930s \cite{terletsky1937limiting}.
It also takes the same form as the work quasiprobability proposed by Allahverdyan in Ref.~\cite{allahverdyan2014nonequilibrium},  whose negativities were however dismissed as an unphysical feature. Here we present a different picture.

Firstly, while $p^{\rm W}$ can show negativities, it can be experimentally accessed in a probing scheme that naturally extends the TPM scheme
(see App.~\ref{app:A} for more details). To estimate $p^{\rm W}$ one has to replace the first  of the two projective measurements of the TPM scheme with a \emph{weak measurement},  which is achieved by very weakly coupling the system to either a continuous variable pointer or a qubit probe  which is then projectively measured. Appealingly, $p^{\rm W}$ can hence be reconstructed in a `minimally disturbing' version of the TPM scheme.

 Secondly, negativity is a precise signature of nonclassicality. In fact, $p^{\rm W}_{i_\Co i_\Ho \veryshortarrow  f_\Co f_\Ho}$ is proportional to a quantity inferred from weak measurements and known as the (generalized) \emph{weak value}~\cite{dressel2015weak}, first introduced by Aharonov \emph{et.al.} \cite{aharonov1988result,dressel2014colloquium} and later generalized to mixed states \cite{wiseman2002weak, dressel2010contextual}. It then follows from the extensions \cite{lostaglio2018quantum, kunjwal2018anomalous} of Pusey's theorem \cite{pusey2014anomalous} that the negativities of $p^{\rm W}$ are witnesses of a strong form of nonclassicality known as generalised contextuality \cite{spekkens2005contextuality}. Specifically, the statistics collected by the weak measurement scheme probing $p^{\rm W}$ cannot be explained by \emph{any} noncontextual hidden variable model, a claim valid even in the presence of noise~\cite{kunjwal2018anomalous}.  This is directly related to the absence of \emph{any} classical stochastic process explaining the relevant statistics~\cite{spekkens2008negativity}.

Thus far we showed that $p^{\rm W}$ is a natural extension of the TPM scheme, both formally and experimentally; that it reproduces the heat flows of $Q$, including strong heat backflows due to initial entanglement; and that its negativity is a precise signature of strong nonclassicality.
However, what is the \emph{thermodynamic} significance of negativity and how does it impact the observed heat flows? We now answer these questions.

\subsection{Thermodynamic role of negativities in the heat fluctuations}
 To single out the meaning of negativities, let us split the quasiprobability in positive and negative components: $ p^{\rm W}_{i_\Co i_\Ho \veryshortarrow  f_\Co f_\Ho} = p^+_{i_{\rm C} i_{\rm H} \rightarrow f_{\rm C} f_{\rm H}} + p^-_{i_{\rm C} i_{\rm H} \rightarrow f_{\rm C} f_{\rm H}}$, where $p^\pm_{i_{\rm C} i_{\rm H} \rightarrow f_{\rm C} f_{\rm H}}=  p^{\rm W}_{i_\Co i_\Ho \veryshortarrow  f_\Co f_\Ho}$ if $ p^{\rm W}_{i_\Co i_\Ho \veryshortarrow  f_\Co f_\Ho}$ is positive (negative), and zero otherwise.
Recall that, with our sign convention, $Q>0$ means backflow. Then, the net flow $Q$ can be decomposed into direct (H $\rightarrow$ C) and back (C~$\rightarrow$~ H) flows:
\begin{equation}
Q = Q^{\rm back} - Q^{\rm direct},
\end{equation}
 where the `back' term includes all contributions to the flow ${\rm C} \rightarrow {\rm H}$ and the `direct' term includes all contributions ${\rm H} \rightarrow {\rm C}$:
\begin{align}
Q^{\rm back} = \sum_{E_{i_{\rm C}}>E_{f_{{\rm C}}}} (p^+_{i_{\rm C} i_{\rm H} \veryshortarrow f_{\rm C} f_{\rm H}} - p^-_{ f_{\rm C} f_{\rm H} \veryshortarrow i_{\rm C} i_{\rm H} }) \Delta E_{i_{\rm C} f_{{\rm C}}},
\label{Qback} \\
Q^{\rm direct} = \sum_{E_{i_{\rm C}}>E_{f_{{\rm C}}}} (p^+_{f_{\rm C} f_{\rm H} \veryshortarrow i_{\rm C} i_{\rm H} } - p^-_{ i_{\rm C} i_{\rm H} \veryshortarrow f_{\rm C} f_{\rm H} })  \Delta E_{i_{\rm C} f_{{\rm C}}}.
\label{Qdirect}
\end{align}
These equations give a very suggestive interpretation of the role of negative probabilities in heat flows (see Fig.~\ref{fig:heatflows}). Since $E_{i_{\rm C}} > E_{f_{{\rm C}}}$, $Q^{\rm back}$ has two positive contributions:
1.~From the transition $i_{{\rm C}} \rightarrow f_{{\rm C}}$, which removes energy from C (as expected).
2.~From the transition $f_{{\rm C}} \rightarrow i_{{\rm C}}$ \emph{when}  $p^-_{ f_{\rm C} f_{\rm H} \veryshortarrow i_{\rm C} i_{\rm H} }<0$.  Classically, the transition $f_{{\rm C}} \rightarrow i_{{\rm C}}$ would add energy to C, yet it contributes to the backflow from C to H when the correspondent quasiprobability turns negative.
Symmetrically, a transition taking energy away from C can contribute to the direct-flow, when $p^-_{ i_{\rm C} i_{\rm H} \rightarrow f_{\rm C} f_{\rm H} }<0$. These counterintuitive heat flow contributions signal the breakdown of a stochastic description and the onset of quantum effects,  as depicted in Fig.~\ref{fig:heatflows}.

\begin{figure*
\centering
\subfloat{\includegraphics[width=0.48\linewidth]{heatflows13a.pdf}}
\qquad
\subfloat{\includegraphics[width=0.48\linewidth]{heatflows13b.pdf}}
	\caption{ Schematic illustration of the classical and quantum contributions to the overall heat flow $Q$ between a cold (C) and a hot (H) system. $Q$ is an average of direct (H $\rightarrow$ C) and back flows (C $\rightarrow$ H). {\bf a}, Classical contributions (positive probability events) corresponds to standard contributions, where the loss of energy from H contributes to the flow towards C and vice versa. {\bf b}, Quantum contributions (negative probability events) effectively act as contributions in which the loss of energy from H contributes to the flow (C $\rightarrow$ H).}
	\label{fig:heatflows}
\end{figure*}

 In particular, we can understand strong flows as the result of negative probabilities, which are destroyed by the measurement stage of the TPM scheme  (a projective energy measurement on C and H forces $p^{\rm W}>0$). Furthermore, Eq.~\eqref{Qdirect} suggests  \textit{negativities may lead to strong direct-flows as well}.

 It is important to bear in mind that negativity will not always result in strong heat flows.
Since only $Q =Q^{\rm back} - Q^{\rm direct}$ is observed, negativites from $Q^{\rm back}$ and $Q^{\rm direct}$ must not cancel each other if we are going to observe an overall effect in the heat flow.

\subsection{Classical limit}

 When discussing quantum effects in heat or work fluctuations, it is important to distinguish between two different forms:
	\begin{enumerate}
		\item Quantumness due to  correlations or entanglement in the initial state, such as that responsible for strong heat backflows.
		\item Quantumness due to the superpositions created by $U$ \emph{during} the dynamics.
	\end{enumerate}

The first form of nonclassicality is the one preventing a straighforward description of quantum heat exchanges by means of a classical stochastic process. Classicality in the first sense  emerges in the expected regimes. The typical mechanism is decoherence. If ${\rm C}$ and ${\rm H}$ are not protected from the external environment, superpositions of different energy states are quickly suppressed ($\chi_{{\rm C} {\rm H}} \rightarrow 0$) and so $p^{\rm W}~\rightarrow ~p^{\rm TPM} \geq 0$, i.e. negativity disappears.
 Decoherence may be also due to the measurement process, which is the mechanism that introduces classicality in the TPM scheme.

  Another mechanism relevant at mesoscopic scales is coarse-graining. Suppose one can only access the heat exchanges between two clusters, each containing $N$ particles, but cannot keep track of the process at the level of the individual constituents. We idealise the situation via an average Hamiltonian $H_{{\rm C} {\rm H}} \approx \sum_{i=1}^N (H_{{\rm C}_i} + H_{{\rm H}_i})/N$ and a dynamics $U \approx \otimes_{i=1}^N U_{{\rm C}_i {\rm H}_i}$.
	 Then, whatever the initial state $\rho_{{\rm C} {\rm H} }$ is, $p^{\rm W} \approx p^{\rm TPM} + O(1/N)$ \footnote{This follows from the fact that $[U^\dag H_{{\rm C} {\rm H}}U, H_{{\rm C} {\rm H}}] \propto 1/N$.}.
	 The ratio between fluctuations and average flows, on the other hand, is suppressed as $O(1/\sqrt{N})$. There is hence an intermediate, ``mesoscopic'' regime in which quantum effects are irrelevant but fluctuations are not yet negligible.

	  The second form of nonclassicality was studied in \cite{lloyd2015no}. In this study it was shown that for initially uncorrelated bipartite system energy cannot flow without generation of quantum correlations in the form of discord during the unitary process. This form of quantumness is the one responsible for the differences between the classical stochastic process predicted by the TPM scheme and the one obtained within a fully classical description (e.g., a classical mechanics model), whenever the latter is well-defined. Note that, since we recover the TPM probability in the regimes discussed above, we can exploit previous results on the emergence of this second form of classicality when the superpositions created by $U$ \emph{during} the dynamics can be neglected~\cite{jarzynski2015quantum}.

	  We conclude that our definition of heat fluctuations appears to provide a clear path to the emergence of classicality, from negative quasiprobability all the way up to a classical mechanical model, passing through the standard TPM scheme that plays the role of a ``quasi-classical'' description.

\section{Heat flows witnessing negativities in a two qubit system}
\label{sec:qubits}
 By attempting to derive a  stochastic interpretation of quantum fluctuations, we arrived at a picture in which quantum effects are incorporated as non-classical features in the heat exchange, such as contributions to the energy flow from C to H associated with transitions to \emph{higher} levels in C. From a practical perspective, however, the next step is to ask if these non-classical features can be directly related to observable effects and if these effects can be accessed in realistic experiments. We answer both questions in the positive.

 We start with the archetypal (and fully solvable)  idealized two-qubit scenario.
Since $U$ injects no work, nontrivial energy exchanges can only happen if the two qubits are resonant.
\begin{theorem} \label{th:qubitspecificbounds}
Let $\rho_{{\rm C} {\rm H}}$ be a two qubit system with thermal marginals ($\beta_{{\rm C}} \neq \beta_{{\rm H}}$), $H_{\rm C} = H_{\rm H} = E \ketbra{1}{1}$ and $U$ a unitary injecting no work, $[U, H_{\rm C} + H_{\rm H}] = 0$. If $p^{\rm W}$ is nonnegative,
	\begin{equation}
	\label{eq:qubitbound}
	|Q| \leq \frac{2 + e^{\beta_{\rm H} E} + e^{\beta_{\rm C} E} }{e^{\beta_{{\rm C}}E} - e^{\beta_{{\rm H}}E}} |Q^{\rm TPM}|.
	\end{equation}
	\end{theorem}
See App.~\ref{subsec:proof1} for the proof.  This `heat flow inequality' shows that strong enough direct or backflows of heat can only occur when negativity is at play in $p^{\rm W}$. Note that to violate the bound of Eq.~\eqref{eq:qubitbound}, we necessarily  witness $|Q| > |Q^{\rm TPM}|$, i.e. a direct or inverse flow bigger than the corresponding flow observed in the TPM scheme. While violations of Eq.~\eqref{eq:uncorrelatedbound},  reported in Ref.~\cite{micadei2019reversing}, can occur in a purely classical setup, due to initial correlations in the energetic degrees of freedom, violations of Inequality~\ref{th:qubitspecificbounds} imply negativity, which is a proxy for genuinely quantum effects.

\subsection{Experimental verification}

Next, we analyze the violations of Inequality \ref{th:qubitspecificbounds}. In  recent NMR experimental setups~\cite{micadei2019reversing,pal2018experimental}, the heat backflow between two local-thermal qubits with initial correlations was measured. Here we focus on the data of Ref.~\cite{micadei2019reversing} and
 show that, while other general approaches (e.g., Ref.~\cite{jennings2010entanglement}) fail to witness quantum effects for this setup,
Inequality~\ref{th:qubitspecificbounds} is violated.

 In the experiment, the heat backflow $Q$ between two nuclear spin $1/2$, prepared in an initial locally thermal state, was measured. We denote by $\sigma_i$, $i=x,y,z$, the Pauli matrices. The initial Hamiltonian is $H_{{\rm H} {\rm C}} = H_{{\rm H}} + H_{{\rm C}}$ with $H_{{\rm C}({\rm H})} = h \nu (1- \sigma^{{\rm C}({\rm H})}_z)$  and $h$ Planck's constant.  A locally thermal state is prepared with the form
\begin{equation}
\label{eq:stateexperiment}
\rho_{{\rm H} {\rm C}}=\rho_{{\rm H}}\otimes\rho_{{\rm C}}+\gamma\ketbra{01}{10}+\gamma^{*}\ketbra{10}{01},
\end{equation}
where $\rho_{{\rm H}({\rm C})} = e^{-\beta_{{\rm H}({\rm C})} H_{{\rm H}({\rm C})}}/ \tr{}{e^{-\beta_{{\rm H}({\rm C})} H_{{\rm H}({\rm C})}}} $. The unitless rescaled inverse temperatures  reported in the experiment are  $\beta_{{\rm H}}=0.9618$ and $\beta_{{\rm C}}=1.13$, with the gap $\nu=1{\rm kHz}$. In the correlated scenario, $\gamma=-0.19$, whereas in the uncorrelated scenario  $\gamma=0$.

During the process, the nuclear spins are coupled via the interaction Hamiltonian \mbox{$H_{\rm int}=J\pi \hbar/2(\sigma_x^{{\rm H}}\sigma_y^{{\rm C}}-\sigma_y^{{\rm H}}\sigma_x^{{\rm C}})$} with $J=215.1{\rm Hz}$. Since $[H_{\rm int}, H_{{\rm H} {\rm C}}] = 0$, this generates an energy-preserving unitary. A simple yet important observation is the following: the measurement of the heat exchange when $\gamma = 0$ corresponds to the measurement of $Q^{\rm TPM}$ for the case $\gamma = - 0.19$. This follows from the fact that the initial energy measurement on ${\rm C}$ and ${\rm H}$ (the first step of the TPM protocol) has the effect of setting $\gamma= 0$. Thanks to this observation, the data collected in the experiment suffices to test Inequality~\ref{th:qubitspecificbounds}.

In Fig.~\ref{fig:qubitslutz} we plot the heat transfer between qubits ${\rm C}$ and ${\rm H}$ for various interaction times. The green line corresponds to the heat flow when the qubits are initially quantum correlated, whereas the orange line corresponds to initially uncorrelated qubits. Equivalently, this corresponds to the measurements of $Q$ (green) and $Q^{\rm TPM}$ (orange).  $Q>0$ corresponds to backflow from ${\rm C}$ to ${\rm H}$, whereas $Q<0$ corresponds to direct flow from ${\rm H}$ to ${\rm C}$.
\begin{figure}[t!]
	\includegraphics[width=0.85\linewidth]{qubitslutz.pdf}
	\caption{ Heat transfer $Q$ and $Q^{\rm TPM}$  (green and orange curves, respectively) for the experiment in Ref.~\cite{micadei2019reversing}. The shaded areas indicate the regime of negative probabilities. The dark blue shaded area indicates the violation of Inequality~\ref{th:qubitspecificbounds}. }
	\label{fig:qubitslutz}
\end{figure}
The shaded areas indicate the interaction times in which negativity in heat transfer quasiprobability  is present.
The dark blue shaded area indicates the interaction times for which Inequality~\ref{th:qubitspecificbounds} successfully witnesses negativity for the relevant experimental parameters.

While in the experiment only the backflow regime was explored, which is sufficient for detecting negativity, here we theoretically extend the interaction time to the regime where direct flow is observed.  We see that also in this regime, the measurement of the heat flows would allow us to detect the negativity of the quasiprobability,  showing that quantum effects are not only present in backflows.
We further extend the investigation of Inequality~\ref{th:qubitspecificbounds} beyond the specific parameters of Ref.~\cite{micadei2019reversing}, see Fig.~\ref{fig:qubitbounds}.
We conclude that negativity can be detected for a wide range of parameters in accessible experiments.

 Using the experimental parameters of Ref.~\cite{micadei2019reversing}, the smallest eigenvalue of the partial transpose $\rho_{{\rm C} {\rm H}}^{T_{\rm H}}$ of $\rho_{{\rm C} {\rm H}}$ is $\approx 0.0014$, hence $\rho^{ T_{{\rm H}}}_{{\rm C} {\rm H}} \geq 0$. By the Peres-Horodecki criterion, $\rho_{{\rm C} {\rm H}}$ is separable.  This shows that the heat fluctuations can have negativities, as detected by Inequality~\ref{th:qubitspecificbounds}, even when $\rho_{{\rm C} {\rm H}}$ is not entangled.

  We mention in passing a nontrivial result that holds specifically for two qubit systems: \textit{the existence of quantum correlations in the energy degrees of freedom is necessary for witnessing a backflow for an appropriate energy-preserving dynamics} (see App.~\ref{app:B}).  This is the reason the orange curve in Fig.~\ref{fig:qubitslutz} is all below the $x$-axis.

\begin{figure}[t!]
\includegraphics[width=0.8\linewidth]{figqubitsnew2.pdf}
	\caption{
	Region of negative probabilities (large  trianguloid) and violations of Inequality~\ref{th:qubitspecificbounds} (darker trianguloid) in the parameter space $(\theta,\eta)$. $\theta$ characterizes the unitary evolution (the rotation in the $\ket{01}$, $\ket{10}$ subspace) and $\eta$ the initial energetic superposition between $\ket{01}$, $\ket{10}$.
In App.~\ref{app:B} we show that for fixed $\beta_{{\rm C}}$, $\beta_{{\rm H}}$, and a simple reparametrization, the parameter space is restricted to $(\theta, \eta,P_{00})$, with
$P_{00}$ the population in the state $\ket{00}$.
In the experiment \cite{micadei2019reversing},  $\beta_{{\rm C}}=1.13$,  $\beta_{{\rm H}}=0.962$, $P_{00} = 0.547$, and $\eta=-0.19$, which is indicated by the dashed line at the bottom of the figure. }
	\label{fig:qubitbounds}
\end{figure}

\subsection{Nonideal heat exchange}

While it is commonplace to assume that $U$ injects an amount of work into the system that can be neglected in calculations, strictly speaking Inequality~\ref{th:qubitspecificbounds} holds only in an ideal case. Next, we extend Inequality~\ref{th:qubitspecificbounds} to realistic heat exchange scenarios, in which (1) the energy levels are not resonant, (2) the dynamics injects some work:

\begin{theorem} \label{th:qubitspecificboundsnonideal}
	Let $\rho_{{\rm C} {\rm H}}$ be a two qubit system with thermal marginals ($\beta_{{\rm C}} \neq \beta_{{\rm H}}$), $H_{\rm C} = E_C \ketbra{1}{1}$, $H_{\rm H} = E_H \ketbra{1}{1}$ and U a unitary dynamics. Suppose
	\begin{enumerate}
		\item $U$ is close to an energy preserving unitary $\tilde{U}$: \mbox{$
			\| U - \tilde{U} \| \leq \epsilon$}.
		\item The gap is not too wide: $\Delta := \frac{|E_{{\rm C}} - E_{{\rm H}}|}{2 \bar{E} }< \frac{1-R}{1+R}$, where  $\bar{E}:= (E_{{\rm C}} + E_{{\rm H}})/2$ is the average energy and
		\mbox{$R = (1+ e^{\beta_{{\rm H}} E_{{\rm H}}})/ (1+ e^{\beta_{{\rm C}} E_{{\rm C}}})$}.
		\item The energy flow is large enough that we can distinguish a heat contribution: $|Q| > 2 \epsilon \bar{E}$.
	\end{enumerate}
	If $p^{\rm W}$ is nonnegative,
	\begin{small}
		\begin{equation}
		\label{eq:qubitboundnonideal}
		|Q| \leq \frac{1+ R + \Delta (1+R)}{1- R - \Delta (1-R)} |Q^{\rm TPM}| + 4 \epsilon \bar{E} \frac{2 + \Delta (1+R)}{1-R - \Delta(1-R)} .
		\end{equation}
	\end{small}
\end{theorem}
When $\epsilon \rightarrow 0$, $\Delta \rightarrow 0$ we recover Inequality~\ref{th:qubitspecificbounds}. The above inequality is suitable to be applied in realistic experimental scenarios. Note that the condition of unitarity is not crucial and Inequality~\ref{th:qubitspecificboundsnonideal} can be extended to quantum channels. The proof of this inequality, stronger bounds, and a case study are presented in App.~\ref{app:nonidealflow}.  This shows that analogue heat flow inequalities can be derived beyond the ideal regime.

\section{A quantum exchange fluctuation theorem}
\label{sec:qxft}
 The above heat flow inequality is limited to two-qubit systems. We now derive a heat flow inequality valid for arbitrary finite-dimensional systems. This is obtained through the derivation of a quantum exchange fluctuation theorem  (QXFT) in the spirit of \cite{allahverdyan2014nonequilibrium, kwon2019fluctuation}.

Denote by  $\Delta I~:=~ I_{f_{{\rm C}} f_{{\rm H}}}~-~I_{i_{{\rm C}} i_{{\rm H}}}$ with
	\begin{equation}
	I_{i_{{\rm C}} i_{{\rm H}}}=\log \frac{\tr{}{\Pi^{i_{\rm C} i_{\rm H}} \rho_{{\rm C} {\rm H}}} }{\tr{}{\Pi^{i_{\rm C}} \rho_{{\rm C}}} \tr{}{\Pi^{ i_{\rm H}} \rho_{{\rm H}}} },
	\end{equation}
	the elements of the classical mutual information. Assuming $E_{i_{\rm C}} - E_{f_{\rm C}} \approx E_{f_{\rm H}} - E_{i_{\rm H}}$ for all nonzero $p^{\rm W}_{i_\Co i_\Ho \veryshortarrow  f_\Co f_\Ho}$, we obtain the XFT in the full quantum regime
\begin{equation}
\label{eq:XFT}
\ave{e^{\Delta I +\Delta\beta \Delta E_{ i_{{\rm C}} f_{{\rm C}}}}}_{ W}=1+\bar{\chi}.
\end{equation}
The r.h.s of Eq.~(\ref{eq:XFT}) consists of 1, contribution from the populations, and $\bar{\chi}$, involving contributions from the coherence terms in the energy eigenbasis
\begin{equation}
\label{eq:chi}
\bar{\chi}=\sum_{l,k\neq m}\frac{\rho_{ll}}{\rho_{kk}}{\rm Re}\lbrace\rho_{km}\bra{l}U\ket{k}\bra{m}U^{\dagger}\ket{l}\rbrace,
\end{equation}
which are absent in the TPM scheme. Note that
\begin{enumerate}
	\item If $\rho_{{\rm C} {\rm H}}$ is classically correlated in the energy basis, $\bar{\chi} =0$, recovering the main result of Ref.~\cite{jennings2012exchange}.
	\item If, furthermore, $\rho_{{\rm C} {\rm H}}$ is uncorrelated, $\Delta I = 0$ and we recover Jarzynski's original result \cite{jarzynski2004classical}.
\end{enumerate}
Eq.~\eqref{eq:XFT} is thus the generalization of the XFT to the full quantum regime. The proof of the quantum XFT is provided in App.~\ref{app:C}.

\subsection{Witnessing negativity in arbitrary dimension}
Using the quantum XFT, we derive a general inequality witnessing negativity in any dimension: :
\begin{theorem} \label{th:xftbounds}
Let $\rho_{{\rm C} {\rm H}}$ be an arbitrary finite-dimensional system with thermal marginals ($\beta_{{\rm C}} \neq \beta_{{\rm H}}$) and $U$ a unitary that injects no work ($E_{i_{\rm C}} - E_{f_{\rm C}} = E_{f_{\rm H}} - E_{i_{\rm H}}$). If $p^{\rm W}$ is nonnegative,
	\begin{equation}
	\label{eq:xftbound}
	Q \leq -\frac{\ave{\Delta	I}_{\rm W}}{\Delta \beta}+\frac{\log(1+\bar{\chi})}{\Delta \beta}.
	\end{equation}
	\end{theorem}
This result can be derived by applying Jensen's inequality to Eq.~(\ref{eq:XFT}). $\ave{\Delta	I}_{W}$ in Eq.~(\ref{eq:xftbound}) is a classical mutual information-like term satisfying $\ave{\Delta	I}_{\rm W}=\ave{\Delta	I}_{\rm TPM}$ if the initial state only has classical correlations in the energy basis, and  $\ave{\Delta	I}_{W}=0$ in the absence of initial quantum and classical correlations \footnote{Note that the term in Eq.~(\ref{eq:chi}) may diverge as $\rho_{kk}\rightarrow 0$. In  App.~\ref{app:C} we derive a different bound for the heat in Eq.~(\ref{eq:xftbound}) that bypasses this problem. Also note that neither $\ave{\Delta	I}_{\rm W}$ nor $\ave{\Delta	I}_{T\rm PM}$ admit a simple interpretation as the change in classical mutual information in the energetic degrees of freedom.}. Violations of the inequality imply that some quasiprobabilities are necessarily negative. Note that the inequality becomes trivial if no initial coherence is present ($\chi_{{\rm C} {\rm H}} =0$), as expected. In App.~\ref{app:C}  we provide a second heat flow inequality valid for general finite-dimensional systems.

For finite-dimensional systems with equal Hamiltonians and nondegenerate energy gaps, we  prove in App.~\ref{app:C} an additional generic property of negativity:
 \textit{If all quasiprobabilities that contribute to direct or backflow are negative, then the flows are necessarily strong, i.e. $|Q|>|Q^{\rm TPM}|$}.

\subsection{Nonideal heat exchange}
Inequality~\ref{th:xftbounds} can be generalized to situations in which the injected work cannot be neglected. If for all nonzero $p^{\rm W}_{i_\Co i_\Ho \veryshortarrow  f_\Co f_\Ho}$ one has $|E_{i_{\rm C}} - E_{f_{\rm C}} -( E_{f_{\rm H}} - E_{i_{\rm H}})| \leq \epsilon$, almost the same derivation allows us to generalise Eq.~\eqref{eq:xftbound} to
\begin{equation}
p^{\rm W} \textrm{ nonnegative} \Rightarrow Q \leq -\frac{\ave{\Delta	I}_{\rm W}}{\Delta \beta}+\frac{\log(1+\bar{\chi})}{\Delta \beta} - \frac{\beta_{{\rm H}}\epsilon}{\Delta \beta}.
\end{equation}
A symmetric inequality can be derived for the energy flowing into ${\rm C}$.

 \subsection{Case study: heat exchange between two qutrits}
\label{subsec:case_study}
In this section we investigate the heat flows between two qutrits and show that Inequality~\ref{th:xftbounds} can witness negativity in these higher dimensional systems as well.
The Hamiltonian of the qutrits is assumed to take the form $H_{{\rm C}} = H_{{\rm H}} = \sum_{n=0}^{2} E_n \ketbra{n}{n}$ with
  $E_0 = 0$, and  no degeneracy of the energy gaps (the `Bohr spectrum').
The state  $\rho_{{\rm C} {\rm H}}=\sum_{i,j =0}^8 \rho_{ij}\ketbra{i}{j}$, can be expressed by the elements  $\rho_{ij}=\rho^{*}_{ji}=\eta_{ij}e^{-i\xi_{ij}}\sqrt{\rho_i \rho_j} \quad \forall \, i\neq j$,  the population $ \rho_{i} \equiv \rho_{ii}$, and the parameters  \mbox{$\eta_{ij} \in [0,1]$} and $\xi_{ij}\in \mathbb{R}$. Here, we used the natural labeling $(00,01,02,10,11,12,\dots) \equiv (0,1,2,3,4,5,\dots)$.

 By imposing the constraint that the marginal states $\tr{{\rm C}}{\rho_{{\rm C}{\rm H}}}$ and $\tr{{\rm H}}{\rho_{{\rm C}{\rm H}}}$  are thermal, and following the same reasoning as in the two qubit scenario, heat flows between two qutrits are determined by: the six parameters $\eta_{ij}$, $\xi_{ij}$ for $(i,j) = \{(01,10) \equiv (1,3), (02,20) \equiv (2,6) ,(12,21) \equiv (5,7)  \}$; and four undetermined populations $\rho_i$ that must comply with the nonnegativity of $\rho_{{\rm C}{\rm H}}$.
 Furthermore, all possible heat flows achievable by energy preserving unitaries can be realised by $U = \oplus_{i=1}^3 U^{(i)}$, where $U^{(i)}$ are real $2$d rotations within the three above-mentioned manifolds, fixed by three angles $\theta_{01} $, $\theta_{02}$, $\theta_{12}$.
In App.~\ref{app:C} we provide the explicit structure of the state $\rho_{{\rm C} {\rm H}}$ and the energy preserving unitary.
In Fig.~\ref{fig:qutrits}, we plot the regions in parameter space $\theta_{01}$ and $\theta_{02}$ in which these probabilities turn negative (all shaded areas).
The yellow shaded area indicates the regions where Inequality \ref{th:xftbounds} is being violated.
 This clearly shows that our method allows to detect negativities also in higher dimension.
A more detailed analysis and  a straightforward generalization to arbitrary dimension is presented in App.~\ref{app:C}.

\begin{figure}[t!]
\includegraphics[width=0.8\linewidth]{figqutrits4.pdf}
	\caption{
 Negativity of the  quasiprobabilities (all shaded areas) and the violation of Inequality~\ref{thm:newbound} (yellow shaded area) for different interaction protocols determined by $\theta_{01} $ and $\theta_{02} $. Parameters: $\beta_{{\rm H}}=0.3$, $\beta_{{\rm C}}=1.3$, $\theta_{12}=\theta_{02}$, $E_1=1$, $E_2=1.15$,  $\xi_{ij} = 0$, $\eta_{ij}= 1$ for $(i,j) = \{(1,3), (2,6), (5,7)  \}$, $\rho_0=0.3$, $\rho_5 =0.03$, $\rho_7 =0.07$, and $\rho_8 =0.06$.}
	\label{fig:qutrits}
\end{figure}

 \section{Witnessing negativity with projective energy measurements in arbitrary dimension}
\label{sec:highD}

Testing violations of Eq.~\eqref{eq:xftbound} is a way to probe the fully quantum regime in arbitrary finite-dimensional systems. However, it has a significant drawback of being  experimentally very demanding, as $\bar{\chi}$ depends on all the entries in the initial state, as well as the details of the dynamics.
Standard fluctuation theorem experiments, on the other hand, involve only projective energy measurements. The drawback of these is the inability to probe quantum effects due to initial quantum correlations and energy coherence.
Can our method be adapted to obtain a scheme satisfying both desiderata, that is 1. probing the fully quantum regime 2. using only the experimental capabilities required by the traditional TPM quantum fluctuation theorems? Here we answer this question in the positive.

We consider two $d$-dimensional systems and for simplicity we assume there are no matching gaps in their energy spectrum (formally, `non-degenerate Bohr spectrum'). The most general energy-preserving unitary is then a direct product of two dimensional rotations. This allow us to derive the following result (see App.\ref{app:C} for the proof)
\begin{theorem}
	\label{thm:newbound}
	Consider two arbitrary finite-dimensional systems with non-degenerate Bohr spectrum undergoing an energy-preserving unitary. If $p^{\rm W}$ is nonnegative then
	\begin{equation}
	\label{eq:newbound}
Q \in \left[Q^{\rm TPM} -2 \lambda_-, Q^{\rm TPM} + 2\lambda_+ \right],
	\end{equation}
	where
	\begin{align}
	\lambda_- = &  \sum_{E_{i_{\rm C}} > E_{f_{\rm C}}} p^{\rm TPM}_{i_{\rm C} i_{\rm H} \rightarrow f_{\rm C} f_{\rm H}} \Delta E_{i_{\rm C} f_{{\rm C}}}, \\ \lambda_+ = &  \sum_{E_{i_{\rm C}} < E_{f_{\rm C}}} p^{\rm TPM}_{i_{\rm C} i_{\rm H} \rightarrow f_{\rm C} f_{\rm H}} \Delta E_{f_{\rm C} i_{{\rm C}}}.
	\end{align}
	and $\Delta E_{nm} := E_n -E_m$.
\end{theorem}

Remarkably, to confirm a violation of Eq.~\eqref{eq:newbound} and witness negativity, it suffices to have the experimental capability of performing projective energy measurements. In fact, $\lambda_-$ and $\lambda_+$ and $Q^{\rm TPM}$ can all be inferred from the data already required to test standard TPM fluctuation theorems. To infer $Q$, one can prepare $N$ samples of the initial state and use half to estimate the average energy before the dynamics and half to estimate the average energy after the dynamics. The difference between the two is the heat $Q$. In other words,  if one can protect the initial state from decoherence,  negativity can be witnessed in any of the platforms used to test TPM fluctuation theorems.  Its worth noting that Inequality~\ref{thm:newbound} does not require local thermality. In other words, it can also be tested measuring energy exchange between two quantum systems in a generic state.

To illustrate this result we consider the two qutrits case study of Sec.~\ref{subsec:case_study}.
Fig.~\ref{fig:qutrits2} presents the negativity of the quasiprobability related to heat exchange (all shaded area), and the violation of the Inequality~\ref{thm:newbound} (green and yellow areas) for different interaction protocols determined by the unitary angles $\theta_{01}$ and $\theta_{02}$.
The different insets represent different initial states, determined by the parameter $\eta=\eta_{ij}$ for  $(i,j) = \{(1,3), (2,6), (5,7) \}$ that quantifies the amount of initial quantum correlations between the states.
The yellow area represents violation of the lower bound of Inequality~\ref{thm:newbound} and corresponds to negativity in the direct flow as illustrated in the
quantum contributions of
Fig.~\ref{fig:heatflows} (upper right process, red arrows). The green area, on the other hand, represents violations of the upper bound of Inequality~\ref{thm:newbound}, that corresponds to negativity in the backflow as illustrated in the quantum contributions of Fig.~\ref{fig:heatflows} (lower right process, blue arrows).
The upper left inset of Fig.~(\ref{fig:qutrits2}) with $\eta=1$ corresponds to the parameters used for Fig.~\ref{fig:qutrits}. This indicates that Inequalities~\ref{th:xftbounds} and \ref{thm:newbound} can witness negativity in different regimes.  As already remarked, however, testing Inequality~\ref{thm:newbound} is much simpler.

\begin{figure}[t!]
\includegraphics[width=1.03\linewidth]{highD.pdf}
	\caption{
	Negativity of the  quasiprobabilities (all shaded areas) and the violation of Inequality \ref{thm:newbound} (yellow shaded area) for different interaction protocols determined by $\theta_{01} $ and $\theta_{02} $.
The insets shows different initial states, determined by the parameter $\eta=\eta_{ij}$ for  $(i,j) = \{(1,3), (2,6), (5,7) \}$
Parameters: $\beta_{{\rm H}}=0.3$, $\beta_{{\rm C}}=1.3$, $\theta_{12}=\theta_{02}$, $E_1=1$, $E_2=1.15$,  $\xi_{ij} = 0$, $\rho_0=0.3$, $\rho_5 =0.03$, $\rho_7 =0.07$, and $\rho_8 =0.06$.}
	\label{fig:qutrits2}
\end{figure}

\section{Outlook}
\label{sec:outlook}
 Recent no-go results~\cite{perarnau2017no, lostaglio2018quantum} strongly suggest that, if we are to probe the fully quantum thermodynamic regime, due to noncommutativity we need to renounce a straightforward statistical interpretation of microscopic fluctuation processes. Moving in this direction, here we introduced a natural quasiprobability for heat fluctuations whose negativity captures nonclassical effects in the strong form of contextuality. This formally implies the impossibility of describing the weak measurement scheme probing $p^{\rm W}$ as the result of any underlying classical stochastic process.  Differently from previous approaches, we showed that negativities in the heat fluctuations have a thermodynamic interpretation and that certain heat flows between two locally thermal states at different temperatures can only happen in their presence.  To prove the latter, we developed families of heat flow inequalities which can be only violated in the presence of negativities. We used these tools to witness nonclassicality by analysing data collected in a recent two-qubit experiment. We also presented heat flow inequalities for generic finite-dimensional systems whose testing only requires the control already needed to test standard TPM fluctuation theorems, with the advantage that one is testing the fully quantum regime.

 Moving forward, we propose that a quasiprobability approach to fluctuations can advance our understanding of the thermodynamic behavior of realistic devices and the study of heat flows in more complex networks. Classically unachievable heat flow configurations can also be used as a witness of nonclassicality in thermal machines and could unlock performance boosts in thermodynamic protocols.  We also expect the same ideas developed here can  be applied  to work fluctuations as well. Building on the seminal paper of Allahverdyan \cite{allahverdyan2014nonequilibrium}; this may provide the most direct experimental test of the framework developed here.

 More broadly, this work also opens up the possibility of applying to nonequilibrium thermodynamics the extensive technical tools based on quasiprobability representations developed to identify quantum effects and advantages in quantum optics, quantum foundations, quantum computing and condensed matter physics.

\begin{acknowledgments}
We thank David Jennings and Karen Hovhannisyan for helpful discussions. ML thanks Eran Rabani for the kind hospitality during a visit to UC Berkeley, where part of this work was completed.  ML acknowledges financial support from the the European Union's Marie Sklodowska-Curie individual Fellowships (H2020-MSCA-IF-2017, GA794842), Spanish MINECO (Severo Ochoa SEV-2015-0522 and project QIBEQI FIS2016-80773-P), Fundacio Cellex and Generalitat de Catalunya (CERCA Programme and SGR 1381). This work is part of the 'Photonics at Thermodynamic Limits' Energy Frontier Research Center funded by the U.S. Department of Energy, Office of Science, Office of Basic Energy Sciences under Award Number DE-SC0019140.
\end{acknowledgments}