---
abstract: 'Having accurate tools to describe non-classical, non-Gaussian environmental fluctuations is crucial for designing effective quantum control protocols. We show how the Keldysh approach to quantum noise characterization can be usefully employed to characterize frequency-dependent noise, focusing on the quantum bispectrum (i.e. frequency-resolved third cumulant). Using the paradigmatic example of photon shot noise fluctuations in a driven bosonic mode, we show that the quantum bispectrum can be a powerful tool for revealing distinctive non-classical noise properties, including an effective breaking of detailed balance by quantum fluctuations. The Keldysh-ordered quantum bispectrum can be directly accessed using existing noise spectroscopy protocols.'
author:
- 'Yu-Xin Wang'
- 'A. A. Clerk'
bibliography:
- 'refFCS.bib'
title: 'Spectral characterization of non-Gaussian quantum noise: Keldysh approach and application to photon shot noise'
---

[*Introduction–* ]{}An accurate description of environmental fluctuations is crucial for quantum information processing and quantum control. While it is common to assume noise that is both classical and Gaussian, there are many physically-relevant situations where these assumptions fail [@Cywinski2014; @Viola2014; @Viola2016; @Oliver2019; @Ramon2019]. Understanding how to usefully characterize non-Gaussian, non-classical noise in a frequency-resolved manner could enable the design of more optimal dynamical decoupling protocols, enhancing qubit coherence. It could also provide fundamental insights into the nature of the underlying dissipative environment.

For classical noise, the frequency-resolved higher noise cumulants (so-called polyspectra [@Rao1984book]) provide a full characterization. Recent work has proposed [@Viola2016; @Ramon2019] and demonstrated [@Oliver2019] measurement protocols to reconstruct polyspectra using a qubit driven by classical non-Gaussian noise. However, the proper generalization of these ideas to truly [*quantum*]{} non-Gaussian noise (as produced by a generic quantum bath) remains an interesting open question; here, the non-commutativity of noise operators at different times poses a challenge as to how best to define polyspectra.

In this paper, we show that the Keldysh approach [@Levitov1996; @Nazarov2003; @BelzigPRL2010; @Clerk2011; @Hofer2017], a method used extensively to characterize low-frequency noise, also provides an unambiguous and practically useful way to describe non-Gaussian quantum bath noise in the frequency domain. It provides a systematic way to construct a quasiprobability distribution to describe the noise, and to assess whether the noise can be faithfully mimicked by completely classical noise processes [@Nazarov2003; @Clerk2011]. It also has a direct operational meaning: the “quantum polyspectra" we introduce are [*exactly*]{} the quantities that contribute to the dephasing of a coupled qubit at each order in the coupling. Moreover, these quantities can be measured using [*the same*]{} non-Gaussian noise spectroscopy techniques designed for classical noise sources [@Viola2016; @Oliver2019; @RBLiu2019]; one does not have to decide in advance whether the noise is classical or quantum to perform the characterization. Note that a recent work presented a method to measure arbitrary quantum bath correlation functions [@RBLiu2019]; in contrast, our work focuses on characterizing the most physically-relevant correlation function at each order and identifying a corresponding quasiprobability.

To highlight the utility of our approach, we apply it to the concrete but non-trivial case of photon shot noise in a driven-damped bosonic mode (a relevant source of dephasing noise in circuit QED systems [@Schoelkopf2006; @Devoret2019] among others). Prior work used the Keldysh approach to study this noise at zero frequency [@Clerk2011; @Clerk2016]; here we instead focus on the behaviour of the frequency-resolved third cumulant, the “quantum bispectrum" (QBS). We show that the QBS reveals important new physics and distinct quantum signatures: at low temperatures, qualitatively new features emerge that would never be present in a classical model with only thermal fluctuations. We also show that the QBS is a generic tool for revealing the breaking of detailed balance and violation of Onsager-like symmetry relations. We find that the photon shot noise QBS violates detailed balance at low temperatures.

[*Keldysh ordering and quantum polyspectra–* ]{} Consider first a classical noise process $\xi(t)$. Its moment generating function (MGF) is defined as $$\begin{aligned}
\label{Eq:MGFcl}
    \Lambda_{\rm class}[F(t); t_f] & = \overline{ \exp \left[ -i \int_0^{t_f} F(t) \xi(t) \right] },\end{aligned}$$ where the bar indicates a stochastic average. Functional derivatives of $\Lambda_{\rm class}$ with respect to $F(t)$ can be used to calculate arbitrary-order correlation functions of $\xi(t)$, while functional derivatives of $\ln \Lambda_{\rm class}$ generate the cumulants of $\xi(t)$ (see e.g. [@JacobsStochasticBook]). Fourier transforming these cumulants yields the polyspectra, which completely characterizes the noise in the frequency domain [@Rao1984book]. In the quantum case, our noise is a Heisenberg picture operator $\hat{\xi}(t)$ whose evolution is generated by the Hamiltonian of some bath; we take $\hat{\xi}(t)$ to be Hermitian for simplicity. Defining correlation functions now has some subtlety, as $\hat{\xi}(t)$ will not in general commute with itself at different times; hence, different time-ordering choices yield different results. Correlation functions at a given order describe both how the bath responds to external perturbations, as well as its intrinsic fluctuations [@Kamenev2011book]. We are interested here in characterizing the latter quantity, and asking whether these fluctuations are equivalent to an effective classical noise process.

The well-developed machinery of Keldysh quantum field theory provides a precise method for accomplishing our task [@Levitov1996; @Nazarov2003; @BelzigPRL2010; @Clerk2011; @Hofer2017]. While this approach is completely general, the simplest derivation is to imagine coupling an ancilla qubit to $\hat{\xi}$, such that the only qubit dynamics is from the interaction picture Hamiltonian $\hat H_{\mathrm{int}} (t) = \frac{1}{2} F(t) \hat{\xi}(t)  \hat \sigma_z$. We then use the dephasing of the qubit to [*define*]{} the MGF of the noise in the quantum case, exactly like we would if the noise were classical: $$\begin{aligned}
\label{Eq:MGFqu}
  &   \Lambda[F( t );t_f] \equiv 
         \langle \hat{\sigma}_{-}(t) \rangle / \langle \hat{\sigma}_{-}(0) \rangle,
         \\
                \label{Eq:MGFK}
         = &
        \textrm{Tr} \left[
            \mathcal{T} e^{- \frac{i}{2} \int_0^{t_f} dt' F(t') \hat{\xi}(t')}
            \hat{\rho}_{\rm B}
            \tilde{\mathcal{T}} e^{-\frac{i}{2} \int_0^{t_f} dt' F(t') \hat{\xi}(t')}
        \right].    
                $$ Here $\hat{\rho}_B$ is the initial bath density matrix, the trace is over bath degrees of freedom, and $\mathcal{T}$ ($\tilde{\mathcal{T}}$) is the time-ordering (anti-time ordering) symbol. Expanding $\Lambda$ in powers of $F(t')$ defines correlation functions at a given order with a particular time-ordering prescription (the so-called Keldysh ordering). We stress that this approach amounts to trying to ascribe the qubit evolution to an effective classical stochastic process; this correspondence then defines cumulants (and implicitly a quasiprobability) for the quantum noise of interest. For truly classical noise, Eq.  reduces to classical MGF in Eq. . The moments and corresponding quasiprobability defined here are intrinsic to the noisy system, and can be used to predict the outcomes of a wide class of schemes designed to measure this noise [@Nazarov2003; @Hofer2017]. They also have a direct role in Keldysh non-equilibrium field theory: they characterize the fluctuations of the “classical" field associated with the operator $\hat{\xi}(t)$. At each order, the Keldysh-ordered correlation function describes the intrinsic fluctuations of the system [@Kamenev2011book]. In contrast, the remaining independent correlation functions at the same order describe how the system responds to external fields which couple to $\hat{\xi}(t)$ (either higher-order response coefficieints, or generalized noise susceptibilities; see Supplemental Material (SM) for details [@SI]). We stress again that at each order, the Keldysh-ordered correlation function is precisely the correlation function “seen" by the qubit.

It follows that the Keldysh-ordered cumulants $C^{(k)}(\vec{t}_{k}) \equiv \langle \langle \hat{\xi}\left(  t_1\right) \cdots 
\hat \xi \left(  t_k\right) \rangle\rangle_\mathcal{K}$ of the noise can be generated from $\chi[F(t);t_f] \equiv \ln \Lambda[ F(t); t_f ] $ via $$\label{eq:ncdef}
    \chi[F(t);t_f]  =  \sum_{\ell=1}^\infty\!\frac{(-i)^\ell }{  \ell !}\prod_{j=1}^\ell
    \left[ \int_0^{t_f}\!\! dt_j F(t_j) \right]\!C^{(\ell)}(\vec{t}_{\ell}),$$ where we define $\vec{t}_n \equiv (t_1,\ldots,t_n)$. Explicit expressions for the first few Keldysh-ordered cumulants are provided in Eqs. (\[SEq:Kcml\_2\]) and (\[SEq:Kcml\_3\]) of the SM [@SI]. The Keldysh-ordered second cumulant is simply a symmetrized correlation function, whereas the third cumulant corresponds to suppressing time-orderings where the earliest operator appears in the middle of an expectation value.

For stationary noise, the $k$-th order cumulant $C^{(k)}(\vec{t}_{k})$ only depends on the $k-1$ time separations $\tau_j \equiv t_{j+1}-t_1$, $j= 1,\ldots, k-1 $. We thus define the quantum polyspectra as Fourier transforms of the Keldysh-ordered cumulants with respect to $\{ \tau_j \}$: $$\begin{aligned}
    &S_{n} [\vec{\omega}_{n}] \equiv \!\int_{\mathbb{R}^n}\! d\vec{\tau}_{n} \,
    e^{-i \vec{\omega}_{n}\cdot\vec{\tau}_{n}} C^{(n+1)}(\vec{\tau}_{n}), \quad n \geq 1. 
\label{eq::polyspectra}\end{aligned}$$ For more discussion on properties of classical polyspectra, see Refs. [@Elgar1994; @Viola2016; @Ramon2019]. The zero frequency limit $\omega_j \to 0 $ of $S_{n} [\vec{\omega}_{n}]$ characterize fluctuations in $\hat m=\int_0^t dt'  \hat \xi ( t' ) $ in the long-time limit (so-called full counting statistics (FCS)). This is the typical setting where the Keldysh approach has found great utility, largely for studying electronic current fluctuations. Here we extend the method to study non-classical, non-Gaussian noise [*at non-zero frequencies*]{} (see also Ref. [@PekolaPRB2006] for an application to frequency-dependent current noise).

![image](S2_10.eps){width="180mm"}

[*Photon shot noise dephasing–* ]{} In what follows, we focus on the physics of the frequency-dependent third cumulant, the so-called quantum bispectrum (QBS) $S_2[\omega_1,\omega_2]$; we drop the subscript $2$ hereafter. The quantum non-Gaussian noise of interest will be the energy fluctuations of a driven-damped bosonic mode. As we will see, the frequency-dependent QBS reveals a host of physics that is not apparent if one only considers the low-frequency behaviour of fluctuations (as has been studied earlier [@Harris2010; @Clerk2011; @Clerk2016]).

Our “bath" here is a driven damped cavity mode $c$ (frequency $\omega_c$, Markovian energy decay rate $\gamma$). Coupling the photon number of this mode to an ancilla qubit, working in a rotating frame at the drive frequency, and letting $\hat{\rho}$ denote the qubit-cavity reduced density matrix, the system dynamics follows the master equation $$\dot {\hat \rho}  =  
    - i [ \hat H_0 + \hat H_{\mathrm{int}} (t) ,\hat \rho ] + \gamma ( {\bar n}_\mathrm{th} + 1 ) \mathcal{D} [ \hat c ]\hat \rho  + \gamma {\bar n}_\mathrm{th} \mathcal{D} [\hat c^\dag ]\hat \rho.
\label{Eq:meq}$$ Here $ \mathcal{D} [ \hat A  ]\hat \rho  = \hat A\rho {\hat A^\dag } -  ( {{\hat A^\dag }\hat A\hat \rho  + \hat \rho {\hat A^\dag }\hat A}  )/2$ is the Lindblad dissipator and ${\bar n}_\mathrm{th} $ the thermal photon number associated with the cavity dissipation. Letting $f$ ($\delta$) denote the drive amplitude (detuning), the cavity Hamiltonian $\hat H_0$ is $$\hat H_0 =- \delta {\hat c^\dag }\hat c - \left( f \hat c + \mathrm{H.c.} \right)  .$$ Finally, to extract the statistics of the photon number $\hat n = \hat c^\dag \hat c$, we couple the driven cavity to an ancilla qubit via $\hat H_{\mathrm{int}} \left( t\right)=    F(t) \hat n \hat \sigma_z /2$. With this setup in place, we solve the master equation in Eq.  to obtain the time-dependent qubit dephasing; from Eq. (\[Eq:MGFK\]), this directly yields the Keldysh-ordered MGF $\Lambda$ for the cavity photon number fluctuations. Even with an arbitrary time-dependent coupling $F(t)$, the qubit dephasing can be solved exactly using an extension of the phase space method in Ref. [@Utami2007] (see also SM [@SI]). An equivalent approach is to calculate correlation functions using standard techniques (e.g. quantum regression theorem, Heisenberg-Langevin equations) [@Gardiner2004], and then apply the Keldysh ordering defined in Eqs.  and . For the remainder of the paper, we will always take the long time limit $t_f \rightarrow \infty$, so that the fluctuations are stationary. One finds the QBS can be written as: $$\begin{aligned}
    S[\omega_1,\omega_2] =
            S_{\rm th}[\omega_1,\omega_2] + 
            S_{\rm dr}[\omega_1,\omega_2],
            \label{eq:QBSDecomposition}\end{aligned}$$ where the first term is completely independent of the drive $f$, and the second term is proportional to $|f|^2$.

[*Quantum bispectrum of drive-independent photon number fluctuations–* ]{} The $f$-independent QBS $S_{\rm th}[\omega_1,\omega_2]$ can be calculated by solving Eq.  with $f=0$. In this case, the cavity relaxes to a thermal steady state with no coherence between different Fock states. Its fluctuations can thus be mapped to a classical Markovian master equation. For such a classical and thermal Markov process, the bispectrum $S_{\mathrm{th} }   [ {\omega}_{1},{\omega}_{2}  ] $ [*must*]{} always be real [@Semerjian2004; @Sinitsyn2016]. In the SM, we also show that in our case, $S_{\mathrm{th} }   [ {\omega}_{1},{\omega}_{2}  ] $ must also be positive semidefinite [@SI]. Letting $\omega_{3} \equiv -\omega_1 -\omega_2$ in all equations that follow, our full calculation for the Keldysh-ordered QBS yields as expected a real, positive function: $$\label{Eq:bispmth}
     S_{\mathrm{th} }   [ {\omega}_{1},{\omega}_{2}  ]  =
        \mathcal{C}_{{\bar n}_\mathrm{th} }  \gamma^2   
        \left( 6 \gamma^2+   \sum\limits_{j=1}^3 { \omega}^2 _j  \right) \Big /  \prod\limits_{j=1}^3 (  \gamma^2  + { \omega}^2 _j   ),$$ with $\mathcal{C}_{{\bar n}_\mathrm{th} }= {\bar n}_\mathrm{th} ( {\bar n}_\mathrm{th} +1  )  ( 2{\bar n}_\mathrm{th} +1) $. We see that the frequency dependence of this contribution to the bispectrum is the same both in the classical high-temperature limit ${\bar n}_\mathrm{th} \to \infty $, and in the extreme quantum limit ${\bar n}_\mathrm{th} \to 0$; the only temperature dependence is in the prefactor. $S_{\mathrm{th} }   [ {\omega}_{1},{\omega}_{2}  ]$ vanishes in the absence of thermal fluctuations (i.e. ${\bar n}_\mathrm{th} \to 0$). In the limit ${\bar n}_\mathrm{th} \to 0$, this expression corresponds (as expected) to the bispectrum of asymmetric telegraph noise (see e.g., [@Sinitsyn2013]), corresponding to fluctuations between the $n=0$ and $n=1$ Fock states. While our general result here suggests that the $\omega$ dependence of the QBS is not sensitive to quantum corrections, we will see that this is not true as soon as a coherent drive is added to the system.

[*Quantum bispectrum of driven photon-number fluctuations–* ]{} We now consider the drive-dependent contribution to the bispectrum, $S_{\mathrm{dr} }   [ {\omega}_{1},{\omega}_{2}  ]$ in Eq. (\[eq:QBSDecomposition\]). This quantity only depends on the drive amplitude $f$ through the overall prefactor ${\bar n}_{ \rm{dr} } =4 |f|^2/(\gamma^2+ 4\delta^2)$ (the intracavity photon number generated by $f$). Note that $S_{\mathrm{dr} }   [ {\omega}_{1},{\omega}_{2}  ]$ remains non-zero at zero temperature, and is the only contribution to the QBS in this limit.

We find that the drive-dependent QBS shows striking quantum signatures. In the classical limit of high temperatures, it is always real and positive (like the purely thermal contribution) [@SI]. However, as temperature is lowered and quantum fluctuations become more important, this quantity can have a negative real part, and even a non-zero imaginary part. These purely quantum features become more pronounced as the magnitude of the drive detuning $\delta$ is increased. The real and imaginary parts of $S_{\mathrm{dr} }   [ {\omega}_{1},{\omega}_{2}  ]$ are plotted for zero temperature in Figs. \[fig:S2\_10\](b) and \[fig:S2\_10\](c) for a large drive detuning ($\delta / \gamma =10$).

We start by discussing the surprising negativity of the real part of the zero-temperature QBS. Negativity in the zero-frequency limit was already discussed in [@Clerk2011; @Clerk2016]. These works showed that this negativity is a purely quantum effect, and that moreover, for large detunings the magnitude of the third cumulant is so great that the fluctuations cannot be described by a positive-definite quasiprobability distribution. Our results show how this striking non-classicality also manifests itself in the [*non-zero*]{} frequency fluctuations. It also demonstrates that this quantum correction has a non-trivial frequency dependence. We find that the quantum bispectrum $ S_{\mathrm{dr}} [ {\omega}_{1},{\omega}_{2}]$ has a different frequency dependence in the quantum limit (${\bar n}_{ \rm{th} }=0$) from its classical limit ${\tilde S}_{ \mathrm{cl}}[ {\omega}_{1},{\omega}_{2}]$. To see this, it is useful to write $$\frac{S_{\mathrm{dr}} [ {\omega}_{1},{\omega}_{2}] }{{\bar n}_{ \rm{dr} }}
     =   ( 2{\bar n}_{ \rm{th} }+1)^2 \!  {\tilde S}_{ \mathrm{cl}}[ {\omega}_{1},{\omega}_{2}] 
     + 
      {\tilde S}_{ \mathrm{q}}[ {\omega}_{1},{\omega}_{2}]  .$$ The first term is the classical contribution which dominates in the high-temperature limit; ${\tilde S}_{ \mathrm{cl}}[ {\omega}_{1},{\omega}_{2}]$ is independent of both $ {\bar n}_{ \rm{dr} } ,{\bar n}_{ \rm{th} }$, and is real and positive for all frequencies. Its form can be found directly from a classical Langevin equation calculation (see Eq.  of the SM [@SI]). In contrast, the second term is the temperature-independent quantum correction. It has a [*completely different*]{} frequency dependence from the classical limit, as described by ${\tilde S}_{ \mathrm{q}}[ {\omega}_{1},{\omega}_{2}]$ $${\tilde S}_{ \mathrm{q}}[ {\omega}_{1},{\omega}_{2}] = -\frac{1}{2 } \! \! \sum\limits_{ \substack{\alpha  \ne \beta \\ \alpha ,\beta  = 1,2,3  }}   \frac{\frac{\gamma }{2} + i{ \omega} _\beta }{ (\gamma  - i { \omega} _\alpha) [( \frac{\gamma }{2} + i { \omega}_\beta  )^2+ \delta ^2]}.
 \label{eq:Sdr2qu}$$ This function can have both a negative real part, and a non-zero imaginary part. In the quantum limit ${\bar n}_{ \rm{th} }=0$, one finds that real part of the QBS only becomes negative above a critical value of the detuning $|\delta|$. Moreover, the initial onset of negativity occurs at $\omega_1 = \omega_2 = 0$. In the large-detuning regime $|\delta| \gg \gamma$, the negative region of the QBS is peaked near a polygon whose shape is defined by the resonance conditions $\omega_j = \pm \delta$ ($j=1,2,3$).

![Frequency dependence of the imaginary parts of the photon-shot noise bispectrum $\mathrm{Im} S [ {\omega}_{1} ,{\omega}_{2} ]$ in the extreme quantum limit ${\bar n}_{ \rm{th} }=0$ at different detunings. Parameters: (a) $\delta / \gamma = 0$, (b) $\delta / \gamma = 1$. Note that for ${\bar n}_{ \rm{th} }=0$ the full quantum bispectrum coincides with the drive-dependent contribution $S_{\rm dr}   [ {\omega}_{1},{\omega}_{2}  ]$.[]{data-label="fig:ImS2"}](ImS2.eps){width="85mm"}

![Time-dependent Keldysh-ordered photon-shot noise third cumulant $C_{\mathrm{dr} }^{(3)}(| t|,|t|)$ in the quantum limit ${\bar n}_{ \rm{th} } =0$ for $t>0$ (red solid lines) and $t<0$ (orange dashed lines). Difference between curves highlights asymmetry under time reversal $t \to -t$. In contrast, the thin blue curves correspond to the same correlator $C_{\mathrm{dr} }^{(3)}(| t|,|t|)$ in the classical limit normalized by thermal photon number, which is symmetric. All correlation functions are normalized by in-cavity drive photon number ${\bar n}_{ \rm{dr} }$. Detunings are: (a) $\delta=0$, and (b) $\delta/\gamma=5$.[]{data-label="fig:timeskn"}](n3K.eps){width="85mm"}

[*Imaginary bispectrum and violations of detailed balance–* ]{} We now turn to another striking feature of the photon shot noise QBS: while in the classical, high-temperature limit it is always real, the quantum correction ${\tilde S}_{ \mathrm{q}}[ {\omega}_{1},{\omega}_{2}]$ has a non-zero imaginary part (see Fig. \[fig:ImS2\]). We stress that this non-trivial imaginary bispectrum can [*only*]{} be probed at finite frequency: by its very definition in Eq. , the imaginary part of the QBS must vanish if either $ \omega_1=0 $ or $\omega_2 =0$. The non-zero imaginary QBS is directly related to the basic symmetries of our quantum noise process, in particular the violation of higher-order Onsager reciprocity relations [@Onsager1931; @Semerjian2004; @Sinitsyn2016]. If a temporal cumulant $C^{(n+1)}(\vec{\tau}_{n})$ is invariant under $\vec{\tau}_{n} \to -\vec{\tau}_{n}$ (i.e. exhibits Onsager reciprocity), then the corresponding polyspectrum must be real [@Brillinger1967]. Further, a classical Markov process obeying detailed balance always respects this symmetry. In our system quantum corrections (as described by ${\tilde S}_{ \mathrm{q}}[ {\omega}_{1},{\omega}_{2}]$) cause a breaking of this symmetry and hence of detailed balance. There is a long history of studying detailed balance in non-equilibrium driven-dissipative quantum systems (see e.g. [@Kubo1966; @Agarwal1973; @Carmichael1976; @Tomita1973; @Tomita1974; @Carmichael2002]); the QBS provides yet another tool for exploring this physics. In the SM, we discuss another related quantum system which exhibits an apparent breaking of detailed balance, namely a cavity driven by squeezed noise [@SI].

For a heuristic understanding of this symmetry breaking, we consider a simpler object, the temporal (Keldysh-ordered) third cumulant $C^{(3)}(\tau_1,\tau_2)$ at $\tau_1 = \tau_2 = t$. The non-zero imaginary QBS implies that this correlator differs for $t$ and $-t$ (see Fig. \[fig:timeskn\]). Using the definition of Keldysh ordering in Eqs. (\[Eq:MGFK\]) and (\[eq:ncdef\]) we find: $$\begin{aligned}
     &  \langle \delta \hat n(0) \delta \hat n(t ) \delta \hat n(t )\rangle_\mathcal{K} =  
    \label{eq:KeldyshOrderedSkewness}  \\
    & \frac{1 }{2} \langle \{\delta \hat n(0) ,[\delta \hat n(t )]^2 \} \rangle
           -\frac{ \Theta(-t)}{4}   \langle [ \delta \hat n(t), [ \delta \hat n(t ),\delta \hat n(0)] ] \rangle
  \nonumber ,\end{aligned}$$ where $\delta \hat{n}(t) = \hat{n}(t) - \langle \hat{n}(t) \rangle$, and $\Theta(t)$ is the Heavisde step-function. One finds that any imaginary quantum correction $ \mathrm{Im}  \tilde{S}_q[\omega_1,\omega_2]$ is entirely due to the second term on the RHS; it is thus completely responsible for the the lack of time-symmetry. What does this mean physically? As we have emphasized, the Keldysh ordering is relevant for any measurement protocol that directly probes $\hat{n}(t)$ [@Nazarov2003; @Hofer2017]. In contrast, the first term on the RHS of Eq. (\[eq:KeldyshOrderedSkewness\]) would be relevant if we correlated a measurement of $\delta \hat{n}$ with a separate, direct measurement of $\delta \hat{n}^2$ (i.e. the Keldysh approach would give this answer for this sort of setup [@Hofer2017]). These protocols are not equivalent: measuring $\delta \hat{n}$ and then squaring the result has a different back-action than if one directly measured $\delta \hat{n}^2$. The latter measurement provides less information (and hence has less backaction), as it provides no information on the sign of $\delta \hat{n}$. This now provides a heuristic way of understanding the second term on the RHS of Eq. (\[eq:KeldyshOrderedSkewness\]) (and the consequent lack of time symmetry). For $t<0$, one is first measuring $\delta \hat{n}^2$. As a result, the two measurement protocols have different backaction effects, and the two correlation functions are distinct. In contrast, for $t>0$, the earlier measurement is the same in both protocols, hence the backaction effect is identical, and the two protocols agree.

While our heuristic explanation here invokes measurement backaction, we stress yet again that the Keldysh-ordered correlation function is an [*intrinsic*]{} property of the driven cavity system [@Kamenev2011book; @Nazarov2003], with a relevance that goes beyond the analysis of just a single measurement setup. Further, this is the ordering that is “chosen" by our qubit: if one simply interprets the qubit dephasing as arising from classical noise, then the Keldysh ordered bispectrum (with its imaginary part) plays the role of the bispectrum of this effective classical noise.

[*Conclusions–* ]{} We have shown how the Keldysh approach to quantum noise provides a meaningful way to define the polyspectra of non-classical, non-Gaussian noise. In the experimentally-relevant case of photon shot noise fluctuations in a driven-damped resonator, the quantum bispectrum reveals distinct quantum features and a surprising quantum-induced breaking of detailed balance. We stress that our approach amounts to interpreting the dephasing of a qubit by quantum noise as arising from an effective classical noise process. As such, the same noise spectroscopy techniques that have been used successfully to measure classical bispectra with qubits [@Viola2016; @Oliver2019] can be directly used (without modification) to measure our quantum bispectra.

This work was supported as part of the Center for Novel Pathways to Quantum Coherence in Materials, an Energy Frontier Research Center funded by the U.S. Department of Energy, Office of Science, Basic Energy Sciences.

**Supplementary Material: Spectral characterization of non-Gaussian quantum noise**

Distinguishing fluctuations from response properties
====================================================

As discussed in the main text, a general $n$-point quantum correlation function describes both the intrinsic fluctuation properties of the system of interest (i.e. quantities that play the role of classical noise), as well as the response properties of the system to external applied fields. The situation is very clear at second order, where the product ${\hat{\xi}}(t) {\hat{\xi}}(t')$ can be decomposed as the sum of a commutator and an anti-commutator. The commutator determines the retarded Green function $$G^R(t) \equiv -i \Theta(t) \left \langle [ {\hat{\xi}}(t), {\hat{\xi}}(0) ] \right \rangle.$$ This describes how the average value $\langle {\hat{\xi}}(t) \rangle$ changes to first order in response to an external perturbing field $V(t)$ entering the Hamiltonian as $$\hat{H}_{\rm ext}(t) = V(t) {\hat{\xi}}.$$ The relevant Kubo formula is: $$\delta \langle {\hat{\xi}}(t) \rangle = 
        \int_{-\infty}^{\infty} dt' G^R(t-t') V(t').$$ In contrast, the anti-commutator describes the symmetrized noise spectral density: $$S[\omega] \equiv \frac{1}{2} \int dt e^{i \omega t} \langle \{ {\hat{\xi}}(t), {\hat{\xi}}(0) \} \rangle .$$ As has been discussed in many places (see e.g. Ref. [@ClerkRMP]), this spectral density plays the role of a classical noise spectral density.

The Keldysh technique provides an unambiguous way of extending this separation between noise and response to higher orders. A full exposition of this method is beyond the scope of this paper; we refer the reader to Ref. [@Kamenev2011book]. We sketch the main ideas needed here. In the path-integral formulation of the Keldysh technique, each operator corresponds to two different fields, the classical field $\xi_{\rm cl}(t)$ and the quantum field $\xi_{\rm q}(t)$. Averages of these fields (weighted by the appropriate Keldysh action describing the system) then correspond to operator averages with a particular time ordering. One finds that:

-   Averages only involving quantum fields are necessarily zero.

-   Averages involving at least one classical field $\xi_{\rm cl}(t)$ and one or more quantum fields $\xi_{\rm q}(t)$ can always be interpreted as response coefficients to an external perturbation of the form $\hat{H}_{\rm ext}(t)$.

-   Averages [*only*]{} involving classical fields $\xi_{\rm cl}(t)$ do not correspond to any kind of response function. Instead, they describe the intrinsic fluctuation properties of the system

Formally, this dichotomy arises because the perturbation $\hat{H}_{\rm ext}(t)$ enters the action of the system as a term that [*only*]{} involves the quantum field, i.e. $S_{\rm ext} = \int dt V(t) \xi_{\rm q}(t)$. Perturbation theory in $V(t)$ thus necessarily introduces powers of the quantum field. For example, at second order we have:

-   The average $\overline{\xi_{\rm cl}(t) \xi_{\rm q}(0)}$ is directly proportional to the retarded Green function $G^R(t)$, and thus describes linear response to the external field.

-   The average $\overline{\xi_{\rm cl}(t) \xi_{\rm cl}(0)}$ is proportional to $\langle \{ {\hat{\xi}}(t), {\hat{\xi}}(0) \} \rangle$ and thus determines the usual symmetrized noise spectral density.

The same decomposition applies at higher orders. Consider third order correlators. The average of three classical fields $\overline{\xi_{\rm cl}(t_1) \xi_{\rm cl}(t_2) \xi_{\rm cl}(0)}$ is precisely the Keldysh ordered correlator discussed in the main text; it cannot be associated with a response coefficient. The remaining non-zero correlators describe different kinds of response:

-   The average $\overline{\xi_{\rm cl}(t) \xi_{\rm q}(t') \xi_{\rm q}(t'')}$ represents a second-order Kubo response coefficient. It determines to second order how $\langle {\hat{\xi}}(t) \rangle$ is modified by $\hat{H}_{\rm ext}(t')$ at earlier times (i.e. how it depends on $V(t')$ and $V(t'')$).

-   The average $\overline{\xi_{\rm cl}(t) \xi_{\rm cl}(t') \xi_{\rm q}(t'')}$ describes a first order noise-susceptibility [@Reulet2007]. It determines how the symmetrized correlator $\langle \{ {\hat{\xi}}(t), {\hat{\xi}}(t') \} \rangle$ is modified to first order by $\hat{H}_{\rm ext}(t'')$

The arguments sketched here provided perhaps the deepest justification for considering Keldysh ordered correlation functions: they provide a clear and unambiguous way to distinguish fluctuation properties from response properties. We stress that an arbitrary correlation function could always be written as a linear combination of the Keldysh-ordered correlator (which describes pure noise) and additional terms describing response properties.

Explicit expressions for the second and third Keldysh-ordered cumulants
=======================================================================

For concreteness, here we provide explicit expressions for the first few Keldysh-ordered cumulants $C^{(k)}(\vec{t}_{k}) \equiv \langle \langle \hat{\xi}\left(  t_1\right) \cdots \hat \xi \left(  t_k\right) \rangle\rangle_\mathcal{K}$ defined by Eqs.  and in the main text. The second order cumulant function $C^{(2)}(\vec{t}_{2}) $ is just the auto-correlation function of $ \hat \xi (  t )$ $$\label{SEq:Kcml_2}
C^{(2)}(\vec{t}_{2})  =    \langle \langle \hat \xi\left( {t_1} \right)  \hat \xi\left( {t_1} \right) \rangle \rangle_\mathcal{K} =  \frac{1}{2}\langle\{ \delta \hat \xi\left( {t_1} \right) , \delta \hat \xi\left( {t_1} \right) \}\rangle,$$ where $\delta \hat \xi=\hat \xi -\langle \hat \xi \rangle  $. However, the third cumulant corresponds to a more complex ordering $$\label{SEq:Kcml_3}
C^{(3)}(\vec{t}_{3})  =    \langle \langle \hat \xi\left(  {t_1}  \right)  \hat \xi\left( {t_2}  \right) \hat \xi\left(  {t_3} \right)\rangle \rangle_\mathcal{K} =  \frac{1}{4}\sum_{\vec{\pi}_3 \in \mathcal{P}_3}  [1-\Theta( t_{\pi _1} - t_{\pi _2}) \Theta( t_{\pi _3} - t_{\pi _2})]  \langle \delta \hat \xi (t_{\pi _1}) \delta \hat \xi (t_{\pi _2}) \delta \hat \xi (t_{\pi _3})  \rangle,$$ where $\mathcal{P}_3$ denotes the set of all possible permutations of $(123)$ indices. Such ordering is given by an average over all permutations of the three displaced operators $\delta \hat \xi (t_j)$, except for the terms where the earliest time appears in the middle position (as implied by the step functions), in agreement with expansion of the operator in Eq.  in powers of coupling $F(t')$. A similar expression of Keldysh-ordered third cumulant has also been derived for current operators in Ref. [@PekolaPRB2006].

Phase space method for computing Keldysh-ordered cumulants of a driven damped cavity
====================================================================================

In this section, we outline the phase space method to calculate Keldysh-ordered cumulants. However, we remark that once we have defined the unique Keldysh ordering for each higher cumulant using Eqs.  and , standard techniques for computing multi-point correlation functions (e.g. Langevin equations of motion, and quantum regression theorem) work equally well for the Keldysh-ordered cumulants.

In the phase space method, we need to solve the time evolution of qubit coherence $\hat \rho _{ \uparrow  \downarrow } (t) \equiv \langle\uparrow \! |\hat \rho (t)|\! \downarrow \rangle$ $$\dot {\hat \rho}_{ \uparrow  \downarrow }  =  
    - i [ \hat H_0  ,\hat \rho_{ \uparrow  \downarrow } ]  - i \{ \hat H_{\mathrm{int}}(t) ,\hat \rho_{ \uparrow  \downarrow } \}+ \gamma ( {\bar n}_\mathrm{th} + 1 ) \mathcal{D} [ \hat c ]\hat \rho_{ \uparrow  \downarrow }  + \gamma {\bar n}_\mathrm{th} \mathcal{D} [\hat c^\dag ]\hat \rho_{ \uparrow  \downarrow },
\label{SEq:meqcoh}$$ which is a direct extension, with a time-modulation $F(t)$ in interaction $\hat H_{\mathrm{int}} \left( t\right)=  {\lambda} F(t) \hat n \hat \sigma_z /2$, of the technique used in Ref. [@Utami2007]. Here we use a constant coefficient $\lambda$ to keep track of orders in expansion on the coupling; by the end of the calculation, one can always set $\lambda=1$. We stress that if we replace the time-independent coupling $\lambda$ with a time-dependent one, the relevant derivations in Ref. [@Utami2007] still hold rigorously, and we refer interested readers to this paper for more detail.

Without loss of generality, the system initial state can be chosen as a product state between the qubit and the cavity, with the cavity in thermal equilibrium. Thus, Wigner function $W(x,p;t)$ of the coherence operator $\hat \rho _{ \uparrow  \downarrow } (t) $ is Gaussian throughout time evolution. Moreover, for the Fourier transform of $W(x,p;t)$, we can assume the following ansatz [@Utami2007] $$W\left[ {k,q;t} \right] =  e^{ - \nu ( t  )} \exp \left(   - i [ {k\bar x (t) + q\bar p (t)}  ] - \frac{1}{2} (k^2+q^2) \sigma _s (t) \right),$$ from which the moment generating function can be computed as $\Lambda[F( t );t_f] = e^{-\nu ( t_f )} $. Substituting this ansatz into the master equation in Eq. , we then need to solve a set of ordinary differential equations for the coefficient functions

\[SEq:cohODE\] $$\begin{aligned}
&\dot \nu_\mathrm{th}  =    i\lambda F ( t ) \left( {\sigma _s- \frac{1}{2}} \right) ,\\
&  {\dot \sigma }_s  = \gamma \left( {\bar n}_\mathrm{th}+ \frac{1}{2}  \right) - \gamma {\sigma _s}   - iF ( t ) \lambda \sigma _s^2 + \frac{ i\lambda F \left( t\right)  }{4} , \\
&\dot \nu_\mathrm{dr}  =     \frac{{i\lambda }}{2} F(t)  ({\bar x}^2 +{\bar p}^2)  ,\\
& \dot {\bar x} =  - \delta \bar p + \sqrt 2 \, \mathrm{Im}f  - iF ( t ) \lambda {\sigma _s}  \bar x- \frac{\gamma }{2}\bar x , \label{SEq:mfGSantcoef-xb}  \\
&  \dot {\bar p}  = \delta \bar x + \sqrt 2 \, \mathrm{Re}f  - i F (t) \lambda {\sigma _s} {\bar p} - \frac{\gamma }{2}\bar p , \label{SEq:mfGSantcoef-pb} 
  \end{aligned}$$

where the exponent $\nu (t)=\nu_\mathrm{th}(t) +\nu_\mathrm{dr}(t)$ can be automatically written as a sum of drive-independent and drive-dependent parts.

The Keldysh-ordered cumulants $C^{(\ell)}(\vec{t}_{\ell})$ can now be extracted using the equation (see Eq.  in the main text) $$\label{SEq:ncdef}
    \nu (t_f)  = - \sum_{\ell=1}^\infty\!\lambda^\ell\frac{ (-i)^\ell }{  \ell !}\prod_{j=1}^\ell
    \left[ \int_0^{t_f}\!\! dt_j F(t_j) \right]\!C^{(\ell)}(\vec{t}_{\ell}),$$ i.e. the cumulants can be obtained by solving Eqs.  perturbatively in orders of $\lambda$, and comparing the results to the integrals above. Since the cumulant functions $C^{(\ell)}(\vec{t}_{\ell})$ must be symmetric over permutations of its variables $\{\vec{t}_{\ell} \}$, such procedure will lead to a unique result. For example, for the photon shot noise in a driven damped cavity discussed in the main text, first few drive-independent contributions to cumulants are given by

$$\begin{aligned}
& C_{\mathrm{th} }^{(1)}( t_1) = {\bar n}_\mathrm{th}  ,\\
& C_{\mathrm{th} }^{(2)}(\vec{t}_2) ={\bar n}_\mathrm{th} ( {\bar n}_\mathrm{th} +1 ) e^{-\gamma | t_1 - t_2 |}, \\
&C_{\mathrm{th} }^{(3)}( \vec{t}_{3}) =  {\bar n}_\mathrm{th} ( {\bar n}_\mathrm{th} +1 )  ( 2{\bar n}_\mathrm{th} +1 ) \exp \left(  - \frac{\gamma }{2} | t_1 - t_2 | - \frac{\gamma }{2}| t_2 - t_3 | - \frac{\gamma }{2}| t_1 - t_3 |  \right) . \label{SEq:cmlth_3} 
  \end{aligned}$$

Taking Fourier transform of Eq.  for the third cumulant, we obtain the drive-independent QBS, as given by Eq.  in the main text.

Proof of non-negative energy shot noise bispectrum in a classical driven damped oscillator
==========================================================================================

In the classical limit ${\bar n}_\mathrm{th} \to \infty $, the cavity mode annihilation operator $\hat c$ in the main text can be described by a classical stochastic variable $c(t)$, describing the amplitude of a driven damped classical harmonic oscillator. The equation of motion is now given by $$d c  =   - ( {\gamma }/{2}-i\delta )c dt + if dt+ \sqrt{ \gamma{\bar n}_\mathrm{eff} }  dW ,$$ where ${\bar n}_\mathrm{eff} ={\bar n}_\mathrm{th} +1/2 \simeq {\bar n}_\mathrm{th}$ (the $1/2$ correction is added so that second-order correlators between $c(t)$, $c^*(t')$ match their symmetrized quantum counterparts), and $dW$ is a complex-valued Wiener increment. The solution to this stochastic differential equation can be written as $c(t) = c_0 + \zeta (t)$, where $c_0 = \langle c(t) \rangle$ is a complex constant number, and $\zeta (t)$ is a complex zero-mean stochastic variable. In the long time limit, $\zeta (t)$ is Gaussian and stationary, satisfying the equation $$\langle\zeta^*(t)\zeta(t') \rangle  ={\bar n}_\mathrm{eff}  \exp \!\left[-i \delta (t-t') -\frac{\gamma}{2} |t-t'|\right],$$ whereas all other second correlators vanish $\langle\zeta(t)\zeta(t') \rangle =\langle\zeta^*(t)\zeta^*(t') \rangle^* \equiv 0 $. The photon number operator $\hat n$ then corresponds to the energy of the classical oscillator $n (t)=|c(t)|^2$, so that its Fourier transform can be expressed using Fourier components of $\zeta(t)$ as $$n [\omega] = \int dt e^{i \omega t} n(t)= |c_0|^2 +\int d\omega' \zeta^* [\omega - \omega']\zeta[ \omega']+c^*_0 \zeta[ \omega] +c_0 \zeta^*[ \omega] .$$

Since the Fourier transform $\zeta [\omega]$ of a Gaussian variable must also be Gaussian, polyspectra of $n(t)$ can be calculated using the expression above by applying Wick’s theorem. Noting that all the anomalous correlators vanish, the only contractions that contribute would be given by terms of the following form $$\langle\zeta^*[\omega]\zeta[\omega'] \rangle  =   \frac{\gamma {\bar n}_\mathrm{eff}}{\left(\omega - \delta\right)^2+\left(\frac{\gamma}{2}\right)^2}   \delta (\omega +\omega')  ,$$ which is always non-negative. It is then straightforward to show that both drive-independent and drive-dependent contributions to polyspectra must also be non-negative for all frequencies. In particular, the frequency dependence ${\tilde S}_{ \mathrm{cl}}[ {\omega}_{1},{\omega}_{2}]$ of the drive-dependent bispectrum in the classical limit (see main text for definition) is real and positive semidefinite, which can be explicitly written as $$\begin{aligned}
{\tilde S}_{ \mathrm{cl}}& [ {\omega}_{1},{\omega}_{2}] = \frac{1}{4}  \sum\limits_{ \alpha  \ne \beta ; \,\alpha ,\beta  = 1,2,3 } \frac{\gamma ^2}{ \left[ {\left( \frac{\gamma }{2} \right)}^2  + {\left( \omega _\alpha  +\delta \right)}^2 \right]\left[ {\left( \frac{\gamma }{2} \right)}^2  + {\left( \omega _\beta - \delta \right)}^2  \right] } .
\label{Seq:QBSdrcl}\end{aligned}$$

Temporal skewness for squeezed bath photon fluctuations
=======================================================

In the main text, we show a violation of higher-order Onsager reciprocity relations solely due to quantum corrections in the temporal third cumulant (skewness), which can be probed by an imaginary part in the QBS. Here we provide an example where the temporal skewness exhibits time asymmetry in both the classical and the quantum limits, and the skewness function also reveals insights into non-equilibrium dynamics in well-defined classical systems. We again consider photon shot noise in a dissipative bosonic mode, but now driven by squeezed noise. The master equation is $$\dot {\hat \rho}  =  - i [ \hat H_0 + \hat H_{\mathrm{int} },\hat \rho ] + \gamma ( {\bar n}_\mathrm{cl} + 1 ) \mathcal{D} [ \hat s_r ]\hat \rho  + \gamma {\bar n}_\mathrm{cl} \mathcal{D} [\hat s_r^\dag ]\hat \rho ,
\label{Seq:meqSQ}$$ where $\hat s_r = \hat c \cosh{r} + \hat c^\dag \sinh r$ denotes the squeezed bath operator. In the rotating frame, the oscillator Hamiltonian is $\hat H_0 =- \delta {\hat c^\dag }\hat c$, and its interaction with the qubit is $\hat H_{\mathrm{int}} (t) = \frac{1}{2} F(t) \hat{n}(t)  \hat \sigma_z$. Such noise model has a well-defined classical limit if we let ${\bar n}_\mathrm{cl} \to \infty$, where the bosonic mode can be equivalently described by a classical stochastic variable $c(t)$. We note that the steady state of the corresponding classical model is not thermal equilibrium, enabling a violation of Onsager-like relations even in the classical limit.

For concreteness, we again consider the temporal third cumulant $C ^{(3)}( t,t)$, which can be written as a sum of classical and quantum contributions as $$C ^{(3)}( t,t)= (2{\bar n}_\mathrm{cl} + 1)^3 f( t )  \left[ {\tilde C}_{\mathrm{cl}}^{(3)}(t) - \frac{1}{(2{\bar n}_\mathrm{cl} + 1)^2}\right],$$ where $f( t )=e^{ - \gamma | t |}\cosh (2r )/4$ is an even function of time $t$ and independent of ${\bar n}_\mathrm{cl}$. The coefficient function ${\tilde C}_{\mathrm{cl}}^{(3)}(t)$ for the classical contribution is given by $${\tilde C}_{\mathrm{cl}}^{(3)}(t) =\cosh^2 (2r)+  \frac{  \gamma ^2 \sinh^2 (2r) }{ \gamma ^2+ 4\delta ^2} [ 1+  2 \cos  ( \delta t +\delta  | t  | )   ].$$ The situation is now reversed: the quantum correction is symmetric under time reversal $t \to -t$, whereas the classical contribution is asymmetric for a generic nonzero detuning $\delta \ne 0$.

The time asymmetry in $C ^{(3)}( t,t)$ has its roots in classical non-equilibrium dynamics: in the classical limit ${\bar n}_\mathrm{cl}  \gg 1$, we can introduce two real quadratures $x$ and $p$ defined by $c=( x + i p)/ \sqrt 2$ to describe the corresponding classical oscillator. Their dynamics satisfies the stochastic differential equations

$$\begin{aligned}
& dx = (-\delta p - \frac{\gamma }{2} x ) dt + e^{r} \sqrt {\gamma {\bar n}_\mathrm{eff}} dW_1  , \\
& dp = ( \delta x - \frac{\gamma }{2} p ) dt + e^{-r} \sqrt {\gamma {\bar n}_\mathrm{eff}} dW_2 ,\end{aligned}$$

where ${\bar n}_\mathrm{eff} ={\bar n}_\mathrm{cl} +1/2 \simeq {\bar n}_\mathrm{cl}$, and $dW_1 $ and $dW_2$ are independent Wiener increments. These equations formally also describe time evolution of a resonantly coupled pair of real harmonic modes, where the interaction strength is given by $|\delta |$, and each oscillator is also coupled to a thermal reservoir with thermal excitations $e^{\pm 2 r} {\bar n}_\mathrm{eff}$. This coupled two-mode system for $r \ne 0$ is a typical example of non-equilibrium system that violates detailed balance, manifested as time asymmetry in cross correlation functions $\langle A(0)B(t)\rangle$ [@Tomita1973; @Tomita1974; @Carmichael2002]. Noting that $n(t)$ corresponds to the total energy in the classical limit, the skewness $C ^{(3)}( t,t)$ can then be viewed as a correlation function between energy fluctuations $\delta n(0)$ and its higher order fluctuations $[\delta  n(t )]^2$ at a different time. Thus, the time asymmetry in $C ^{(3)}( t,t)$ is again a signature of detailed balance violation, which in turn is due to the imbalanced thermal baths set by the nonzero $r$.
