Quasiprobabilities in quantum thermodynamics and many-body systems

In this tutorial, we present the definition, interpretation and properties of some of the main quasiprobabilities that can describe the statistics of measurement outcomes evaluated at two or more times. Such statistics incorporate the incompatibility of the measurement observables and the state of the measured quantum system. We particularly focus on Kirkwood-Dirac quasiprobabilities and related distributions. We also discuss techniques to experimentally access a quasiprobability distribution, ranging from the weak two-point measurement scheme, to a Ramsey-like interferometric scheme and procedures assisted by an external detector. We illustrate the use of quasiprobabilities in quantum thermodynamics to describe the quantum statistics of work and heat, and to explain anomalies in the energy exchanges entailed by a given thermodynamic transformation. On the one hand, in work protocols, we show how absorbed energy can be converted to extractable work and vice versa due to Hamiltonian incompatibility at distinct times. On the other hand, in exchange processes between two quantum systems initially at different temperatures, we explain how quantum correlations in their initial state may induce cold-to-hot energy exchanges, which are unnatural between any pair of equilibrium nondriven systems. We conclude the tutorial by giving simple examples where quasiprobabilities are applied to many-body systems: scrambling of quantum information, sensitivity to local perturbations, and quantum work statistics in the quenched dynamics of models that can be mapped onto systems of free fermions, for instance the Ising model with a transverse field. We meticulously present derivations of essential concepts alongside straightforward examples, aiming to enhance comprehension and facilitate learning.

In this paper we provide a tutorial that delves into the use of quasiprobabilities and their associated distributions within the realm of quantum science and with specific applications in quantum thermodynamics and manybody quantum systems.
As evident from numerous studies [1][2][3][4][5][6][7][8][9][10][11][12][13][14][15][16][17][18][19], quasiprobabilities have garnered significant interest within the quantum community.In fact, they are a proper tool to describe the statistics of the outcomes resulting from consecutive events in several areas of quantum mechanics and associated technologies, from foundations to quantum devices.The most appealing feature of a quasiprobability is its capability to embody the incompatibility of quantum observables, due to their nonzero commutator, that are evaluated at different times of a given experiment.This incompatibility is singled out by the fact that the distribution of the measurement outcomes admits nonclassical traits.The latter ones, expressed in the form of 'negative' and 'nonreal probabilities', reflect the effects entailed by the Heisenberg uncertainty principle that indeed concerns the impossibility of concurrently measuring two complementary and incompatible properties of a quantum system in contiguous times.
In the following, we are going to introduce the concept of quasiprobability from theoretical arguments, and we then outline their main properties, especially those concerning the loss of positivity.After that, we will address measurement procedures that allow observation of the 'negative' and 'nonreal probabilities' of pairs of measurement outcomes, by showing that such genuinely quantum features have a clear physical interpretation inherently related to quantum coherence and correlations [20,21].We are particularly interested in pedagogically explaining those physical interpretations that can be linked to thermodynamic quantities under outof-equilibrium conditions (work, heat and their distributions) [5,8,12,15,[22][23][24][25][26][27], and in the rich arena of manybody quantum systems [4,6,11,14,16,17,28].In this regard, interpretations of the well-known Loschmidt echo and out-of-time-ordered correlators (OTOCs) for manybody systems in terms of quasiprobabilities are discussed and illustrated with step-by-step worked examples.Finally, the tutorial is concluded with a discussion containing some perspectives on possible theoretical studies and experimental observations of quasiprobabilities in current quantum platforms [12,19,25,29].
This tutorial is targeted at graduate students and researchers, both theoretical and experimental, with a basic knowledge of quantum mechanics.Specifically, we aim to provide both the formal definitions that lay the foundations for quasiprobabilities in quantum science, and simple analytical examples helping to understand the fundamental concepts and applications of the method-ology.We would also give some hints on new directions on the use of quasiprobabilities that have not yet been explored, especially in quantum thermodynamics and many-body quantum systems.

II. QUASIPROBABILITIES
In this section, we introduce the concept of quasiprobabilities in quantum mechanics, following the seminal papers by Kirkwood [30] and Dirac [31] in the 1930s and 1940s, respectively.There have been several approaches that brought to light the notion of a 'negative probability' and crucially even probabilities represented by a complex number.In this tutorial, we choose to approach this topic from the fundamental standpoint of joint measurability in quantum mechanics.
Under conditions we are going to detail, Kirkwood-Dirac quasiprobabilities (KDQ) can take both negative and imaginary values.'Negative probability' and even 'probabilities represented by a complex number' can be explained from the fundamental standpoint of joint measurability in quantum mechanics.
For this purpose, let us set the theoretical framework.First, consider a quantum state preparation that generates a generic density operator ρ that, by definition, is a Hermitian, semi-definite operator with trace 1.Then, we define two quantum observables, associated to two Hermitian operators, O 1 (t 1 ) and O 2 (t 2 ), that we measure at two distinct times t 1 and t 2 with t 1 < t 2 .The observables can be generally expressed using their spectral decomposition as in terms of their eigenvalues o si (t i ) and the associated projectors Π si (t i ) onto the corresponding eigenspace.
In the general case, ρ, O 1 (t 1 ) and O 2 (t 2 ) do not commute with each other.Moreover, in the time interval [t 1 , t 2 ]-after the first measurement of O 1 (t 1 ) and before the second measurement of O 2 (t 2 )-the state of the quantum system can be subject to a generic quantum process described by a completely positive trace preserving (CPTP) quantum map Φ : ρ(t 1 ) → Φ[ρ(t 1 )] = ρ(t 2 ), which operates on and returns density operators.For simplicity, we omit the explicit time dependence of Φ to avoid making the notation too heavy.We introduce its Kraus representation such that Φ[ρ] = α K α ρK † α and α K † α K α = I, where I is the identity operator.The reader interested in the main properties of quantum maps can refer to Refs.[32,33].
In this section, we discuss how to characterize the statistics of measurements outcomes from the two-time evaluation of O 1 (t 1 ) and O 2 (t 2 ), by also taking into account the noncommutativity of the involved operators, i.e., the initial density operator ρ and the two observables.We explicitly wrote "two-time evaluation" in order to clearly distinguish it from the wording sequential measurements.In fact, as we will explain in a while, a procedure based on sequential measurements is necessarily invasive.As a result, (a) the measured system is perturbed; (b) initial quantum coherence in ρ (with respect to the basis decomposing O 1 (t 1 )) is destroyed; (c) the statistics of the outcome pairs (o s1 , o s2 ) resulting from sequentially measuring O 1 (t 1 ) and O 2 (t 2 ) changes.
Because of these reasons, we are going to consider an approach that, even in the general noncommutative case, can return a statistics of (o s1 , o s2 ) that is exempt from the invasiveness of the measurement apparatus, and is not affected by the first measurement of a sequential procedure, at least in some statistical moments [34].This aspect should not be surprising since it is known that [35] when noncommuting observables are taken into account, there is no unique formula to describe the joint probabilities of (o s1 , o s2 ) that has both a correspondence with the classical theory of probability and at the same time is returned by a noninvasive measurement routine; see Sec.II B for more details on this aspect.Now let us discuss in greater depth the invasiveness under the joint measurability problem, by introducing the celebrated two-point measurement (TPM) scheme [36].
A. Sequential projective measurements and the TPM scheme The TPM scheme is the procedure to characterize the statistics of (o s1 , o s2 ) by means of sequential measurements, which has a correspondence with the classical theory of probability [37].Albeit the process underlying the measurement outcomes has quantum traits, the results from such a measurement procedure can also be described by a probability distribution that obeys Kolmogorov's axioms of probability.The Kolmogorov's axioms are as follows: (i) the probability of getting a measurement outcome is a nonnegative real number; (ii) the probability to measure at least one of the outcomes is 1; (iii) the probability to measure any countable sequence of mutually exclusive measurement outcomes is equal to the sum of the probabilities for each outcome.
Operationally, the two quantum observables O 1 (t 1 ) and O 2 (t 2 ) are measured at times t 1 and t 2 , respectively, and a pair of outcomes (o s1 , o s2 ) is obtained.The probability distribution associated to any outcome pair is determined by repeating the sequential measurement procedure several times.Formally, the joint probability to get (o s1 , o s2 ), according to the TPM scheme, is where ρ is the initial quantum state (prepared at time t 1 ), and Π H s2 (t 2 ) = Φ † [Π s2 (t 2 )] = α K † α Π s2 (t 2 )K α is one of the projectors of O 2 (t 2 ) evolved in the Heisenberg picture.For unitary dynamics, Π H s2 (t 2 ) = U † Π s2 (t 2 )U , where we have introduced the unitary evolution operator of the system U ≡ T exp −(i/ℏ) t2 t1 H(t)dt with ℏ denoting the reduced Planck's constant, set to 1 for simplicity from now on, T the time-ordering operator and H(t) the Hamiltonian of the system at time t.
A procedure based on sequential projective measurements is invasive as it violates the no-signaling in time condition [38].To better single out this aspect, we introduce the generic joint probability d(s 1 , s 2 ) at two times, depending on the outcome pairs (o s1 , o s2 ).We do not explicitly provide a specific expression for d(s 1 , s 2 ), as we aim to study the mathematical properties that any joint probability needs to satisfy assuming a (non)invasive measurement procedure.Thus, if applied to our casestudy, the no-signaling in time condition states that the statistics of (o s1 , o s2 ) that is returned by a noninvasive measurement apparatus must fulfill the condition In other terms, the requirement for the noninvasiveness is that, marginalizing the distribution over the outcomes s 1 of the first observable O 1 (t 1 ) at time t 1 , we recover the unperturbed single-time probability p s2 (t 2 ) associated to the outcomes of the second observable O 2 (t 2 ) at t 2 .
In this respect, noninvasiveness can be considered as a synonymous of unperturbed marginals.The validity of Eq. ( 3) is a necessary condition for macrorealism [39] and, interestingly, Eq. ( 3) can be violated even in situations where no violation of Leggett-Garg inequalities is allowed [38].The violation of the no-signaling in time condition marks the main consequence of the joint measurability problem due to the incompatibility of the involved quantum operators ρ, O 1 (t 1 ) and O 2 (t 2 ).
The question that now arises is: "What information is erased by using the TPM scheme in the attempt of attaining the statistics of (o s1 , o s2 )?" or equivalently "How invasive is a procedure based on sequential measurements?"We are going to show that the TPM scheme is noninvasive if and only if [ρ, Π s1 (t 1 )] = 0 or [Π s1 (t 1 ), Π H s2 (t 2 )] = 0. Otherwise, Eq. ( 3) would be violated.In order to determine what information is erased by the TPM scheme, let us compare the final-time probability p s2 (t 2 ) = Tr[Π H s2 (t 2 )ρ] and s1 p(s 1 , s 2 ).The latter equals where the superoperator denotes the dephasing channel, which is defined over the eigenbasis of the quantum observable O 1 (t 1 ) and is applied to the initial density operator ρ.In Eq. ( 5), p s1 (t 1 ) = Tr[Π s1 (t 1 )ρ].Hence, if we compare p s2 (t 2 ) and s1 p(s 1 , s 2 ), we can see that the first measurement of the TPM scheme erases the quantum coherence contained in ρ once projected onto the eigenbasis of O 1 (t 1 ).
It is worth noting that the initial density operator can always be linearly decomposed in the basis of O 1 (t 1 ) as where is the operator containing the off-diagonal elements of ρ, with Tr[χ] = 0.

B. No-go theorem for sequential outcomes statistics
Previously, we have outlined that a procedure based on sequential measurements fails in recovering the marginal distribution p s2 (t 2 ) from the joint TPM distribution p(s 1 , s 2 ) of the pairs (o s1 , o s2 ) with respect to the measurement outcomes o s1 of O 1 (t 1 ), see Eq. ( 4) and the related discussion.This directly violates the no-signaling in time condition Eq. ( 3), and establishes that the measurement procedure is invasive.The origin of such violation lies in the fact that at least one of the commutators [ρ, Π s1 (t 1 )] and [Π s1 (t 1 ), Π H s2 (t 2 )] is different from zero.This conclusion is related to a deeper statement summarised by the no-go theorem reported in Ref. [11], which is less restrictive than the formulation firstly proven in Ref. [41].The no-go theorem states that the following three properties cannot be valid simultaneously for any initial density operator ρ if and only if i) The probability distribution of the pairs (o s1 , o s2 ), defined by the generic joint probabilities d(s 1 , s 2 ) in the time interval [t 1 , t 2 ], obeys the Kolmogorov's axioms of the classical theory of probability.
ii) The joint probabilities d(s 1 , s 2 ) lead to unperturbed marginals: iii) The joint probabilities d(s 1 , s 2 ) are linear functions of the initial density operator ρ.Formally, this means that, given a linear combination ρ = The three properties i)-iii) are all simultaneously satisfied under the assumption of the commutative condition [Π s1 (t 1 ), Π H s2 (t 2 )] = 0, and the probability distribution that fully characterizes the statistics of the pairs (o s1 , o s2 ) is the one returned by the TPM scheme.
In the following and throughout the tutorial, we will give up the property i).As a consequence we can no longer employ sequential projective measurements to characterize the statistics of (o s1 , o s2 ).Avoiding the direct application of the TPM scheme on the quantum system under scrutiny may completely eliminate the invasiveness of the measurement procedure and allow us to recover unperturbed marginal distributions [property ii)].Such a requirement for the generic joint probabilities d(s 1 , s 2 ) is well justified if we want that our knowledge on the fluctuations of the pairs (o s1 , o s2 ) is not decreased by the quantum measurement back-action.We therefore require that the no-signaling in time condition is fulfilled.Furthermore, we demand that the probability distribution of (o s1 , o s2 ) exhibits linearity, in conformity with the property iii).In this way, for any variation of ρ, one does not need to repeat from scratch the experimental procedure (which, as noted earlier, should not be sequential) to determine d(s 1 , s 2 ).
As a result, by linearly decomposing ρ as in Eq. ( 6), in terms of its diagonal and off-diagonal parts with respect to the eigenbasis of O 1 (t 1 ), we recover the results of the TPM scheme whenever χ = 0, with 0 denoting the matrix with all zeros.In addition, another consequence of the linearity property is that the procedure for measuring d(s 1 , s 2 ) can be independent on the initial density operator ρ, as it is customary in the classical case.
C. Beyond the two-point measurement scheme: quasiprobability approach Under the noncommutativity hypothesis [Π s1 (t 1 ), Π H s2 (t 2 )] ̸ = 0 for some pair (o s1 , o s2 ), dropping the property i) of the no-go theorem mentioned in Sec.II B allows for the introduction of a quasiprobability distribution (QD), whose terms can be nonpositive (i.e., negative real numbers or even complex numbers), albeit still summing to 1.In general, there is not a unique QD due to ordering ambiguities in how the QD is defined (see Refs. [9-11, 13, 35]).As a consequence, there are, in principle, infinite QD that are linear in the initial state ρ and lead to unperturbed marginals, at both the initial and final times t 1 and t 2 .
Let us introduce quasiprobabilities.We start from the expression for the generic joint probabilities d(s 1 , s 2 ) and assign a linear operator M (s 1 , s 2 ) to each pair (o s1 , o s2 ) of measurement outcomes.Without loss of generality, we can write: From classical probability theory, in the case of the TPM scheme [see Sec.II A], we find that which returns Eq. ( 2).Multiplying the observable Π H s2 (t 2 ) by the projector Π s1 (t 1 ) on both the left and right sides is equivalent to performing a projective measurement entailing the collapse of the measured quantity.In order to overcome this, the minimal change is to remove one projector Π s1 (t 1 ).We can set either or Substituting the linear operators M KDQ 1 (s 1 , s 2 ) or M KDQ 2 (s 1 , s 2 ) in Eq. ( 11) gives two perfectly valid KDQ.The difference of applying M KDQ 1 or M KDQ 2 on ρ is that they operate on off-diagonal terms of ρ with ex- thus meaning that the KDQ Tr[M KDQ 1 (s 1 , s 2 )ρ] and Tr[M KDQ 2 (s 1 , s 2 )ρ] differ by their imaginary parts that are opposite in sign.
Other quasiprobabilities have been considered in the literature; for example, for systems weakly interacting with a detector in order to avoid the invasiveness of the first measurement of the TPM scheme [22,24,[43][44][45][46][47].Given their affinity with MHQ and KDQ, we will discuss them in detail in Sec.III.Moreover, in Sec.II D, we will also give a short overview of alternative formulations of quasiprobabilities considered in the literature.
The quasiprobabilities defined in Eq. ( 16) and Eq. ( 19) fulfil properties ii) and iii) of the no-go theorem in Sec.II B, meaning that the no-signaling in time condition is fulfilled and the measurement procedure that allows to get a QD is independent of the initial state ρ.Moreover, the TPM statistics is recovered in the case in which ρ, O 1 (t 1 ), O 2 (t 2 ) all commute with each other.
These characteristics are important requisites to build a consistent thermodynamic theory via quasiprobabilities.In particular, the property iii) of the no-go theorem, regarding the linearity on the initial state, is meaningful when drawing a correspondence with classical thermodynamics.In fact, any nonlinear dependence of thermodynamic quantities-work, heat and entropy-on the initial state seems to be in contradiction with their standard definition in classical thermodynamics that does not change depending on the way the phase-space distribution taken as the input ensemble is split.

nonpositivity
KDQ naturally encode temporal correlations between the measurement outcomes of the quantum observables O 1 (t 1 ) and O 2 (t 2 ).As explained in Sec.II B, in relation to the no-go theorem, the quasiprobabilities q(s 1 , s 2 ) can be nonpositive, although they are still subject to the normalisation constraint In the case of KDQ, nonpositivity can mean the following two facts: I) the real part of q(s 1 , s 2 ) is negative; II) q(s 1 , s 2 ) is a complex number with a nonzero imaginary part.
From a mathematical point of view, the onset of nonpositivity in KDQ is a consequence of the fact that the product of two quantum observables (say, A and B) can always be decomposed as the linear combination of selfadjoint operators as AB = {A, B}/2+i[A, B]/(2i).Thus, for the product Π H s2 (t 2 )Π s1 (t 1 ), we have The first term on the right-hand-side of Eq. ( 21) gives rise to the MHQ-the real part of the KDQ-when evaluated (i.e., averaged) with respect to the initial density operator ρ.Then, looking at the second term of Eq. ( 21), it is evident that a necessary condition for the KDQ to have an imaginary part is the noncommutativity of Π s1 (t 1 ) and Π H s2 (t 2 ).In this regard, we stress , and the KDQ coincides with the TPM probabilities.In order to ensure the validity of Eq. ( 20), the imaginary parts of KDQ must cancel each other out.
In Sec.II C 2, we will show a simple case study that directly connects the imaginary parts of the KDQ with the presence of imaginary coherence terms in the initial density operator ρ, with respect to the eigenbasis of O 1 (t 1 ).Hence, if we ignored the imaginary parts of q(s 1 , s 2 ), we would exclude some information stemming from quantum coherence and correlations that may emerge in the quantum statistics of (o s1 , o s2 ).
The occurrence of nonpositivity can be considered as a nonclassical feature in the statistics of the measurement pairs, underlining the presence of genuinely quantum features due to the interplay of quantum dynamics and measurement.From here on, we will refer to this with the term nonclassicality.The formal conditions and experimental routines allowing to identify nonclassicality take also the name of quantum contextuality [1].In a scenario involving a quasiprobability distribution describing the occurrence of a given pair of measurement outcomes at two times, an experimental protocol is defined contextual if it is able to yield nonpositive values, provided incompatibility of noncommuting operators plays a role.No classical stochastic process can explain such a behaviour.
a. noncommutativity is a necessary condition.For the MHQ q MHQ , see Eq. (19), it has been recently shown that the pairwise noncommutativity of the initial density operator and the quantum measurement observables is only a necessary but not sufficient condition for nonpositivity (i.e., negativity of the MHQ).This means that there are counterexamples where [ρ, Π s1 (t 1 )] ̸ = 0 and/or [Π s1 (t 1 ), Π H s2 (t 2 )] ̸ = 0, but still Re [q(s 1 , s 2 )] ≥ 0 [48].A detailed analysis of this aspect can be found in Ref. [9], where it has been formulated with the wording negativity is stronger than noncommutativity.
b. Direct link with weak values.The nonpositivity of KDQ can find a physical interpretation from their direct connection with weak values [49][50][51] that, indeed, are conditional KDQ [11] (24) We recall that the weak values can be obtained via a weak measurement that is performed on both a properly chosen preselected quantum state and a post-selected one.Weak values are called anomalous when ⟨O 1 (t 1 )⟩ WV lies outside the spectrum of O 1 (t 1 ) (see Refs. [52][53][54][55]).In order to guarantee such anomaly, it is required that the (possibly mixed) pre and postselection states have quantum coherence with respect to the eigenbasis of O 1 (t 1 ) and that the corresponding KDQ exhibits negativity [56].Moreover, the generalization of weak values to mixed density operators (instead of pure quantum states), can be a complex number, and this is evidently in a one-toone correspondence with complex KDQs [52,54,57,58].Overall, the occurrence of an anomalous weak value is identified by the nonpositivity of a KDQ, implying quantum contextuality [53,54].
c. nonpositivity functional.We conclude this subsection by introducing the nonpositivity functional [9,11,12,59] that quantifies the 'amount' of nonclassicality in the statistics of the outcome pairs (o s1 , o s2 ).It is worth noting that both the real and imaginary parts of the KDQ contribute to the nonclassicality, whereby if present one has that As noticed in Ref. [9], one could quantify the negativity and nonreality of a KDQ distribution by using the functionals Re [q(s 1 , s 2 )] ( 27) that act separately on the real and imaginary parts of KDQ, respectively.The condition ℵ = 0 occurs when all the KDQ are positive real numbers [60].

Comparing KDQ and TPM probabilities
In this subsection we are going to compare the KDQ q(s 1 , s 2 ) with the joint probabilities p(s 1 , s 2 ) returned by applying the TPM scheme.In this regard, notice that is the projector orthogonal to Π s1 (t 1 ).Interestingly, as discussed in Ref. [35], Eq. ( 29) quantifies the interference patterns between the two different sequential pairs of projectors, also known in the literature as quantum histories, (Π s1 (t 1 ), Π H s2 (t 2 )) and (Π ⊥ s1 (t 1 ), Π H s2 (t 2 )).Moreover, the right-hand-side of Eq. ( 29) is also recovered from the so-called nondemolition quasiprobability (NDQP) [22,24,61] with The NDQP is evidently defined over three indexes: two, s 1 and s ′ 1 (different each other), refer to two possible measurement outcomes of the quantum observable O 1 (t 1 ) at time t 1 , while s 2 refers to O 2 (t 2 ) as it holds for q(s 1 , s 2 ).Thus, by marginalizing over s ′ 1 ̸ = s 1 , one directly obtains the difference between the KDQ and TPM (joint) probabilities: It can be easily observed that if Re [q(s 1 , s 2 )] < 0, then necessarily s ′ 1 ̸ =s1 Re [q(s 1 , s ′ 1 , s 2 )] < 0; moreover, when the KDQ q(s 1 , s 2 ) is a complex number, also ) is a complex number with the same imaginary part of q(s 1 , s 2 ).
Let us now exemplify these concepts with a simple case study.We consider a spin-1/2 particle, first initialized in the generic density operator ρ.Then, the spin of the particle is consecutively measured along two orthogonal axes, the z and x axis, respectively, i.e.O 1 (t 1 ) = σ z and O 2 (t 2 ) = σ x .Moreover, we assume U = I, so that the system does not evolve between t 1 and t 2 .Note that, under specific conditions [62], this setup is representative of the physics underlying the well-known Stern-Gerlach experiments.
At the end of this quantum process, the state of the system collapses onto one of the eigenstates of σ x − = −1, respectively.The quantum process is inherently probabilistic, due to the stochastic nature of any quantum measurement.We thus need to calculate the probabilities of obtaining the pairs of measurement outcomes (z k (t 1 ), x j (t 2 )) measured at times t 1 and t 2 with k ∈ {0, 1} and j ∈ {−, +}, for the initial density operator ρ.
As previously anticipated, if [ρ, O 1 (t 1 )] ̸ = 0, then the application of the TPM scheme (i.e., sequential projective measurements) does no longer suffice.This fact is confirmed by the direct computation of the differences q(z k , x j ) − p(z k , x j ): where k,j p(z k , x j ) = k,j q(z k , x j ) = 1 by construction.
From Eqs. ( 33)- (36), it is also apparent that at least two of the differences q(z k , x j )−p(z k , x j ), among all four, exhibit negative real parts whenever the initial state ρ does not commute with the quantum observable σ z at time t 1 .Notably, such a negativity is preserved from applying a second measurement of the observable σ x immediately after the first.
Of course, in the case ρ 0,1 = 0, the KDQ q(z k , x j ) reduces to the TPM joint probabilities p(z k , x j ), and the no-signaling in time condition is fulfilled.In contrast, in the case ρ 0,1 ̸ = 0, the first measurement of σ z required by the TPM scheme turns out to be invasive for the joint statistics of the measurement outcomes (z k (t 1 ), x j (t 2 )).
We now provide the average of the difference of outcomes ∆o = x(t 2 ) − z(t 1 ) (thus, ∆o j,k = x j (t 2 ) − z k (t 1 )) that is evaluated with respect to the KDQ q(z k , x j ).We have By setting ρ = I/2, it holds that ⟨∆o⟩ = 0 that stems from having all the KDQ equal to 1/4.This finding is in accordance with the classical theory of probability applied to our case study.In fact, if the initial density operator of the spin-1/2 is mixed with both elements equal to 1/2 (i.e., the spin of the particle is initially up or down with equal probability 1/2), then the sequence of measurement outcomes ±1 obtained from applying two mutually uncorrelated operations (i.e., the sequential projective measurement of σ z and σ x ) is naturally equiprobable.As a result, on average the difference of the measurement outcomes ∆o is zero.
Let us now add quantum coherence to the initial state ρ with respect to the eigenbasis of σ z , by taking with χ defined in Eq. ( 8).Hence, from Eq. (37), we obtain ⟨∆o⟩ = 2Re [ρ 0,1 ], meaning that a correction to the "classical" result ⟨∆o⟩ = 0 has to be included.In this case study, such a correction is directly proportional to the quantum coherence of ρ.
From the previous discussion, it is evident that the nonclassicality due to the incompatibility of non commuting operators is very fragile and can be erased by external noise and dissipation in an open quantum system.To show this, let us assume that the initial state of the qubit at time t 1 is ρ and that before the second measurement the qubit undergoes a pure-dephasing (PD) channel.The latter is described by the following map (superoperator acting on ρ): which realizes a phase flip occurring with probability p.If we repeat the calculations leading to Eqs. ( 33)- (36), then the quantum coherence in the resulting new equations is reduced by a factor 1 − 2p as an effect of the action of the noisy channel.In particular, one gets:

Distribution and characteristic function of KDQ
As mentioned in the previous sections, the KDQ q(s 1 , s 2 ) describes the joint probability of the outcomes' pairs (o s1 , o s2 ) from measuring the quantum observables O 1 (t 1 ) and O 2 (t 2 ) at times t 1 and t 2 , initial and final times of the quantum process in analysis, with t 1 < t 2 .The individual outcomes s 1 and s 2 correspond to the eigenvalues of the observables in Eq. ( 1).
Let us introduce the generic difference of outcomes such that ∆o s1,s2 = o s2 (t 2 ) − o s1 (t 1 ).The number of values that ∆o can take depends on the combinations of all possible measurement outcomes at t 1 , t 2 .Therefore, the KDQ distribution of ∆o is defined by where δ(•) is the Dirac δ function.We remark that the KDQ distribution P [∆o] is not unique due to ordering ambiguities entailed by the noncommutativity of ρ, O 1 (t 1 ) and O 2 (t 2 ), as discussed in Sec.II C. We also note that the distribution of ∆o provided by the TPM scheme is where, as before, p(s 1 , s 2 ) denotes the TPM joint probabilities.
All the information about the statistics of the outcome pairs (o s1 , o s2 ) is also encoded in the characteristic function of P [∆o] defined as its Fourier transform: While in principle for a Fourier transform the variable u is real, it may be useful to extend Eq. ( 47) with u as a complex number, as we will see in Sec.IV B. Interestingly, both the KDQ q(s 1 , s 2 ) and the characteristic function G(u) are quantum correlation functions, namely they can be obtained as the expectation value of the product of two operators (not necessarily Hermitian, but defined at two times) on the initial density operator ρ.In general, the distribution P [∆o] depends on the time duration of the quantum system dynamics.Of course, the time dependence of P [∆o] is mirrored in a time-dependent characteristic function G(u).Both for P [∆o] and G(u), the time-dependence is omitted, unless specified, to enhance the clarity of the presentation.
The characteristic function can be used to detect complex values in the KDQs.In fact, a violation of the equality G(−u) = G(u) * implies that the positive-semidefinite condition for q(s 1 , s 2 ) does no longer hold [11].The equality G(−u) = G(u) * is violated only when Im[q(s 1 , s 2 )] ̸ = 0, and such a violation serves as a witness of complex values in the KDQs, which are identified by the nonpositivity functional.

D. Alternative formulations
We conclude this section of the tutorial, dedicated to the introduction of quasiprobabilities, by mentioning other possible formulations, alternative to the ones already discussed before, able to identify the presence of nonclassical temporal correlations.
The first formulation to consider (even for historical reasons) are phase-space distributions.In quantum mechanics, particularly in quantum optics, phase-space distributions represent the state of a quantum system, e.g., a light mode or a spin ensemble, using quasiprobability distributions like the Wigner function [63].These distributions combine position and momentum, providing a comprehensive framework to analyze quantum states, coherence, and entanglement in optical systems [64][65][66][67][68][69].It is important to remark that phase-space distributions describe the state of a quantum system rather than a quantum process involving two measurements at distinct time, as we have described in this tutorial.Negativity in phase-space distributions signals quantum coherence in the state and has been used as a quantifier of nonclassical states.
It is also worth mentioning the Keldysh-ordered full counting statistics (FCS) [70][71][72][73].It has been originally introduced to study the fluctuations of a time-integrated quantum observable, like current fluctuations in quantum electronic conductors averaged over a given time interval.Formally, if O H (t) is an observable in the Heisenberg picture at time t, the FCS corresponds to the probability distribution P τ [I] of the time integral Notice that the probability distribution of I depends on the time interval τ .In a similar fashion of a KDQ or NDQP distribution, P τ [I] can be obtained by means of the inverse Fourier transform of the moment generating function for I, namely where H λ ≡ H + λI/2 with H denoting the Hamiltonian of the quantum system that is initialized in ρ.The quasiprobability distribution P τ [I] is not always positive, signaling nonclassicality [73].In fact, from the unravelling of the FCS in the spirit of Feynman's path-integral approach, the negativity in P τ [I] results from the interference of the products of probability amplitudes that comprise the distribution P τ [I].In this context, each product of probability amplitudes forms a discrete trajectory that is weighted by the elements of the initial density operator corresponding to the initial position of each trajectory.If this kind of interference between pairs of trajectories is absent, then P τ [I] is positive definite.The unravelling of the probability distribution P τ [I] [73] has marked similarities with the derivation of the nondemolition quasiprobability distribution [22,61] in Sec.II C 2 that, we recall, identifies the interference between pairs of sequential measurement projectors.Moreover, the unravelling of the FCS is also at the basis of the so-called Keldysh quasiprobabilities [74,75].In fact, the Keldysh formalism applied for investigating fluctuations of noncommuting operators extends both the phase-space distributions in quantum optics and the FCS of a timeintegrated quantum observable.This is evident from the fact that the Keldysh quasiprobability distribution coincides with the Wigner function when position and momentum operators are considered, and it reduces to the FCS when we are interested in an observable integrated over time.Interestingly, what allows for this extension is to account for the back-action exerted by a detector measuring the two-commuting observables.Also this feature is shared by nondemolition quasiprobabilities that can be attained using a detector-assisted measurement scheme, as we will detail in Sec.III C. We conclude by noting that, similarly to the KDQ, the negativity of a Keldysh quasiprobability distribution requires the noncommutativity of the initial density operator with the first measurement observable, or the noncommutativity of the operators measured at distinct times.

III. MEASURING QUASIPROBABILITIES
In this section, we are going to present two approaches that allows the reconstruction of a QD: the first is based on performing only projective measurements [12,76], while the second (based on interferometry or assisted by a detector) is aimed at measuring directly the characteristic function of the QD under scrutiny.More than these two approaches have been formulated so far to achieve such a reconstruction [77][78][79][80][81][82]; the reader can find more details in Ref. [11] where alternative methods have been surveyed and some of them extended.Moreover, it is also worth mentioning Ref. [83] that investigates the use of quantum circuits for the measurement of weak values and KDQ distributions.

A. Weak two-point measurement scheme
The real part of the KDQ distribution, defined in Sec.II C as the Margenau-Hill (MH) distribution, can be determined by resorting only to a scheme entirely based on projective measurements.We have already proved that a procedure directly using sequential projective measurements cannot carry out this task.Instead, the combination of projective measurement schemes accomplishes the task.This is indeed enabled by the weak two-point measurement (wTPM) scheme for the measurement of quantum time correlators [11,76].The main feature of the wTPM scheme is to cancel the measurement back-action, thus attaining the back-action-free limit and restoring a condition of no measurement invasiveness [3].
As noticed in Ref. [12], the wTPM scheme can be effectively seen as a probabilistic error cancellation technique, a technique largely employed in quantum computing from sampling noisy circuits [84].
Let us consider the MHQ q MHQ (s )ρ] and the wTPM probability: where Π ⊥ s1 (t 1 ) has been defined in Eq. ( 30).The wTPM probability has a clear physical meaning and can be obtained via a proper measurement procedure.In fact, the transformation is associated to the events "the outcome o s1 is recorded" or "the outcome o s1 is not recorded", both at the initial time t 1 .For this reason, being given by a binary measurement result, the transformation Eq. ( 50) is denoted as nonselective measurement, and applies to a given projector of the quantum observable of interest-in this case, the projector Π s1 (t 1 ) of O 1 (t 1 ).We introduced the wTPM probability because one can infer the MHQ from w(s 1 , s 2 ).To see this, we just need to substitute Eq. ( 30) in Eq. ( 49), and write the explicit expression of w(s 1 , s 2 ) as a function of Π s1 (t 1 ); we get: (51) with the result that Eq. ( 52) is the way the MHQ can be experimentally reconstructed via the wTPM scheme, as done e.g. in Ref. [12] where a pictorial representation of the scheme is provided.In fact, the TPM joint probability p(s 1 , s 2 ) can be obtained via a procedure of sequential projective measurements, and p s2 (t 2 ) is the unperturbed single-time probability to measure one of the outcomes o s2 (t 2 ) of O 2 (t 2 ) at the final time t 2 .Finally, the wTPM probability w(s 1 , s 2 ) is returned via a procedure based on nonselective projective measurements, as already explained above.
Notice that the probability p s2 (t 2 ) also enters the socalled end-point measurement (EPM) scheme [85,86] that, by construction, singles out the presence of quantum coherence in the initial state ρ by performing single measurements at the end of the quantum process under scrutiny.A discussion about the conceptual difference of the KDQ and the joint probabilities stemming from the EPM scheme can be found in Ref. [11].
We conclude this subsection by observing that, for qubits, the expression of w(s 1 , s 2 ) simplifies.This is because where D 1 [ρ] is the dephasing superoperator defined in Eq. ( 5).As a result, the wTPM probability reduces to the marginal of the TPM joint probability p(s 1 , s 2 ) over the outcomes s 1 of the initial observable, i.e., such that where χ, defined in Eq. ( 7), contains the quantum coherence in ρ.

B. Interferometric scheme
Another approach for the inference of the KDQ distribution P [∆o], which we consider in this tutorial, is an interferometric scheme.This method consists in encoding on an auxiliary system A the real and imaginary parts of the characteristic function G(u) of P [∆o] for a given quantum system S.The use of an auxiliary system allows one to infer both the real and imaginary parts of the KDQ by implementing a unique scheme.As explained in Sec.III A, if our aim is to reconstruct only the real part of a KDQ, we can resort to a procedure that is only based on projective measurements.
The interferometric scheme we are going to present here is a simplified variant of the theoretical proposals discussed in Refs.[11,[87][88][89][90], and has similarities with the experimental schemes employed in Refs.[91,92].It has been recently realized in Ref. [93] using an electronnuclear spin system associated with a nitrogen-vacancy center in diamond.However, all these interferometers lead to the same result, namely the direct measurement of the characteristic function G(u).Notably, the observed G(u) can belong to both a probability distribution stemming from a procedure of sequential projective measurements (thus, a TPM distribution), and a quasiprobability one.
When the auxiliary system A is taken as a qubit, the real and imaginary parts of G(u), can be extracted from the expectation values of two Pauli matrices with respect to the state of A at the end of the scheme.As it will be clearer later, in order to implement the interferometer, u is taken as a real number with the dimension of a time t.By collecting several values of the pairs (Re[G], Im[G]) for different u, we can reconstruct the (quasi)probability distribution P [∆o] by applying the inverse Fourier transform to G. The Fourier transform is performed numerically, and hence is subject to finitetime and finite-resolution constraints; see for example Ref. [94].
Let us present the interferometric scheme for quantum systems subject to unitary dynamics by assuming A is a qubit.The extension to open quantum systems, i.e., nonunitary dynamics, is straightforward through the substitution of the unitary operator with a CPTP map Φ, as long as the environment does not affect the auxiliary system.
As pictorially represented in Fig. 1, the working principle of the scheme is to initialize the auxiliary system A in the state |0⟩ A ⟨0| where we denote with |i⟩ A the eigenstates of the Pauli matrix σ z A for the auxiliary system A (i = 0, 1).We then perform a Ramsey interferometric scheme on A [95,96].Between the application of both the Hadamard U Had = 2 −1/2 (σ x A + σ z A ) and σ x A gates, and the final projective measurement of A (with respect Pictorial representation of the interferometric scheme to directly access the characteristic function G(u) of a KDQ distribution P [∆o].The scheme encodes the information on the real and imaginary parts of G(u) associated to the quantum system S of interest that is initialized in the generic density operator ρ.Such an encoding is operated on the auxiliary system A via the two conditional gates F C t 1 and F C t 2 , which applies the operations Et 1 and E † t 2 on S whether A is in |0⟩A (white dot in the figure) or |1⟩A (black dot) respectively.The other gates involved are the Hadamard gate U Had , the system's evolution operator U and the gate σ x A for the auxiliary system.A detector-like box on A denotes its final measurement.
to the observables σ x A , σ y A ), the auxiliary system interacts with the quantum system S via the conditional (C) gates where The latter can be thought as unitary evolution operators corresponding to the effective Hamiltonian O j (t j ) for a time u.Moreover, between the two conditional gates F C t1 and F C t2 , the quantum system S undergoes the actual unitary dynamics of the process ruled by the evolution operator U .
The result of implementing the interferometric scheme of Fig. is that the reduced state ρ ′ A of A, before the final measurement of σ x A and σ y A , is In this way, we can obtain: In Sec.III D, we will illustrate how the interferometric scheme works in a simple qubit case.

C. Detector-assisted measurement of quasiprobabilities
Here we provide a more general framework for the measurement of quasiprobabilities considering a quantum model for a detector coupled to the system of interest.This framework extends the TPM and the Ramsey scheme, by realizing that the observation of the change of an observable O(t) can be attained not through von Neumann projective measurements, but rather via a generalized measurement or, more precisely, a positive operatorvalued measure (POVM).This was first introduced by Roncaglia et al. in Ref. [43] to assess the thermodynamic work done on a quantum system (see Sec. IV), while a proposal for its measurement in cold atoms was reported immediately after in Ref. [44].Moreover, an experimental realization of the POVM to measure quantum work done on a Bose-Einstein condensate is in Ref. [46].In a series of papers, Solinas and his collaborators formalized this approach establishing a profound connection between fluctuations of quantum observables, quasiprobabilities and the full counting statistics [22,24,45,47,61].
We will describe two possible schemes: one that provides access to the characteristic function of the distribution associated to ∆o, and the other providing directly the (quasi)probability distribution [22].
Let us imagine a system coupled to a detector that is modelled as a quantum free particle moving in one dimension.The detector is described by the canonical position X and momentum P operators.We assume that the system-detector interaction Hamiltonian is where the time-dependent coupling constant b(t ] is such that the detector is impulsively coupled to the system only at times t 1 and t 2 with strength κ.The operator O(t) is the observable to be measured and, without loss of generality, can be thought of being O 1 (t 1 ) at time t 1 and O 2 (t 2 ) at time t 2 , as formalized in Sec.II.
In the same spirit of what we discussed above, we consider the initial state of the system to be ρ and that of the detector to be pure, i.e., where |p⟩ are eigenstates of the momentum operator with eigenvalue p and G(p) is the initial momentum distribution of the detector.Moreover, for simplicity, let us assume the system to evolve with the unitary operator U between times t 1 and t 2 .The extension to the case of a non unitary evolution described by a CPTP map is straightforward.The detector may also evolve during these times, but the effect of this free evolution can be very small or compensated, and we will therefore ignore it [22].The final state of the detector, after tracing out the state of the system, can be cast in the following two equivalent forms, with ∆o s1,s2 = o s2 (t 2 ) − o s1 (t 1 ) [see Eq. ( 44)]: where we have used the expression q of the NDQP, see Eq. ( 31).The last two equations are written in the detector momentum and position representations respectively, whereby we have introduced the detector position distribution g(x) as the Fourier transform of G(p).Some considerations are in order: since the detector momentum P is a conserved quantity (as it commutes with the interaction Hamiltonian), the sole effect of the evolution is to induce phase shifts in the momentum eigenstates |p⟩.By measuring these phase shifts we access information about the observable change ∆o.On the other hand, the evolution operator exp{−iκP ∆o} associated with the system-detector interaction is effectively a translation or displacement operator of the detector position.Therefore, the quantities ∆o can be also observed by measuring the detector position distribution.In the following, we will describe two schemes that follow these two ideas, respectively.
In the first scheme, the detector is initially prepared in a superposition of two states with opposite momenta of magnitude p 0 : where A is a real normalization constant [97].This corresponds to a momentum distribution After the evolution, the state of the detector can be written as where information about the dynamics is contained in the phase ϕ given that P is a conserved quantity (see also considerations above).If we now measure the phase ϕ accumulated during the whole evolution between the eigenstates |p 0 /2⟩ and |−p 0 /2⟩, we have access to a modified characteristic function: which resembles the characteristic function G(u) defined in Eq. ( 47), but this time computed over the NDQP of Eq. ( 31).Hence, to all effects, Eq. ( 68) represents the characteristic function of a nondemolition quasiprobability distribution.As already outlined at the level of quasiprobabilities in Sec.II C 2, G is a symmetric version of a KDQ characteristic function, where the symmetrization is done over the indexes s 1 and s ′ 1 labelling the outcomes of the initial observable O 1 (t 1 ).It is interesting to notice that Eq. ( 68) can also be obtained by the trace of Eq. ( 26) in Ref. [98].Thus, the inverse Fourier transform of G returns the NDQP distribution that is real (since q(s 1 , s ′ 1 , s 2 ) + q(s ′ 1 , s 1 , s 2 ) is real) but can assume negative values due to the non commutativity of ρ, O(t 1 ) and O(t 2 ).When the initial state of the system ρ does not have any coherence in the basis of eigenstates of O(t 1 ), then ρ s1s ′ 1 = 0 for s 1 ̸ = s ′ 1 and the inverse Fourier transform of G reduces to the TPM probability distribution, see Eq. ( 2).
In the second scheme, we analyze directly Eq. ( 65), which leads to the position probability distribution for the detector: that is real and never negative as it derives from the expectation value of a Hermitian and positive semi-definite density operator.If we assume an initially localized detector position, g(x) = δ(x), then )δ s1,s ′ 1 , and Eq. ( 70) reduces to that corresponds to the TPM probability distribution for x = κ∆o.In contrast to the case in which g(x) is delocalised, in this case there is a unique relation connecting x and ∆o allowing perfect reconstruction of the TPM probability distribution for ∆o from the statistics of the detector's position X.
Even though P (x) in Eq. ( 70) is real and positive semi-definite, effects due to initial quantum coherences can manifest themselves when the detector's initial wave function g(x) is not localized and has a width comparable or larger than the typical changes κ∆o s1,s ′ 1 .In fact, imagine that ρ s1,s ′ 1 ̸ = 0, then in Eq. ( 70) the functions g(x−κ∆o s1,s2 ) and g * (x−κ∆o s ′ 1 ,s2 ) may have an overlap that results in a modification of the position probability distribution if compared to the one provided by the TPM scheme.

D. Examples
We conclude this section by providing examples of the quasiprobability distributions obtained using the schemes presented above.

Weak-TPM scheme
Let us start with an application of the weak two-point measurement scheme to a three-level system (or a spin-1) that is initialized in a generic density operator ρ, and whose spin is consecutively measured along the orthogonal axes z and x.In particular, we take x ⟩, |ϕ (0) x ⟩, |ϕ (1)  x ⟩ = 1 2 respectively.As in the qubit example in Sec.II C 2, no quantum dynamics occurs in between the projective measurements of O 1 (t 1 ) and O 2 (t 2 ) i.e.U = I.We now write: the expression of the MHQ q MHQ (s 1 , s 2 ), the joint probability p(s 1 , s 2 ) of the TPM scheme, the unperturbed final-time probability p s2 (t 2 ) and the wTPM probability w(s 1 , s 2 ).We recall that p(s 1 , s 2 ), p s2 (t 2 ), w(s 1 , s 2 ) can be all experimentally measured via a procedure based on single or sequential measurements.Moreover, by linearly combining them together according to Eq. ( 52), any MHQ can be fully reconstructed [12].Thus, for the example under scrutiny, and w(s 1 , s 2 ) is given by Eq. ( 49) with Π H s2 (t 2 ) = |ϕ ) are projectors with rank-2 and describe the collapse of the spin-1 state onto a twodimensional subspace.For completeness, the explicit expressions of Π ⊥ s1 (t 1 ) with o s1 (t 1 ) = −1, 0, 1 are: Let us now take the initial density operator ρ = |ψ⟩⟨ψ| with In the following, we provide the analytical expressions of q MHQ (s 1 , s 2 ), p(s 1 , s 2 ), p s2 (t 2 ) and w(s 1 , s 2 ) for a single pair (s 1 , s 2 ): s 1 = −1, s 2 = 1.This choice ensures that q MHQ (−1, 1) < 0. In doing this, we will show a specific example of how Eq. ( 52) effectively works: From direct calculations, one can find that ⟨ϕ x |ϕ x | ρ |ϕ (1)  x . Therefore, Moreover, w(−1, 1) = Tr |ϕ (1)  x ⟩⟨ϕ (1)  x | ⟨ϕ that validates Eq. ( 76).

Interferometric scheme
We consider the Ramsey interferometric scheme applied to a spin-1/2 particle that is sequentially measured along the z and x axis, as in the example reported in Sec.II C 2. We recall that, in the case the initial density operator ρ does not commute with O 1 (t 1 ) = σ z , any procedure based on sequential projective measurements is invasive and unavoidably cancels out the quantum coherence contained in ρ, with respect to the eigenbasis of O 1 (t 1 ).As a result, the statistics of the outcome pairs (z k (t 1 ), x j (t 2 )) is distorted.The unperturbed statistics of the outcome pairs is provided by the KDQ distribution that, as shown in Sec.II C 2, can exhibit negative and imaginary quasiprobabilities.
The interferometric scheme finds application to such a case study by making the substitution E t1 = exp(−i σ z u), U = I, and E t2 = exp(−i σ x u).In this way, by measuring the expectation values ⟨σ x A ⟩(u) and ⟨σ y A ⟩(u), we can recover the real and imaginary parts of the KDQ characteristic function providing the statistics of outcome pairs (z k (t 1 ), x j (t 2 )).In this case, the characteristic function reads as [see also Eq. ( 47)] Using the initial state defined in Eq. ( 38), we obtain: such that the analytical expressions of the expectation values for the auxiliary system are Taking the inverse Fourier transform, we get the following QD for ∆o: which may be complex depending on the form of ρ 0,1 .
It is instructive to point out that, by defining the effective Hamiltonian operators H 1 ≡ ωσ z and H 2 ≡ ωσ x , then Eq. ( 81) takes the more general expression with u = ωt.From this, we observe that the characteristic function of the KDQ distribution can be identically equal to the so-called Loschmidt echo [99,100].Hence, thanks to the link with the Loschmidt echo, H 1 and H 2 can be interpreted as the Hamiltonian operators governing, respectively, the forward and backward evolution of a perturbed quantum system [101], and t as the time instant at which the reversal operation takes place.In Sec.V B we will show that the connection between the characteristic function of a KDQ distribution and Loschmidt echos does not hold only in specific examples, but it is valid in general for any quantum system.In this respect, condensed-matter physics and quasiprobabilities are deeply related.

Detector-assisted scheme
Let us consider the modified characteristic function G in Eq. (68).By choosing the initial state qubit to be Eq.( 38) as in III D 2, we find: which clearly depends on the presence of the off-diagonal element ρ 0,1 of the system's density operator ρ.
Consequently, taking the Fourier transform, we obtain the NDQP: which contains peaks at ∆o = ±1 that are absent in the KDQ (or for an incoherent initial state) and that can be negative depending on the sign of Re [ρ 0,1 ].Another difference with the KDQ extracted from the Ramsey scheme, Eq. ( 85), is that Eq. ( 88) is strictly real.
Let us now consider the signal observed in the position representation of the detector assuming for concreteness the detector's wave function g(x) = (2πσ 2 ) 1/4 exp −x 2 /(4σ 2 ) , albeit any localized function would be suitable.From Eq. ( 70), we obtain P (x) = P inc (x) + P coh (x), where the incoherent and coherent parts of the distributions read as The probability density P inc (x) would always appear in the expression of P (x), even in the absence of initial coherence.It represents a coarse-grained version of the TPM probability distribution; see Eq. ( 46) for the general definition.On the other hand, P coh (x) is proportional to Re [ρ 0,1 ] and can be negative (although P (x) can never be negative).We show P (x) in Fig. 2 and compare the cases when ρ 0,1 = 0 and ρ 0,1 ̸ = 0.The latter case brings an asymmetry to the function P (x) and, for sufficiently large σ, a slight shift in the position of the peaks and deeper dips in between peaks.This effect increases for larger ρ 0,1 and κ and is due to interference between probability amplitudes for the different measurement outcomes.

IV. QUANTUM THERMODYNAMICS
In Sec.II, we have argued that, in the general case of ρ and O 1 (t 1 ) arbitrary noncommuting operators, the first measurement of the TPM scheme is invasive.Specifically, the TPM scheme does not allow one to recover the unperturbed single-time probability p(s 2 ) by marginalizing the joint probability p(s 1 , s 2 ) of the TPM scheme over the outcomes s 1 of the first measurement.As a result, applying the TPM scheme breaks the no-signaling in time condition, failing to capture noncommutativity in the statistics of the measurement outcomes taken at times t 1 and t 2 [5,23].This evidence is proven to be true for arbitrary quantum observables O 1 (t 1 ) and O 2 (t 2 ).Consequently, the same considerations shall hold in a thermodynamic context, where the measured observables are Hamiltonian operators.
Kirkwood-Dirac quasiprobabilities can be employed to investigate nonclassical energetic processes, where here 'nonclassical' indicates the presence of negative and imaginary values in the quasiprobability distribution of the thermodynamic quantity of interest, for instance, work, heat or entropy.In the current literature, Margenau-Hill quasiprobability distributions [5,8,12,35,102,103], the real parts of Kirkwood-Dirac ones, have been already discussed and employed to characterize nonclassical work distributions [12,104], as well as the statistics of anomalous heat exchanges due to quantum correlations [8].
In this section, we discuss the KDQ approach to characterize the statistics of internal energy fluctuations in a generic quantum system, close or open.Then, we will focus on the ways the presence of nonclassicality is responsible for anomalous energy exchanges [25] (see below for a proper definition) that can be identified, e.g., in the average and variance of work and heat distributions.
As we are going to show in Sec.IV D, negative probabilities in a KDQ distribution of work find an interpretation as nonclassical energy transitions that make use of quantum coherence to transform absorbed energy in extractable work.In this regard, we will show how KDQ can take into account genuinely quantum features in energy-change fluctuations, and outline thermodynamic advantages.For example, this becomes evident when noting that, without quantum coherence, stochastic work processes can generate a lower amount of extractable work and thus be less performing in operating a quantum device.
Before proceeding, it is worth stressing that possible thermodynamic advantages enabled by quantum features are inherent to any quantum system subject to a thermodynamic transformation that does not spoil the state of the system with respect to the expected output.For example, a measurement procedure based on projective measurements, such as the TPM scheme that strongly perturbs the system's state, cannot be considered as a thermodynamic transformation preserving quantum features, in the same way as any decoherence process.Approaches using quasiprobabilities, instead, have been thought to characterize the quantum statistics originating from measuring the energy of a given quantum system at multiple times.Of course, to attain this characterization experimentally and possibly in real time, one has to resort to a measurement procedure that is the least invasive possible.The use of a detector weakly interacting with the measured system (see Sec. III C) could represent a solution.
Overall, acting directly on a quantum system while a given thermodynamic transformation is active leads to the spoil of functionalities entailed by quantum coherence or correlations.However, a proper quasiprobability distribution can retain this kind of information, and its knowledge could assist to build up a minimally invasive measurement protocol that is able to return results beyond average values while still characterizing the statistics of energy records.

A. Quantum internal energy distribution
From a stochastic thermodynamic point of view, any internal energy difference of a quantum transformation is a stochastic process.This holds even for isolated systems, since energy-change fluctuations-that however average to zero-are induced by the measurement apparatus.We recall that the latter irreversibly perturbs the measured (thermodynamic) system in any procedure se-quential projective measurements that are directly performed on the system.
In the following, we are going to introduce the concepts of stochastic internal energy and stochastic quantum work in a generic quantum scenario with arbitrary density operators and time-dependent thermodynamic transformations.Let us identify the time-dependent quantum observables O 1 (t 1 ) and O 2 (t 2 ) with the Hamiltonian operators H(t 1 ) and H(t 2 ) at the initial and final times of the transformation under scrutiny.The Hamiltonian operators admit the spectral decomposition: where i, f denote the indexes on the initial and final energies, respectively.From Eqs. ( 91)-( 92), the definition of the stochastic internal energy ∆U if follows.∆U if is defined within the time interval [t 1 , t 2 ], and it is given by the differences ∆U if ≡ E f − E i .Notice that ∆U if depends only on the eigenvalues of the Hamiltonian H at the initial and final times of the thermodynamic transformation described by the CPTP map Φ which the system is subject to, and not directly on Φ itself.On the other hand, what is dependent on Φ is the probability distribution ruling the occurrence of any value of ∆U if .In order to describe such a distribution, we introduce the Kirkwood-Dirac quasiprobability q if ≡ q(E i , E f ) defined as with ρ denoting the initial density operator, such that the QD P [•] of ∆U if is Notice that, here, we have adopted a simplified notation for the KDQ q if using as subscript the labels for the initial and final energies, respectively.Following what we previously stated in Sec.II C about the properties of a KDQ, all the information about the statistics of the stochastic internal energy ∆U if is also contained in the characteristic function obtained by the Fourier transform of P [∆U if ].As such, the KDQ distribution of the internal energy variation can be directly evaluated by means of the interferometric scheme discussed in Sec.III B. As in the general case, the characteristic function G U (u)-as well as each KDQ q ifis formally a quantum correlation function that takes the form For the case of time-dependent unitary dynamics (possibly leading to work fluctuations), where H H (t 2 ) = U † H(t 2 )U is the evolution of the system's Hamiltonian at the final time of the work protocol in the Heisenberg picture.
It is worth stressing that, when [ρ, Π i (t 1 )] ̸ = 0 or [Π i (t 1 ), Π H f (t 2 )] ̸ = 0, the corresponding KDQ q if can be a complex number, with possibly a negative real part.We recall that we have denoted this circumstance as being nonclassical.Even in this quantum thermodynamics case, the nonclassicality of the internal energy distribution P [∆U if ], in the time interval [t 1 , t 2 ], can be measured via the functional ℵ computed over the KDQ q if .

B. Quantum work & KDQ-correction to Jarzynski equality
In any closed quantum system that is driven by a timedependent Hamiltonian H in the time interval [t 1 , t 2 ], the internal energy difference corresponds to the exerted work W .This means that ∆U = W , being the dissipated heat equal to zero in such a case.
If the initial state of the system is a an equilibrium thermal state at inverse temperature β: where Z(t k ) ≡ Tr[e −βH(t k ) ] is the system's partition function, then the TPM probability distribution of work fulfils the celebrated Jarzynski equality (JE) [37,105] ⟨e Eq. ( 99) relates a fluctuating physical quantity (the work W ) measured for an out-of-equilibrium system in a given time during the work protocol, and the equilibrium freeenergy difference The equilibrium free-energy difference is achieved asymptotically by the driven quantum system under the assumptions that, once the work protocol is over, (i) the Hamiltonian of the system is assumed constant and equal to H(t 2 ); (ii) the system is put in contact with a thermal bath at inverse temperature β.Moreover, in Eqs. ( 99)-(100), it is implicitly assumed that the quantum system is connected to the thermal bath at inverse temperature β also before the work protocol is applied.
As surveyed in Ref. [106][107][108][109], the symmetries allowing for the JE in Eq. ( 99) to hold are generally maintained as long as (I) the initial density operator is thermal at inverse temperature β, namely ρ = ρ th (t 1 ) ≡ e −βH(t1) /Z(t 1 ); (II) the dynamics of the quantum system is unital, meaning that the identity I is a fixed point of the quantum map to which the system is subject: Φ[I] = I.Notice that unitary dynamics are a subgroup of such more general family of maps.Therefore, a requirement for the validity of the JE is that ρ is a thermal state, i.e., both (a) [ρ, H(t 1 )] = 0, and (b) the diagonal elements of ρ (with respect to the eigenbasis of H(t 1 )) follow a Boltzmann distribution.
As a result, the JE in Eq. ( 99) can be obtained by applying the TPM scheme, which returns the following characteristic function for any given work distribution: ] the joint probabilities of the TPM scheme, and u complex number.Hence, by setting u = iβ, we get Moreover, if one applies the Jensen inequality on both sides of Eq. ( 102), we directly get the inequality that is one of the formulation of the second law of thermodynamics in relation with the Clausius theorem.
Let us now start connecting these results with the discussion undertaken in the previous sections.In this regard, we already know from Sec. II that, if the density operator ρ at the beginning of the work protocol does not contain quantum coherence χ with respect to the eigenbasis of H(t 1 ) (ρ = D 1 [ρ]), then the first energy measurement of the TPM scheme is not invasive.Hence, in such a case, ⟨W ⟩ TPM − ∆F denotes, without any ambiguities, the dissipated work that is the amount of work that cannot be converted in extracted work.
However, the JE breaks down if ρ is not a thermal state.The failure of the JE also occurs when ρ is thermal but the dynamics of the quantum system is nonunital, possibly leading to heat dissipation [110][111][112][113].As a consequence, one ends up with an expression similar to the JE that exhibits a correction that is not a state function and depends on both the dynamical map to which the quantum system is subject and its initial state.We stress that, in the general case, a correction is present also in the case no quantum coherence χ is contained in ρ, i.e., ρ = D 1 [ρ].Accordingly, by setting ρ = D 1 [ρ], the characteristic function of the work distribution provided by the TPM scheme reads as [114][115][116][117][118][119] ⟨e −βW ⟩ TPM = e −β∆F γ , where with ρ th (t 2 ) ≡ e −βH(t2) /Z(t 2 ).Being expressed as a function of the dynamical transformation applied to the system, the efficacy γ ≥ 0 depends on the time t 2 , making γ generally a time-dependent quantity.Of course, γ = 1 ∀t 1 , t 2 if ρ = ρ th (t 1 ) and Φ is a unital map.We recall that the difference D 1 [ρ] − ρ th (t 1 ) is commonly known as athermality as it quantifies the nonthermal contributions in the diagonal of ρ with respect to the eigenbasis of H(t 1 ).The athermality can be a significant thermodynamic resource if the quantum system undergoes dynamics with feedback [118,120,121].
If we now abandon the use of the TPM scheme and we consider in the more general framework of a quasiprobability distribution, how is the average exponentiated work ⟨e −βW ⟩ modified when [ρ, H(t 1 )] ̸ = 0, namely when quantum coherences are present in the initial density operator ρ?It is indeed clear that, if ρ is not thermal, the JE in Eq. ( 99) is no longer valid and a further correction has to be included to attain a modified JE expression.Notice that a different correction has to be considered for all the protocols that go beyond the TPM scheme [85,[122][123][124].
The use of KDQ to describe quantum work fluctuations leads to the relation where is the KDQ-correction to the JE that holds for any CPTP map Φ.In conformity with the results shown in Sec.II, Γ = γ when ρ = D 1 [ρ], i.e., under the commutative condition [ρ, H(t 1 )] = 0. Similar to the efficacy γ, the KDQ-correction Γ to the JE is not a state function, and therefore depends on the specific thermodynamic transformation that is performed on the system.However, in contrast to the TPM result, Γ is in general a complex number, whose real part can take both negative and positive values.Consequently, as a possible application, if one measured Re [Γ] < 0 or Im [Γ] ̸ = 0, it would imply the presence of nonclassicality, since the nonpositivity functional ℵ in Eq. ( 25) would necessarily be greater than zero.The thermodynamic meaning of the KDQcorrection Γ is still lacking, and further investigations are thus needed.

C. nonclassical work exerted by qubits: a case study
In this section, we discuss a simple example to analyze the KDQ distribution of work done on a single qubit that is driven by a work protocol described by a unitary operator U .Assuming the system does not interact with any external bath, the internal energy change can be fully identified as work.Albeit simple, this model can be solved analytically and finds experimental applications in nuclear magnetic resonance (NMR) spin systems [124] and nitrogen-vacancy (NV) centers [12,86] (point defects in the diamond lattice), where experiments of quantum thermodynamics beyond TPM have been recently performed.
Let us assume the Hamiltonian of the qubit to be corresponding to a spin-1/2 particle subject to an effective magnetic field rotating around the z-axis.In the rotating frame, described by the unitary operator U rot = e iδσ z t/2 , the effective Hamiltonian governing the dynamics of the qubit becomes time-independent and reads so that the system's evolution operator (in the original frame) is To find the statistics of work done between times t 1 = 0 and t 2 = t, we use the spectral decomposition of the time-dependent Hamiltonian, i.e., where we have defined a generalized Rabi frequency ∆ ≡ √ δ 2 + Ω 2 .Moreover, we assume the system to have quantum coherence in the eigenbasis of H(0), so that where 0 ≤ p ≤ 1 corresponds to the populations of the initial eigenstates, and c is the quantum coherence that we have chosen to be real for simplicity.Using the definition of the quasiprobability distribution for work, see Eq. ( 94), we find: where the KDQ q are [see Eq. ( 93)]: First, we notice that the imaginary parts of q if are always proportional to the coherence c.Second, the real parts of q if may become negative.To see this, we specify, for simplicity, initial conditions and take the maximum possible coherence: p = c = 1/2.In this case, the real parts Re[q if ] become: Re It is possible to see that the minimum value of Re The time-dependence of the quasiprobabilities q if is shown in Fig. 3 for the two cases Ω = ( √ 2 ± 1)δ.It is interesting to see that only one of the Re[q if ] may become negative for each value of Ω.Moreover, choosing c complex may lead to another of the Re[q if ] to become negative, but the minimum value is still (1 − √ 2)/4.

D. Enhancement of extractable work
In this section, we are going to explain the meaning of nonclassical work extraction and anomalous energy exchange or variation.
In any system that is subject to a work protocol, the extractable work is defined as the amount of energy that is left over, with respect to the energy of the system at the beginning of the transformation.Accordingly, if a protocol admits non zero extractable work, then the average energy at the end of the work protocol, ⟨H(t 2 )⟩ = Tr[U ρ U † H(t 2 )], is smaller than the average energy at the beginning, ⟨H(t 1 )⟩ = Tr[ρH(t 1 )], so that the extra energy amount can be used by a work reservoir [125] or stored in a battery [126].Hence, the requirement for work extraction is that Recently, it has been discussed whether the negativity of the terms composing a quasiprobability work distribution may correspond to an enhancement of work extraction, and whether this circumstance can be witnessed by violating an inequality that is instead valid under the commutative conditions [ρ, H(t 1 )] = 0 and [H(t 1 ), H(t 2 )] = 0, i.e., when ℵ = 0.The answer to both these questions is positive [12].
In order to see this, at the level of energy transitions, let us consider the fact that an excitation process ∆U if ≡ E f − E i > 0 (indexes i, f label the initial and final energies, respectively) occurring in a quantum process with negative quasiprobability q if (not neccessarily KDQ) is equivalent to a de-excitation process ∆U if < 0 in a classical work transformation with probability |q if |.
During an excitation (stochastic) process, the system absorbs energy and uses this energy to carry out a transition between the energy levels.On the other hand, any de-excitation process that is operated by a thermodynamic transformation contributes to increase the amount of the extractable work.Therefore, the presence of negative quasiprobabilities can be effectively exploited as a resource to enhance work extraction, beyond what any classical stochastic process can achieve.Such enhancement is deemed as 'nonclassical', and the internal energy variations ∆U if associated to negative probabilities, enabling it, are called 'anomalous'.Thus, anomalous energy variations denote energy exchanges that are inherently quantum mechanical, and heralded by nonpositivity.
Let us see how the enhancement of work extraction occurs.If one uses a work protocol (such as the TPM scheme) returning positive joint probabilities p if of the work distribution, the work extraction is maximized if we minimize (with sign) ⟨W ⟩ clas = i,f p if ∆U if , where the subscript 'clas' here specifically stands for 'classical' in the sense of positive joint probabilities.Without specifying anything about the thermodynamic transformation, the necessary condition to achieve the largest extractable work is to set such that leads to extractable work (in absolute value).
If instead the statistics of the internal energy variations ∆U if are described by some quasiprobabilities q if (for example when [ρ, H(t 1 )] ̸ = 0, as shown in Sec.II), then the extractable work can be enhanced beyond what is obtained by a work protocol returning positive joint probabilities.For such a purpose, one can set In this way, the magnitude of the extractable work can be effectively maximized in order to satisfy the inequality To achieve this, any excitation process ∆U if > 0 has to be associated to a negative Re[q if ], while any deexcitation process ∆U if < 0 should occur with a positive quasiprobability.It is worth observing that, in the case a KDQ distribution of work is taken into account, the imaginary parts of KDQ do not play an effective role for the task of work extraction, since they do not affect the average work.What really matters to get enhanced work extraction is to ensure nonclassical behaviours in the time-distribution of negativity, namely the distribution over time of the quasiprobabilities with negative real part.In fact, it is not necessarily crucial for the nonpositivity functional ℵ to take a large value, but that a significant negativity is associated to a positive 'anomalous' energy variations ∆U if > 0. At the same time, work extraction is enhanced when negative values of ∆U if occur with the largest possible positive quasiprobability q if .The interplay of all these conditions is model-dependent and depends on the specific parameters that rule the dynamics of the work process.Therefore, it is evident that attaining enhanced work extraction stems from an optimization routine that makes Eqs.(127) valid in a given time interval of the work protocol.
In Ref. [12], the electronic spin of an NV center in bulk diamond at room temperature was considered as the system to which a work protocol would be applied.Work extraction was observed, and its maximum values were associated to negative MHQ fulfilling Eq. ( 129).The work extraction enhancement observed in Ref. [12] originates from a sub-optimal solution for the optimization of work extraction against the time duration of the work protocol.In fact, due to the experimental constraints, only one internal energy change ∆U, corresponding to the largest possible value, was associated with a negative MHQ.At the same time, the smaller internal energy variation −∆U occurred with positive quasiprobability, with all other MHQ being negligible.

Enhanced extractable work from violating a classical inequality
We are going to show that fulfilling Eq. ( 129) implies the violation of an inequality for work extraction that holds as long as the commutativity condition [ρ, H(t 1 )] = 0 is obeyed; as in [12], we consider MHQ.The violation of such an inequality cannot occur in any experiment realizing a work protocol yielding positive work joint probabilities, as the TPM scheme.
Let us thus consider that the projectors Π i and Π f of the Hamiltonian at the initial and final times t 1 and t 2 of the work protocol are rank-1 operators.This means that Π Moreover, we assume, for simplicity, the initial density operator ρ = |ψ⟩⟨ψ| to be pure.Under these assumptions, the MHQ takes the form Interestingly, all the terms ⟨ψ|E i ⟩, ⟨E i |U † |E f ⟩, ⟨E f |U |ψ⟩ are complex numbers whose real parts are linked with a standard probability amplitude, either defined at a single time or measurable by means of the TPM scheme.In particular, for the probability p i to measure the initial energy of the system, one has where ϕ i is a phase factor.Then, in the same spirit, we can write and where φ if , θ f are the corresponding phase factors.
In Eq. ( 132), p f |i is the conditional probability (associated to the TPM scheme) of measuring the energy E f at time t 2 conditioned to have measured E i at time t 1 .Instead, in Eq. ( 133), p f is the probability to measure the energy E f at the end of the work protocol, by initializing the system in ρ = |ψ⟩⟨ψ|.Notice that, by construction, the probability p f encodes information on the quantum coherence that is initially present in ρ; for this reason, p f is a key element of the EPM scheme [85,86,127].
Overall, combining Eqs. ( 131)-( 133) we arrive at where, by definition, p if = p f |i p i is the joint probability returned by the TPM scheme, and Λ if ≡ cos(ϕ i + φ if + θ f ) that is named activity [12].The latter brings information on the quantum interference fringes among the eigenbasis of ρ, Π i (t 1 ) and Π H f (t 2 ).It is indeed the activity Λ if that is responsible for the negativity of Re If we substitute Eq. ( 134) into the expression of the work extraction in (128), we find whenever Λ if ≥ 0 ∀ i, f .The inequality in Eq. ( 135) gives an upper bound, dependent on the work protocol, to the amount of extractable work when ℵ = 0, and applies also to initial quantum mixed states.Hence, a violation of this bound, as experimentally tested in Ref. [12], is a witness of the presence of negativity, as well as of nonclassical work extraction.In Fig. 4, we show an example of the enhancement of work extraction aided by negativity using the work protocol introduced in IV C, applied to a driven qubit.In √ 2)δ.In particular, we plot the average ⟨W ⟩ of the KDQ distribution of work (blue solid line), the average work ⟨W ⟩TPM returned by the TPM scheme (orange dashed line)-equal to zero for any time-and the classical bound from Eq. ( 135) (green dotted line) taken with opposite sign.particular, we plot the average work of the TPM and MHQ probability distributions using the energies and Hamiltonian projectors in Eqs. ( 112)-( 113), as well as the quasiprobabilities in Eqs. ( 116)- (119), with p = 1/2, c = −1/2 (c = 0 to get the work statistics of the TPM scheme), and δ = Ω/( √ 2 + 1).Interestingly, if the initial state of the qubit [see Eq. ( 114)] is fully mixed (c = 0), then the average work is zero for any value of the final time t 2 , see Fig. 4. On the other hand, turning on the quantum coherence in ρ and making use of quasiprobability to attain the work distribution P [W ], the energy injected by the driving field is transformed into extractable work, beyond the classical bound [right-hand-side of Eq. ( 135)] in the interval (Ωt)/π ∈ [0.6, 1.4] approximately.In this case study with a qubit the classical bound amounts to with U given by Eq. ( 110).

E. Work variance in the KDQ setting
In the previous section, we have shown how using KDQ to describe the work fluctuations makes the average work ⟨W ⟩ = i,f q if (E f − E i ) equal to the corresponding value that is unperturbed by the measurement disturbance.Even though the KDQ q if are complex numbers, the average work ⟨W ⟩ is always a real number, with a clear interpretation with classical physics, as shown above.
In the following, we analyze how the fact that q if are complex numbers affects the second moment of the KDQ distribution of work, P [W ], i.e., the work variance (∆W ) 2 .As noticed in Ref. [128], the variance of work is in general not equal to the variance of the operator H H (t 2 ) − H(t 1 ) calculated for the initial state ρ.Instead, this is formally defined by where, as before, all the averages ⟨•⟩ are performed with respect to ρ, and the second statistical moment ⟨W 2 ⟩ reads as The quantity Tr H H (t 2 )H(t 1 )ρ in the right-hand-side of Eq. ( 137) is a two-time quantum correlation function for the Hamiltonian H and is generally complex, making ⟨W 2 ⟩ also complex.This means that the imaginary part of ⟨W 2 ⟩ is equal to the imaginary part of Tr[H H (t 2 )H(t 1 )ρ], whose meaning lies in the presence of phase correlations in the scalar products of the eigenvectors of ρ and H at the times t 1 , t 2 .For this reason, the quantum correlation function for H preserves information about the quantum coherence contained in ρ, and this feature is transferred to the work variance (∆W ) 2 .In this regard, we are going to show that the imaginary part of (∆W ) 2 , Im (∆W ) 2 , is directly linked with the noncommutativity between ρ and H. From this, following the time-energy Schrödinger-Robertson uncertainty relation [129][130][131], Im (∆W ) 2 is bounded by the product of the uncertainties of H(t 1 ) and H H (t 2 ) respectively.
A first expression of the work variance is obtained by combining Eqs. ( 136)- (137), so that: where and parts of Tr H H (t 2 ) − ⟨H H (t 2 )⟩ H(t 1 ) − ⟨H(t 1 )⟩ ρ are equal respectively to [132]: and Eq. ( 141) defines the quantum covariance of H(t 1 ) and H H (t 2 ).Instead, in Eq. ( 142), Tr[ρ [H H (t 2 ), H(t 1 )]] is a purely imaginary number and, by definition, is the expectation value of the commutator [H H (t 2 ), H(t 1 )] with respect to the initial density operator ρ.
This derivation demonstrates that the work variance has both a real and an imaginary part.The real part has a clear correspondence with the thermodynamics of classical systems, as In addition, the fact that the commutator [ρ, H(t 1 )] ̸ = 0 or [H(t 1 ), H H (t 2 )] ̸ = 0 may lead to a decreased work variance, namely to Re (∆W ) 2 ≤ (∆W TPM ) 2 .We show this for the driven qubit of Sec.IV C and report the results in Fig. 5 where we assume the same values used for Fig. 4. Interestingly, apart from Ωt/π = 0, 2 where both the work average and variances are zero, the real part of the KDQ work variance Re (∆W ) 2 has a local minimum at Ωt/π = 1 that is the time instant with maximum negativity.
On the other hand, the imaginary part of the work variance is that quantifies the noncommutativity of H(t 1 ) and H H (t 2 ).The magnitude of Im (∆W ) 2 can be bounded from above by making use of the time-energy Schrödinger-Robertson uncertainty relation [129][130][131].In fact, the latter states that, for any quantum observables O 1 , O 2 and density operator ρ, where with Therefore, by applying the inequality in Eq. ( 145) to our setting, we obtain F. Heat fluctuations in the quantum regime In this section, we no longer deal with work distributions, and we will focus on heat fluctuations.For this purpose, we consider the paradigmatic model that consists in placing into contact a cold and a hot quantum system, which globally undergo a unitary quantum dynamics.Depending on the initial quantum states of the cold and hot systems, different results as well as thermodynamic interpretations can be drawn.In this context, a first relevant result is in Ref. [133] and goes under the name of Jarzynski-Wójcik exchange fluctuation theorem.In Ref. [133], two quantum systems B c and B h with finite Hilbert-space dimension are prepared in two equilibrium thermal states at different temperatures β c and β h with β c > β h .Then, they are made weakly interacting with one another for a given time interval.Under this assumption, one gets that where Q is the stochastic heat exchanged by the two bodies, and ∆β = β c − β h denotes the difference of the inverse temperatures of the initial thermal states for the two bodies.
If the initial global state of the two systems is a product state then the average ⟨•⟩ in Eq. ( 148) can be performed also with respect to the TPM distribution of the exchanged heat.To find this, it is sufficient to measure the sum of the local Hamiltonian operators of the two bodies, i.e., H = H Bc + H B h (a time-independent Hermitian operator).Furthermore, throughout this section, we also implicitly assume the energy-preserving condition for the unitary operator U that describes the quantum dynamics of the two bodies: Equation ( 149) physically entails that, at any time t, the average energy variation in a body is minus the corresponding average energy variation in the other body.
Such symmetry allows one to study fluctuations of exchanged energy between the two bodies by just measuring one of them.
In the literature, it has been considered also the cases of an initial quantum state that is locally thermal as in Ref. [133], or classically correlated [134].This kind of correlation makes nonthermal the diagonal of the initial density operator ρ for the two bodies taken individually, but does not add off-diagonal elements in ρ with respect to H.As shown in Ref. [134], a generalized exchange fluctuation relation, extending Eq. ( 148), can be still obtained, as we will discuss next.
Let us now introduce the spectral decomposition of the local Hamiltonians H Bc and H B h for each of the two bodies: with k = c, h and ℓ = i, f .This implies that the projectors of the total Hamiltonian H are Π ici h = Π ic ⊗ Π i h .
For the initial state ρ, we require that the reduced states of the each body is in equilibrium at inverse temperature β k : where Z k ≡ Tr e −β k H B k are the local partition functions.We hence have: ρ = ρ th,Bc ⊗ ρ th,B h .While the reduced states are diagonal in the eigenbasis of H B k , in general the global state ρ may contain off-diagonal elements, with respect to the local energy eigenbasis, that may be the signature of the presence of quantum correlations.
We are now in the position to define the average heat flow that, due to the energy-preserving condition, can be inferred from the energy change of either the cold or the hot body.Without loss of generality, we choose to measure it through the cold system as in Ref. [8].The average heat flow at the final time t 2 of the thermodynamic transformation is where ρ ′ = U ρU † denotes the evolved density operator of the two bodies.According to Eq. ( 153), ⟨Q⟩ ≤ 0 denotes heat flowing on average from the hot to the cold body, as naturally requested by the second law of thermodynamics with the intervention of no external drive.On the other hand, by resorting to additional resources, it can also occur that ⟨Q⟩ > 0 meaning that on average heat flows from the cold body to the hot one, as in a refrigerator.Summarising, ⟨Q⟩ ≤ 0 =⇒ hot-to-cold heat flow, ⟨Q⟩ > 0 =⇒ cold-to-hot heat backflow.
Moreover, if the amount of exchanged heat from the cold to the hot body exceeds in magnitude the value of (∆β) −1 ln d, with d the Hilbert-space's dimension of each body, then the heat backflows are called strong.Notably, the observation of strong backflows indicates that the actual quantum state of the bipartite system, on which heat fluctuations are evaluated, is entangled [135].
Let us introduce the quasiprobabilities q ici h fcf h associated to the energy variations ∆E ici h fcf h , eigenvalues of the total Hamiltonian H of the two bodies.Using again the definition of the QD in Eq. ( 94), one has that where Π H fcf h = U † Π fcf h U and is the corresponding joint probabilities returned by the TPM scheme.As before, in Eq. ( 155) we have employed the dephasing operator 154)-( 155) the initial quantum state ρ of the two bodies is linearly decomposed as ρ = D 1 [ρ] + χ [see Eq. ( 5)], where both the diagonal and off-diagonal parts of ρ are considered with respect to the eigenbasis of H.
Based on this framework, we describe an exchange fluctuation relation that is also valid in the noncommutative regime of [ρ, H] ̸ = 0, due to the presence of quantum correlations or entanglement in the initial state.For this purpose, let us introduce the stochastic mutual information I with elements where j = i, f , such that ∆I ≡ I fcf h − I ici h .We also recall that the energy variation in the cold body is assuming the energy-preserving condition for U and a resonant interactions between the bodies.Hence, we find [8] ⟨e ∆I+∆βQ ⟩ = 1 + Υ, (157) where the average ⟨•⟩ in Eq. ( 157) is made with respect to the quasiprobabilities q ici h fcf h , and The correction to the exchange fluctuation theorem, Υ, is equal to zero if [ρ, H] = 0, which is equivalent to ρ = D 1 [ρ] and χ = 0.In this case the application of the TPM scheme suffices.The exchange fluctuation relation of Eq. ( 157) reduces to the Jarzynski-Wójcik identity in Eq. ( 148) in the case ρ = ρ th,Bc ⊗ ρ th,B h , whereby ∆I = 0. Instead, if the diagonal elements of ρ are not thermally distributed-due to classical correlations in the H eigenbasis-and χ = 0, then one recovers the exchange fluctuation relation in Ref. [134], i.e., We conclude this theoretical analysis about heat fluctuations in the quantum regime, by providing the thermodynamic interpretation of the fluctuation profiles associated to the quasiprobability distribution of heat exchanges.Previously, we have seen that the presence of quantum correlations in the initial state can enhance the amount of heat backflows, such that ⟨Q⟩ ≥ 0 according to the used convention.The explanation of this phenomenon lies in the possibility to associate negative quasiprobabilities Re [q ici h fcf h ] to positive heat exchanges Q = E ic −E fc > 0 corresponding to heat flowing from the cold body to the hot one.Such a process, which needs an external energy source for its activation, is triggered by quantum correlations.
Notice that, in order for quantum correlations to be effectively considered as a resource for heat backflows, it is required that negative heat exchanges Q ≤ 0 (i.e., energy fluxes from the hot to the cold bodies) occur with positive quasiprobabilities Re [q ici h fcf h ], similarly to what happens to work extraction in Sec.IV D. When the enhancement induced by quantum correlations in ρ allows for strong cold-hot heat backflows, then one can state that the corresponding energy exchange process takes nonclassical traits.

Example: two-qubit system
We now apply the theoretical framework introduced above to a pair of interacting qubits, at inverse temperatures β c and β h , respectively, and local Hamiltonians H B k = Ωσ z k , k = c, f .The two qubits are initialized in a global state ρ containing off-diagonal elements with respect to H = H Bc + H B h .
As proven in [8], a general form of the initial state for a two-qubit system fulfilling the requirements ( 151)-( 152) is where p ∈ [0, 1] is a population term, [8], it is also shown that the energy-preserving condition [Eq. (149)] is responsible to set a minimal form for the unitary operator U : with θ ∈ [0, 2π], equivalent to a partial swap transformation.Under these assumptions, the analytical expression of the heat exchange (153) between the two bodies reads as [8] ⟨Q⟩ = −η cos(ξ) sin(2θ) + ⟨Q⟩ TPM , where is the average heat flow obtained by applying the TPM scheme, or by setting η = 0 (no quantum coherence) in Eq. (162).For any value of θ, β c and β h , it can be easily found that ⟨Q⟩ TPM ≤ 0, which means that no cold-to-hot heat backflow is possible.This is evident in Fig. 6 where we plot Eq. ( 162) against θ for η = 0 and for η ̸ = 0 for fixed p, ξ, β c , β h .From the figure it can be observed that, for some values of θ, ⟨Q⟩ > 0 (cold-to-hot heat backflows) is possible when η ̸ = 0.The parameter η also affects the magnitude of the heat exchanged between the cold and hot systems.

V. QUANTUM MANY-BODY SYSTEMS
Quasiprobability distributions play an instrumental role in the understanding of quantum many-body systems.In this section, we will present their applications in the context of quantum scrambling and out-oftime-ordered correlators (Sec.V A), the Loschmidt echo (Sec.V B) and a technique for efficiently calculating the KDQ distribution of free fermion systems (Sec.V C).
In this section, we review the basic definition of OTOCs and show how they relate to quasiprobabilities.This connection has been recently noted in various works [4,28,59,150].
Following Ref. [150], let us consider a system that is initially prepared at time t = 0 in the pure state |ψ⟩ and initially perturbed by a unitary operator V , acting locally on a part of the system.For instance, if the system is made of qubits, we may consider the spin flip V = σ x j that acts on qubit j.The system is then evolved for a time interval t following a time evolution described by the unitary operator U .At the end of the dynamics, the system is perturbed by the application of another unitary operator Y [151], acting locally on another disjoint part of the system.Finally, the system is evolved backward in time through the operator U † .At the end of this protocol, the state of the system is: Now, suppose that we perform an alternative protocol in which the perturbation V is applied, not after the initial preparation of ρ, but after the backward evolution.This results in the state |ψ ′′ ⟩ = V U † Y U |ψ⟩.The overlap between these two states equals the OTOC where we have defined Y t = U † Y U that denotes the perturbation operator Y evolved in the Heisenberg picture.The definition in Eq. ( 164) can be extended to initial mixed states ρ and from now on we will assume: While the OTOC is in general a complex number, one can consider a real quantity by introducing the OTO commutator: Its interpretation is the following.Initially, at time t = 0 (U = I), the operators Y t = Y and V commute as they have support on spatially separated parts of the system: [Y, V ] = 0.However, as time progresses, the effects of the perturbation Y t may reach the support region of V and their commutator might become nonzero: [Y t , V ] ̸ = 0.The quantity C(t) measures the magnitude of this commutator.If both operators V and Y are unitary, the OTO commutator is related to the OTOC by the relation This shows that the growth of the commutator C(t) is associated with a decay of the real part of F (t). Next, we are going to prove that an OTOC is equal to the characteristic function of a KDQ.To this end, let us follow Ref. [152], which introduces the wing-flap protocol and express the unitary operator V in exponential form as V (u) = e iuO , with O a Hermitian operator and u a real scalar.The spectral decomposition of the observable O reads: O = m o m Π m , in terms of its real eigenvalues o m and the corresponding orthogonal projectors Π m .Using these definitions, V can be expressed as The wing-flap protocol consists of the following steps: 1. Prepare the system in the state ρ.
3. Evolve the system forward in time with U .
4. Apply the perturbation Y .
5. Evolve the system backward in time with U † .
Recalling the definition in Eq. ( 16), we can define the KDQ for the change ∆o nm = o m − o n in the eigenvalues of O when two measurements of O are performed at times t 1 = 0 and t 2 = t, respectively (steps 2 and 6 of the wing-flap protocol).In Eq. ( 168), the operator Y t plays the role of the complete evolution operator Y t = U † Y U that combines the steps 3-5 of the protocol between the two measurements of O. Notice that, in order to denote a quasiprobability, we continue using the simplified notation adopted in Sec.IV, whereby the subscript in q nm contains the indexes labelling the measurement outcomes at the initial and final times of a two-time procedure.As a result, the quasiprobability distribution to observe a change ∆o(t) at time t is given by and the characteristic function of P [∆o, t] is its Fourier transform (see Sec. II C 3), such that that can thus be expressed as an OTOC.
Similar to Secs.II-IV, we can define the corresponding MHQ distribution that can be associated with an OTOC.Such a distribution is the real part of the corresponding KDQ distribution, that is r nm = Re [q nm (t)], and its characteristic function reads as Interestingly, using the equality we obtain Let us now consider a practical example and let In Fig. 7(a)-(b), we show the real and imaginary parts of q nm .The quasiprobability q 00 is the only one whose real part becomes negative.Since the probabilities need to fulfil the normalization condition, it means that q 11 is amplified due to the presence of quantum coherences in the initial state ρ, compared to a case in which the initial state commutes with the observable O.
In Fig. 7(c), the nonpositivity functional ℵ is plotted against the time for different values of J.In Fig. 7(d), we show a contour plot of the minimum value of ℵ in the time interval 0 ≤ t ≤ 20 as a function of β and J. First, as expected, we confirm that larger values of the interaction strength J always lead to a stronger nonpositivity.Second, the nonclassicality reduces as the temperature increases.

B. Link with the Loschmidt echo
Another interesting connection of quasiprobability distributions with the irreversibility of many-body systems arises in the context of quantum chaos and decoherence.Consider a system of many particles, classical or quantum, that evolves in time for a period t according to a time-independent Hamiltonian H 0 .If we were to invert all momenta and run the evolution backward we should be able to recover the initial state.However, little imperfections in the inverted evolution or decoherence induced by an external environment may cause some deviations.For an initial pure state |ψ⟩, one can define the Loschmidt echo (LE) L(t) = |G(t)| 2 as the absolute square value of the complex amplitude whose generalization to mixed states and non unitary evolutions is straightforward.Mathematically, L(t) represents the fidelity, in terms of the overlap, between the initial state |ψ⟩ evolved with the unperturbed Hamiltonian H 0 and the state |ψ⟩ evolved with the perturbed Hamiltonian H δ .Peres transferred the LE idea in the quantum domain [99], while Ref. [153] used the LE to analyze the decoherence of a many-body spin system and the relation to chaos.For the quantum version of systems with a classically chaotic Hamiltonian (for instance a particle moving in a driven double well, see Ref. [153]) the rate at which the information about the initial state is destroyed by the environment, within a range of couplings to the environment, is set by the classical maximal Lyapunov exponent.Under these assumptions, the LE decays exponentially in time with the Lyapunov exponent, thus revealing the underlying classical chaotic behaviour; see also Ref. [100].
Moreover, the LE was beneficial to uncover a new type of phase transition occurring in time.In this regard, if we assume that the initial state |ψ⟩ is the ground state of H 0 with zero energy, then the LE amplitude (176) reduces to which looks like the partition function of the Hamiltonian H δ but with an imaginary inverse temperature it.Since classical phase transitions arise because of singularities in the partition function, Heyl and coworkers discovered dynamical quantum phase transitions as those that give rise to singularities in the LE at specific instants of time, see Refs.[154,155].
The LE is also strongly connected with the statistics of work as mentioned in Sec.IV and described in detail in Ref. [156].In fact, Eq. ( 176) can be interpreted as the characteristic function of the work done on a quantum system initially in the state |ψ⟩, whose Hamiltonian is subject to a quench dynamics that instantaneously changes H 0 to H δ .
Let us now formalize the connection between the LE and the KDQ.First, we write the spectral decomposition of the two Hamiltonian operators: m .With these assumptions, let us write an expression for the LE for a generic mixed initial state: where we have introduced the KDQ q nm of the random variable W = E m − E n , defined as Thus, the inverse Fourier transform of G(t) with respect to time t, i.e., can be interpreted as the quasiprobability distribution for the work m − E n done on the quantum system, which initially is in the state ρ and whose Hamiltonian is suddenly changed from H 0 to H δ .It is worth noting that in contrast to the case of the characteristic function of work obtained using the TPM scheme [156], in the general case the initial state ρ may not commute with any of the two Hamiltonian operators H 0 and H δ .This fact, as we have seen in other examples earlier, may give rise to a distribution P [W ] with nonpositive values.
We now proceed to illustrate these concepts with a simple example.Let us consider a qubit in the pure initial state and let us choose the Hamiltonian operators H 0 and H δ as The eigenstates of H 0 are simply |0⟩ and |1⟩, with eigenvalues ±B respectively, and we define the projectors on these states as Π i = |i⟩⟨i| with i = 0, 1.Similarly, for H δ , the eigenstates are where the possible values of the quasimomenta are k = 2πm/N with m = −N/2 + 1, . . ., N/2 (assuming for simplicity N even).Let us thus define the fermionic operators γ k that are obtained by applying the following Bogoliubov rotation to the operators c k : The fermionic operators γ k depend on the angles θ k , implicitly given by and satisfy the canonical anticommutation relations In terms of the operators γ k , the Hamiltonian Eq. ( 189) reduces to a diagonal form: where are the single-particle eigenenergies.
In what follows, we consider a sudden change of the Hamiltonian H(λ) in which the magnetic field λ is changed instantaneously from the initial value λ 0 to the final value λ 1 .To calculate the quasiprobabilities defined in Eq. ( 93), we need to express the fermionic operators γ (1) k , which diagonalize H(λ 1 ) with eigenenergies ϵ (1) k = ϵ k (λ 1 ), in terms of the fermionic operators γ (0) k diagonalizing H(λ 0 ) with eigenenergies ϵ (0) k = ϵ k (λ 0 ).This is possible thanks to the linear Bogoliubov transformation [158]: where ∆ k ≡ θ Here, the need for a quasiprobability approach arises whenever one chooses an initial state ρ that has coherences in the eigenbasis of H(λ 0 ).In our example, we choose the following state: Let us now calculate the KDQ distribution of the work done by suddenly change the Hamiltonian from H(λ 0 ) to H(λ 1 ).Since ϵ k = ϵ −k , we can rewrite the initial state as: From Eq. (198), we see that the eigenstates of H(λ 0 ) with quasimomenta ±k are transformed into the superposition of eigenstates of the H(λ 1 ) with the same pair of quasimomenta.Therefore, we can compute the work done for all possible transitions from |m k , n −k ⟩ to m ′ k , n ′ −k that correspond to the work instances W mn,m ′ n ′ (k).The only transitions that have nonzero quasiprobabilities k + ϵ (0) k q 00,00 (k) = e βϵ (0) k − ϵ (0) k q 11,00 (k) = e −βϵ (0) where, for each process, we have included the value of the stochastic work instances and the associated quasiprobability.As expected, the quasiprobabilities with non vanishing work also depends on the mixture parameter p, which weighs the contribution of the initial state ρ containing quantum coherence in the eigenbasis of H(λ 0 ).menta k > 0: A coarse-grained histogram of P [W ] is shown in Fig. 9.For p = 0 the distribution is always nonnegative, while for p = 1 some negative parts appear and are associated with positive values of W .As a consequence P [W < 0], leading to work extraction, tends to be enhanced, so that ⟨W ⟩ < 0, as explained in Sec.IV D and shown explicitly in Fig. 10.Notice that, in order for P [W ] to exhibit negativity, the temperature entering ρ G (λ 0 ) must be high enough for the two-body processes 00 ↔ 11 to be significant.In contrast, if the initial state is close to the ground state these processes are suppressed and P [W ] is non negative everywhere.Initial state coherence also leads to a reduction of the work fluctuations as measured by its variance, see Fig. 10.

VI. DISCUSSION
Quasiprobabilities have been quite elusive quantities so far, due to the difficulty for their experimental inference.As stressed in Sec.II, procedures based on sequential projective measurements cannot reconstruct the quasiprobability distribution of a physical quantity that is defined over two times, as well as its corresponding statistical moments.
Recently, however, we have witnessed a resurgence of quasiprobabilities, thanks to their direct link with two-point quantum correlation functions of the form ⟨O 1 (t 1 )O 2 (t 2 )⟩, with O 1 (t 1 ), O 2 (t 2 ) quantum observables, and the average ⟨•⟩ performed with respect to a generic density operator ρ.Quantum correlation functions are a powerful tool to describe phase changes in quantum statistical mechanics.Hence, the possibility to express them in terms of quasiprobabilities opens the door for building a microscopic, nonequilibrium description of phenomena that naturally includes genuinely quantum resources as quantum coherence and correlations.
Beyond theoretical arguments, quasiprobabilities may turn out to be pivotal also for revealing advantages in quantum technology applications.In this tutorial we have seen several examples in quantum thermodynamic applications, particularly in experiments, including work extraction [12].This is also relevant for the energetic assessment of quantum computation, where the energy exchange of a qubit with its environment can be continuously monitored through weak measurements [25].
The contextual nature of quantum systems facilitates the emergence of anomalous weak values of the energy exchanges due to quantum coherence which manifest themselves through negative quasiprobability distributions.In the context of metrological applications, negative values of quasiprobability distributions may enhance parameter estimation increasing the precision of quantum sensing protocols [7,19,29].
We conclude the tutorial by mentioning some possible future perspectives of the topics treated here.By now, it should be apparent how quasiprobabilities have connections with all main theoretical and experimental aspects of quantum theory.Here, we explored the direct link with quantum measurement theory, fluctuation theorems, work and heat in quantum systems led by genuinely quantum resources, and the scrambling of information in many-body systems.Thus, further investigations on the following subjects could be considered: (i) To determine how the thermodynamic entropy production in a nonequilibrium quantum process is expressed in terms of a quasiprobability distribution.A starting point could be taking Ref. [160], where the entropy production as a quantifier of irreversibility is extended to a regime with noncommuting conserved quantities.Afterwards, one might investigate the link with quantum information theory and feedback mechanism naturally including the so-called Maxwell's demon [161].In this regard, a quasiprobability formulation of quantum trajectories [162] could be taken into account.
(ii) The extension of two-time quasiprobability distribution to access multitime statistics in open quantum systems (see for instance Refs.[69,74] for similar attempts for other QDs).This could help investigating non-Markovianity arising because of memory effects in the environment.
(iii) To define to what extent the quasiprobability distribution underlying an OTOC can be a proper quantum sensing toolbox.In fact, given a quantum many-body system, different perturbations may scramble differently the state of the global system [149], and the corresponding quasiprobability distribution could give access to this information, measurable by means of an interferometric procedure.
We hope that the curious and interested reader can find new, fascinating ideas from this tutorial, and can develop some of the perspectives listed here, by opening in turn further open problems.

DATA AVAILABILITY STATEMENT
The codes implementing the computations for the figures of the tutorial are available as Supplemental Material [163].

FIG. 4 .
FIG.4.Average work in units of Ω for the spin-1/2 model described in Sec.IV C as a function of the final protocol time t2 = t.We consider the parameters p = 1/2, c = −1/2 and Ω = (1 + √ 2)δ.In particular, we plot the average ⟨W ⟩ of the KDQ distribution of work (blue solid line), the average work ⟨W ⟩TPM returned by the TPM scheme (orange dashed line)-equal to zero for any time-and the classical bound from Eq. (135) (green dotted line) taken with opposite sign.

FIG. 5 .
FIG. 5. Plot of the work variances of the KDQ (real part, blue solid line) and TPM (orange dashed line) distributions, respectively.The parameters are the same as in Fig. 4.

1 and O = σ z 2 .
be the Hamiltonian of two qubits initially in the thermal state ρ = e −βH Tr [e −βH ] (175) with inverse temperature β.Then, we choose Y = σ z Notice that neither Y nor O commute with either H or ρ.Therefore, this is an ideal setting to test the possible presence of nonpositivity in the KDQ (168).For this purpose, we write the measurement observable as O = |0⟩⟨0|−|1⟩⟨1|, from which Π 0 = |0⟩⟨0| and Π 1 = |1⟩⟨1| with the corresponding eigenvalues o 0 = 1 and o 1 = −1.We stress that the projectors Π 0 , Π 1 of O act locally on qubit 2. One can see that the evolution operator U = exp(−iHt) is periodic with the frequency ω

FromFIG. 10 .
FIG. 10.Average work (absolute value) and work variance for the quantum Ising model as a function of the weight p introducing quantum coherence in ρ.The other parameters are the same as in Fig. 9.