Correlators in simultaneous measurement of non-commuting qubit observables
Juan Atalaya Shay Hacohen-Gourgy Leigh S. Martin Irfan Siddiqi Alexander N. Korotkov













    Department of Electrical and Computer Engineering, University of California, Riverside, CA 92521, USAQuantum Nanoelectronics Laboratory, Department of Physics, University of California, Berkeley CA 94720, USACenter for Quantum Coherent Science, University of California, Berkeley CA 94720, USA.Quantum Nanoelectronics Laboratory, Department of Physics, University of California, Berkeley CA 94720, USACenter for Quantum Coherent Science, University of California, Berkeley CA 94720, USA.Quantum Nanoelectronics Laboratory, Department of Physics, University of California, Berkeley CA 94720, USACenter for Quantum Coherent Science, University of California, Berkeley CA 94720, USA.Department of Electrical and Computer Engineering, University of California, Riverside, CA 92521, USA
    
    We consider the simultaneous and continuous measurement of qubit observables σ_z and
    σ_zcosφ + σ_xsinφ, focusing on the temporal correlations of the two output signals. Using quantum Bayesian theory, we derive analytical expressions for the correlators, which we find to be in very good agreement with experimentally measured output signals.
    We further discuss how the correlators can be applied to parameter estimation, and use them to infer a small residual qubit Hamiltonian arising from calibration inaccuracy in the experimental data.
    Correlators in simultaneous measurement of non-commuting qubit observables
        Alexander N. Korotkov
        December 30, 2023
    ==========================================================================
    
    
    
    Continuous quantum measurements (CQMs) have become a unique platform to explore fundamental aspects of quantum phenomena and have potential applications to quantum information science.
    They have been discussed theoretically
    in various contexts
    (e.g., <cit.>), and in the past decade superconducting qubits have become the main experimental system for the realization of CQMs <cit.>. CQMs are shedding new light on our understanding of the quantum measurement process, and there is also a growing interest in CQM applications, including quantum feedback <cit.>, rapid state purification <cit.>, preparation of entangled states <cit.>, and continuous quantum error correction <cit.>.
    
    While a simultaneous measurement of non-commuting observables is impossible with projective measurements, nothing theoretically forbids such a measurement using CQMs. Aside from new physics, such a protocol may open up new areas of applications, inaccessible to projective measurements.
    The theoretical discussion of a simultaneous measurement of incompatible observables has a long history <cit.>. For the measurement of non-commuting observables of a qubit, statistics of time-integrated detector outputs and fidelity of state monitoring directly via time-integrated outputs has been analyzed in Ref. <cit.>. The evolution of the qubit state due to simultaneous measurement of incompatible variables has been described theoretically in Ref. <cit.>, and has been recently demonstrated experimentally in Ref. <cit.>.
    
    
    In this letter, we focus on the temporal correlations of the output signals from two linear detectors measuring continuously and simultaneously the qubit observables σ_z and σ_φ≡σ_zcosφ + σ_xsinφ, where σ_x and σ_z are the Pauli matrices and φ is an angle between the two measurement directions on the Bloch sphere (Fig. <ref>). The experimental setup is described in detail in Ref. <cit.>; it is based on a Rabi-rotated physical qubit, which is measured stroboscopically <cit.> using symmetric sideband pumping of a coupled resonator, so that σ_z and σ_φ for an effective rotating-frame qubit are being measured. Description of such a measurement based on the theory of quantum trajectories <cit.> has been developed in Ref. <cit.>. In this letter we will use a simpler approach based on quantum Bayesian theory <cit.>. The quantum Bayesian description of the rotating-frame experiment <cit.> is developed in the Supplemental Material <cit.>. The goal of this letter is calculation of the time-correlators for the output signals measuring σ_z and σ_φ, and their comparison with experimental data.
    As we will see, these correlators may be a useful tool for sensitive parameter estimation in an experiment. These correlators are also important in the analysis of quantum error detection and correction based on simultaneous measurement of non-commuting operators <cit.>. We note that the analyzed output signal correlators are different from qubit-state correlators <cit.>.
    
     Quantum Bayesian theory.—A simultaneous continuous measurement of the qubit observables σ_z and σ_φ by two linear detectors (Fig. <ref>) produces noisy output signals I_z(t) and I_φ(t), respectively <cit.>,
    
        I_z(t) =     Tr[σ_zρ (t)] + √(τ_z) ξ_z(t), 
        
        I_φ(t) =     Tr[σ_φρ (t)] + √(τ_φ) ξ_φ(t),
    
    where ρ (t) is the qubit density matrix and τ_z, τ_φ are the “measurement” (collapse) times needed for an informational signal-to-noise ratio of 1 for each channel. Note the chosen normalization for I_z and I_φ.
    In the Markovian approximation, the noises ξ_z and ξ_φ are uncorrelated, white, and Gaussian with two-time correlators
    
        
        ⟨ξ_z(t) ξ_z(t')⟩ = ⟨ξ_φ(t) ξ_φ(t')⟩ = δ(t-t')
    
    and ⟨ξ_z(t) ξ_φ(t')⟩ = 0. The qubit state is characterized in the Bloch-sphere representation as ρ (t)≡[ + x(t) σ_x + y(t) σ_y + z(t) σ_z]/2. We assume phase-sensitive amplifiers in the experimental setup, amplifying the optimal (informational) quadratures, so that the qubit evolution due to measurement is not affected by the phase backaction related to fluctuations in the orthogonal (non-informational) quadrature <cit.>. Then there is only the quantum informational backaction, which for measurement of σ_z and σ_φ is described <cit.> by the evolution equations (in the Itô interpretation)
    
        
        ẋ=     - Γ_z x - Γ_φcosφ(xcosφ - zsinφ)  - τ_z^-1/2xz  ξ_z  
              - τ_φ^-1/2[xzcosφ - (1-x^2)sinφ] ξ_φ , 
        ẏ=     - (Γ_z +Γ_φ )  y - τ_z^-1/2yz  ξ_z
                   
              - τ_φ^-1/2 y [zcosφ + xsinφ] ξ_φ , 
        ż=    Γ_φsinφ(xcosφ - zsinφ) +
                      τ_z^-1/2(1-z^2) ξ_z    
             + τ^-1/2_φ[(1 - z^2)cosφ -xzsinφ]ξ_φ .
    
    Here Γ_z and Γ_φ are the ensemble dephasing rates due to measurement, so that the quantum efficiencies <cit.> for the two channels are
    η_z = 1/(2τ_zΓ_z) and η_φ = 1/(2τ_φΓ_φ). In the experiment η_z =0.49 and η_φ=0.41.
    
    
    Equations (<ref>)–(<ref>) describe qubit evolution only due to measurement. We also need to add terms due to unitary evolution and due to decoherence not related to measurement. We assume the Hamiltonian H=ħΩ̃_ Rσ_y/2, describing Rabi oscillations about y-axis with frequency Ω̃_ R. In the experiment, Ω̃_ R=Ω_ R-Ω_ rf is a small (kHz-range) undesired mismatch between the physical qubit Rabi frequency Ω_ R and rotating frame frequency Ω_ rf defined by detuning of sideband pumps <cit.>. Decoherence of the effective qubit arises from the decoherence of the physical qubit with energy relaxation time T_1 and dephasing time T_2 [the pure dephasing rate is then T_ pd^-1=T_2^-1-(2T_1)^-1]. Averaging the decoherence over fast rotation Ω_ R≫ T_2^-1 and adding unitary evolution, we obtain <cit.>
        ẋ = Ω̃_ R z - γ x,   ẏ = -T_2^-1 y,    ż = -Ω̃_ R x -γ z,
         
           γ = (T_1^-1+T_2^-1)/2,
    
    Evolution of the effective qubit is described by adding terms from Eqs. (<ref>)–(<ref>) and (<ref>).
    
    
    
    
     Correlators.—Our next goal is to calculate the two-time correlators, K_ij(τ), for the output signals,
    
        
        K_ij(τ) ≡⟨ I_j(t_1+τ)  I_i(t_1)⟩,     τ>0,       i,j ∈{z,φ}.
    
    Self- and cross-correlators correspond to i=j and i≠ j, respectively. The averaging in Eq. (<ref>) is over an ensemble of measurements with the initial qubit state ρ _ in prepared at time t_ in≤ t_1. We will see, however, that  the result does not depend on ρ_ in, t_ in, and t_1, so Eq. (<ref>) can also be understood as averaging over time t_1. We assume that the parameters in Eqs. (<ref>)–(<ref>) do not change with time. By assuming τ>0, we avoid considering the trivial zero-time contribution to the self-correlators, Δ K_ii(τ )=τ_i δ (τ).
    
    As shown in the Supplemental Material <cit.>, calculation of the correlators from Eqs. (<ref>)–(<ref>) is equivalent to the following  recipe<cit.>: we replace an actual continuous measurement at the (earlier) time moment t_1 with a projective measurement of σ_i, so that the measurement result I_i(t_1) is ± 1 with probability {1± Tr[σ_i ρ(t_1)]}/2, and the qubit state collapses correspondingly to the eigenstate |1_i⟩ or |0_i⟩ of σ_i (σ_i|1_i⟩ =|1_i⟩, σ_i|0_i⟩ =-|0_i⟩). This gives the correlator
        
        K_ij(τ) =  Tr[σ_j  ρ_ av(t_1+τ|1_i)] 1+  Tr[σ_i ρ(t_1)]/2
           -   Tr[σ_j  ρ_ av(t_1+τ|0_i)]  1-  Tr[σ_i  ρ(t_1)]/2 ,
    
    where ρ_ av(t_1+τ|1_i) is the ensemble-averaged density matrix at time t_1+τ with the initial condition ρ_ av(t_1|1_i)=|1_i⟩⟨ 1_i|; similarly,  ρ_ av(t_1+τ|0_i) starts with ρ_ av(t_1|0_i)=|0_i⟩⟨ 0_i|. The evolution of ρ_ av is given by Eqs. (<ref>)–(<ref>) without noise, ξ_z=ξ_φ=0 (because of the Itô form), so that
     
        ẋ_ av=     - Γ_z x_ av - Γ_φcosφ( x_ avcosφ - z_ avsinφ)
            
            +Ω̃_ R z_ av -γ x_ av,
            
        ẏ_ av=     - (Γ_z +Γ_φ )  y_ av -T_2^-1 y_ av  ,
            
        ż_ av=    Γ_φsinφ( x_ avcosφ - z_ avsinφ) - Ω̃_ R x_ av -γ z_ av .
    
    These equations have an analytical solution presented in <cit.> (note that the evolution of the y-coordinate is not important in our analysis). Thus we obtain the following correlators (alternative methods for the derivation are also discussed in <cit.>):
        
        K_zz(τ) = 1/2[ 1+Γ_z+cos(2φ) Γ_φ/Γ_+ -Γ_-]e^-Γ_-τ
             + 1/2[1 - Γ_z+cos(2φ) Γ_φ/Γ_+ -Γ_-] e^-Γ_+τ,
            
            K_zφ(τ) = (Γ_z + Γ_φ) cosφ  + 2Ω̃_ Rsinφ/2(Γ_+ - Γ_-)(e^-Γ_-τ -
        e^-Γ_+τ)
            
            +cosφ/2( e^-Γ_-τ + e^-Γ_+τ),
            
           Γ_± =Γ_z +Γ_φ±[ Γ_z^2 + Γ_φ^2 +2Γ_zΓ_φcos(2φ) -4 Ω̃_ R^2]^1/2/2
             +  (T_1^-1+T_2^-1)/2.
    
    Because of the rotational symmetry, the results for the correlators K_φφ(τ) and K_φ z(τ) can be obtained from Eqs. (<ref>) and (<ref>) by exchanging Γ_z↔Γ_φ and φ→ -φ. The rotational symmetry also makes the correlators insensitive to a y-rotation in both measurement directions, z→φ_ add, φ→φ+φ_ add, by any angle φ_ add.
    
    
    We emphasize that the obtained correlators do not depend on the qubit state ρ (t_1) and therefore on ρ_ in and t_ in (this property would not hold in the presence of phase backaction). We also emphasize that the correlators depend on Γ_z and Γ_φ, but do not depend on τ_z and τ_φ and therefore on the quantum efficiencies η_z and η_φ. Physically, this is because non-ideal detectors can be thought of as ideal detectors with extra noise at the output <cit.>, which only affects the zero-time self-correlators K_ii(0).
    
    Let us discuss some special cases for the results (<ref>)–(<ref>). (i) At small times, τ→ +0, we obtain correlators
    
        K_zz(+0)= 1,    K_zφ(0) = K_φ z (0)= cosφ.
    
    (ii) For |φ|≪ 1 and sufficiently small T_2^-1 and Ω̃_ R, we have Zeno pinning near the states |0⟩ and |1⟩ with rare jumps between them with equal rates Γ_ jump. This produces cross-correlator <cit.>K_zφ(τ)≈exp(-2Γ_ jumpτ) with jump rates
        Γ_ jump = φ^2Γ_zΓ_φ + Ω̃_ R^2/2(Γ_z+Γ_φ )  +  (T_1^-1+T_2^-1)/4 .
        
    (iii) In the case Ω̃_ R=T_1^-1=T_2^-1=0, we have full correlation for φ=0,  K_zφ(τ) = K_zz(τ ) = 1, full anticorrelation for φ=π, K_zφ(τ) = -K_zz(τ ) = -1, and no correlation for φ=π/2, K_zφ(τ) = 0, while K_zz(τ )=e^-Γ_φτ and K_φφ(τ )=e^-Γ_z τ. (iv) In the case Ω̃_ R=0, the cross-correlator is symmetric, K_zφ(τ )=K_φ z(τ), for any φ.
    
    
    
     Comparison with experimental results.—Experimental data have been taken in the same way as in Ref. <cit.> (see also <cit.>). Experimental  parameters correspond to well-separated frequency scales, as needed for the theoretical results, (T_1^-1,T_2^-1,|Ω̃_ R|)≪ (Γ_z, Γ_φ) ≪ (κ_z,κ_φ)≪Ω_ R, with T_1=60  μs, T_2=30 μs, Γ_z^-1=Γ_φ^-1= 1.3 μ s, damping rates of the two measurement resonator modes  κ_z/2π=4.3 MHz and κ_φ/2π=7.2 MHz, and Ω_ R≈Ω_ rf= 2π× 40 MHz.
    For this work we use 11 values for the angle φ between the Bloch-sphere directions of simultaneously measured qubit observables: φ_n=nπ/10, with integer n between 0 and 10. While φ_n is determined by well-controlled phases of applied microwaves <cit.>, the effective φ includes a small correction δφ =(κ_φ -κ_z)/2Ω_ R≈ 0.036  (see <cit.>), so that φ =φ_n+δφ. We have used about 200,000 traces per angle for the output signals Ĩ_z(t) and Ĩ_φ (t), each with 5 μs duration and 4 ns sampling interval.
    The traces are selected by heralding the ground state of the qubit at the start of a run and checking that the transmon qubit is still within the two-level subspace after the run. The recorded signals Ĩ_i(t) are linearly related to the normalized signals I_i(t) in Eqs. (<ref>) and (<ref>) as  Ĩ_i(t)=(ΔĨ_i/2)  I_i(t)+Ĩ_i^ off, where responses ΔĨ_i have been calibrated using ensemble-averaged ⟨Ĩ_i(t)⟩ (see details in <cit.>), giving in arbitrary units ΔĨ_z=4.0 and ΔĨ_φ=4.4. The offsets Ĩ_i^ off are approximately zeroed individually for each trace by measuring the non-rotating physical qubit after each run. Additional offset removal, |Ĩ_i^ off|≈ 0.15-0.20, for all traces with the same φ is done using ⟨Ĩ_i(t)⟩<cit.>. For calculating the correlators, we average over the ensemble of ∼200,000 traces and additionally average over time t_1 in Eq. (<ref>) within the 0.5 μs range 1 μ s≤ t_1 ≤ 1.5 μ s (first 1 μs is not used to avoid transients in the experimental procedure, and longer averaging reduces the range for τ; we also used averaging over 1 μs duration with similar results). Note that in the experiment the applied microwave phases in the two measurement channels actually correspond to angles ±φ_n/2; however, because of rotational symmetry, we still label the first measured operator as σ_z and the second operator as σ_φ. Also note that we use subscripts z and φ in various notations (Ĩ_i, κ_i, etc.) simply to distinguish the first (“z”) and second (“φ”) measurement channels.
    
    
    Figure <ref>(a) shows the agreement between the theory and the experimental data, where the solid lines show the symmetrized cross-correlator [K_zφ(τ) +K_φ z(τ)]/2 calculated from the experimental traces for 11 values of the angle φ, while the dashed lines correspond to the theoretical result, Eq. (<ref>).  For the analytics we used Ω̃_ R=0; however, there is practically no dependence on Ω̃_ R for the symmetrized cross-correlator, since the dependence comes only via Eq. (<ref>). Note that because of the Markovian assumption, our theory is formally valid only for τκ_i^-1∼ 30 ns; however, the experimental results agree with the theory even at τ < κ_i^-1. Figure <ref>(b) shows the same symmetrized cross-correlator at τ =0 as a function of φ. The agreement between the theory (cosφ, line) and the experiment (crosses) is also very good.
    
    The self-correlator K_zz(τ) as a function of τ is shown in Fig. <ref>(c) for 11 values of φ (results for K_φφ are similar).
    The agreement between the theory (dashed lines) and experiment (solid lines) is good, except for small τ (discussed below). A significant discrepancy for values of φ close to π/2 is probably caused by remaining offsets Ĩ_i^ off, which vary from trace to trace. Note that the lines in <ref>(c) come in pairs, corresponding to angles φ_n and π-φ_n. The separation of the analytical lines in the pairs is due to δφ, while separation of experimental lines is smaller, probably indicating a smaller value of δφ (partial compensation could be due to imperfect phase matching of applied microwaves or their dispersion in the cable).
    
    Looking at the experimental self-correlators K_zz(τ) and K_φφ(τ) at small τ for φ_n=π/2 [Fig. <ref>(d)], we see that in contrast to the theoretical results, there is a very significant increase of K_ii (τ) at τ 0.1  μ s. The discrepancy is due to the assumption of delta-correlated noise in our theory, while in the experiment the amplifying chain has a finite bandwidth (the Josephson parametric amplifiers have a half-bandwidth of 3.6 MHz and 10 MHz for σ_z and σ_φ channels, respectively), and the output signals Ĩ_i(t) are also passed through analog filters with a quite sharp cutoff at ∼25 MHz (this cutoff produces clearly visible oscillations with ∼40 ns period). Therefore, the theoretical delta-function contribution τ_i δ(τ) to K_ii(τ) becomes widened in experiment. It is interesting to note that, somewhat counterintuitively, a finite bandwidth of measurement resonator modes does not produce a contribution to K_ii(τ) at 0<τκ_i^-1 when Γ_i≪κ_i<cit.> (κ_z^-1≈ 37 ns, κ_φ^-1≈ 22 ns). This can be understood by considering a resonator without a qubit; then a finite bandwidth κ_i does not affect the amplified delta-correlated vacuum noise, so that only classical fluctuations of the resonator field (e.g., due to parameter fluctuations or elevated resonator temperature) will produce output fluctuations with 2/κ_i time scale. We have checked that the lines in Fig. <ref>(d) do not contain noticeable exponential contributions with decay time of 2/κ_i (small expected contributions with amplitude on the order of Γ_j/κ_i<cit.> are below experimental accuracy).
    
    
     Estimation of residual Ω̃_ R.—We now show that the antisymmetrized cross-correlator is a useful tool and can be used to estimate small residual Rabi oscillations frequency Ω̃_ R in the experiment. From Eq. (<ref>) we find
        
        K_zφ(τ)-K_φ z(τ )=2Ω̃_ Rsinφ/Γ_+-Γ_-(e^-Γ_-τ - e^-Γ_+τ) .
    
    Since in the case |Ω̃_ R|≪Γ_z,φ we can neglect Ω̃_ R in Eq. (<ref>) for Γ_±, Eq. (<ref>) gives a direct way to find Ω̃_ R from the experimental antisymmetrized cross-correlator. The solid line in Fig. <ref> shows K_zφ(τ)-K_φ z(τ ) from the experimental data for φ=π/2. Fitting this dependence on τ with Eq. (<ref>) (dashed line), we find the value Ω̃_ R/2π≈ 12 kHz, which is within the experimentally expected range of frequency mismatch between Ω_ R and Ω_ rf. Note that the overall shapes of the solid and dashed lines agree well with each other. Estimation of Ω̃_ R via the antisymmetrized cross-correlation is a very sensitive method and can be used to further reduce |Ω̃_ R| in an experiment, in which a direct measurement of 40 MHz Rabi oscillations with a few-kHz accuracy is a difficult task.
    
     Conclusion.—Using the quantum Bayesian theory for a simultaneous measurement of non-commuting qubit observables, we obtained analytical results for the self- and cross-correlators of the output signals from the measurement. Their comparison with experimental results shows a very good agreement. The correlators can be used for sensitive parameter estimation, in particular, to estimate and eliminate the mismatch between the Rabi oscillations and the sideband frequency shift used for measurement.
    
     Acknowledgements.—We thank Justin Dressel and Andrew Jordan for useful discussions. The work was supported by ARO grant No. W911NF-15-0496. L.S.M acknowledges support from the National Science Foundation Graduate Fellowship Grant No. 1106400.
    
    
    
    
    99Kraus1983 K. Kraus,  States, effects and operations: fundamental notoins of quantum theory (Springer, Berlin, 1983).
    
    Caves1986 C. M. Caves,  Phys. Rev. D 33, 1643 (1986).
    
    Menskii1998 M. B. Menskii, Phys. Usp. 41, 923 (1998).
    
    Belavkin1992 V. P. Belavkin, J. Multivariate Anal.  42, 171 (1992).
    
    Braginsky-book V. B. Braginsky and F. Ya. Khalili,  Quantum measurement (Cambridge University Press, Cambridge, UK, 1992).
    
    Aharonov1998 Y. Aharonov, D. Z. Albert, and L. Vaidman, Phys. Rev. Lett.
     60, 1351 (1988).
    
    Dalibard1992 J. Dalibard, Y. Castin, and K. Mø lmer, Phys. Rev. Lett.  68, 580 (1992).
    
    Carmichael1993 H. J. Carmichael,  An open systems
      approach to quantum optics (Springer, Berlin, 1993).
    
    Wiseman1993 H. M. Wiseman and G. J. Milburn, Phys. Rev. A 47, 642 (1993).
    
    Korotkov1999 A. N. Korotkov, Phys. Rev. B  60, 5737 (1999).
    
    Katz2006 N. Katz, M. Ansmann, R. C. Bialczak, E. Lucero, R. McDermott,
    M. Neeley, M. Steffen, E. M. Weig, A. N. Cleland, J. M. Martinis,
    and A. N. Korotkov, Science  312, 1498 (2006).
    
    Laloy2010 A. Palacios-Laloy, F. Mallet, F. Nguyen, P. Bertet, D. Vion, D.
    Esteve, and A. N. Korotkov, Nat. Phys.  6, 442 (2010).
    
    Vijay2012 R. Vijay, C. Macklin, D. H. Slichter, S. J. Weber, K. W. Murch,
    R. Naik, A. N. Korotkov, and I. Siddiqi, Nature (London)  490, 77 (2012).
    
    Hatridge2013 M. Hatridge, S. Shankar, M. Mirrahimi, F. Schackert, K.
    Geerlings, T. Brecht, K. M. Sliwa, B. Abdo, L. Frunzio, S. M. Girvin, R. J. Schoelkopf, and M. H. Devoret, Science  339, 178 (2013).
    
    Murch2013 K. W. Murch, S. J. Weber, C. Macklin, and I. Siddiqi, Nature
    (London)  502, 211 (2013).
    
    deLange2014 G. de Lange, D. Riste, M. J. Tiggelman, C. Eichler, L. Tornberg,
    G. Johansson, A. Wallraff, R. N. Schouten, and L. DiCarlo, Phys. Rev. Lett.  112, 080501 (2014).
    
    Campagne2014 P. Campagne-Ibarcq, L. Bretheau, E. Flurin, A. Auffeves, F.
    Mallet, and B. Huard, Phys. Rev. Lett.  112, 180402 (2014).
    
    Wiseman1993fb H. M. Wiseman and G. J. Milburn, Phys. Rev. Lett.  70,
    548 (1993).
    
    Ruskov2002  R. Ruskov and A. N. Korotkov, Phys. Rev. B  66, 041401
    (2002).
    
    Sayrin2011 C. Sayrin, I. Dotsenko, X. Zhou, B. Peaudecerf, T. Rybarczyk,
    G. Sebastien, P. Rouchon, M. Mirrahimi, H. Amini, M. Brune, J.-M. Raimond, and S. Haroche, Nature (London)  477, 73 (2011).
    
    Jacobs2003 K. Jacobs, Phys. Rev. A  67, 030301(R) (2003).
    
    Ruskov2003 R. Ruskov and A. N. Korotkov, Phys. Rev. B  67, 241305(R) (2003).
    
    Riste2013 D. Risté, M. Dukalski, C. A. Watson, G. de Lange, M. J.
    Tiggelman, Y. M. Blanter, K. W. Lehnert, R. N. Schouten, and L. DiCarlo, Nature (London)  502, 350 (2013).
    
    Roch2015 N. Roch, M. E. Schwartz, F. Motzoi, C. Macklin, R. Vijay, A. W. Eddins, A. N. Korotkov, K. B. Whaley, M. Sarovar, and I. Siddiqi, Phys. Rev. Lett.  112, 170501 (2014).
    
    Ahn2002  C. Ahn, A. C. Doherty, and A. J. Landahl, Phys. Rev. A  65, 042301 (2002).
    
    Sarovar2004 M. Sarovar, C. Ahn, K. Jacobs, and G. J. Milburn, Phys. Rev. A  69, 052324 (2004).
    
    Arthurs1965 E. Arthurs and J. L. Kelly, Bell System Technical Journal
     44, 725 (1965).
    
    Busch1985 P. Busch, Int. J. Theor. Phys. 24,  63 (1985).
    
    Stenholm1992 S. Stenholm, Ann. Phys. 218,  233 (1992).
    
    Jordan2005 A. N. Jordan and M. Büttiker, Phys. Rev. Lett.  95,
    220401 (2005).
    
    Wei2008 H. Wei and Yu. V. Nazarov, Phys. Rev. B  78, 045308 (2008).
    
    Ruskov2010 R. Ruskov, A. N. Korotkov, and K. Mø lmer, Phys. Rev. Lett.  105, 100506 (2010).
    
    Shay2016 S. Hacohen-Gourgy, L. S. Martin, E. Flurin, V. V. Ramasesh,
    K. B. Whaley, and I. Siddiqi, Nature (London)  538, 491 (2016).
    
    Averin2002 D. V. Averin, Phys. Rev. Lett.  88, 207901 (2002).
    
    Wiseman-book H. M. Wiseman and G. J. Milburn,  Quantum measurement and
    control (Cambridge Univ. Press, 2010).
    
    Gambetta2008 J. Gambetta, A. Blais, M. Boissonneault, A. A. Houck, D. I.
    Schuster, and S. M. Girvin, Phys. Rev. A  77, 012112 (2008).
    
    Korotkov2001-16 A. N. Korotkov, Phys. Rev. B  63, 115403 (2001); arXiv:1111.4016; Phys. Rev. A  94, 042326 (2016).
    
    
    supplement See Supplemental Material.
    
    Atalaya2016 J. Atalaya, M. Bahrami, L. P. Pryadko, and A. N. Korotkov,
    arXiv:1612.02096.
    
    Chantasri2015 A. Chantasri and A. N. Jordan, Phys. Rev. A 92, 032125 (2015).
    
    Korotkov2001sp A. N. Korotkov, Phys. Rev. B  63, 085312 (2001).
    
    Korotkov2011  A. N. Korotkov, Phys. Rev. B  83, 041406 (2011).
    
    
    
    
    
    
     Supplemental material for “Correlators in simultaneous measurement of non-commuting qubit observables”
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    § EXPERIMENTAL SETUP
    
    
    The experimental setup is the same as the one used in the experiment <cit.>, where full details can be found. For clarity we briefly describe the experimental apparatus for simultaneously applying and controlling two measurement observables. We use a transmon qubit placed inside an aluminum cavity, such that it is dispersively coupled to the two lowest modes of the cavity. The cavity has two outputs, each primarily coupled to a different mode. The outputs of these modes are amplified using two lumped-element Josephson parametric amplifiers (LJPA) operated in phase sensitive mode. Each mode is then used to measure an observable of the qubit, as described below. The apparatus is cooled to 30 mK inside a dilution refrigerator.
    
    We drive Rabi oscillations Ω_ R/2π = 40  MHz on the qubit by applying a resonant microwave tone modulated by an arbitrary waveform generator.
    
    In the frame rotating with Ω_ R, this produces an effective low frequency qubit. To couple the effective qubit to the cavity modes for measurement, we apply a pair of microwave sidebands to each mode. The sidebands are detuned above and below the two cavity modes by Ω_R, which leads to a resonant interaction between the qubit Rabi oscillations and the mode. This coupling may be understood as a stroboscopic measurement of the qubit oscillations. The relative phase of the sidebands determines which quadrature of the qubit oscillations is measured. This coupling causes the cavity mode state to displace in a way that depends on the state of the qubit.  We couple to the internal cavity field using a small antenna that protrudes into the cavity, allowing read out the cavity state as described above. Quantum trajectory reconstructions are validated using post-selection and tomographic measurements.
    
    
    
    
    
    § QUANTUM BAYESIAN APPROACH TO QUBIT MEASUREMENT IN RABI-ROTATED FRAME 
    
    
    In this section we develop the quantum Bayesian theory of the stroboscopic qubit measurement in the Rabi-rotated frame, used in the experiment <cit.> and briefly described above. We start with measurement of one effective observable σ_φ=σ_z cosφ +σ_x sinφ, then adding the second measurement in the same way and deriving Eqs. (4)–(6) of the main text. In this derivation we assume that the qubit Rabi frequency Ω_ R is exactly equal to the sideband frequency shift Ω_ rf (which defines the rotating frame), while a small mismatch between Ω_ R and Ω_ rf is added later via Eq. (7) of the main text (also discussed in this section). The focus is on the simple physics of the qubit measurement in the rotating frame.
    
    
    
     §.§ Measurement of one observable σ_φ
    
    
    The  physical qubit is Rabi-rotated with frequency Ω_ R about the y-axis, so that its Bloch coordinates rotate as
        
        x_ ph(t)= r_0 sin (Ω_ R t +ϕ_0),
            
                    y_ ph(t)=y_0,
            
               z_ ph(t) =r_0 cos (Ω_ R t+ϕ_0),
    
    where the radius r_0(t) within the xz-plane, the rotation phase ϕ_0(t), and the coordinate y_0(t) slowly change in time (e.g., due to measurement). The oscillations of the qubit z-component lead to a small change of the effective resonator frequency,
        ω_ r(t) = ω_ r^ m +χ  r_0 cos (Ω_ R t+ϕ_0),
        
    where  χ  is the (small) dispersive coupling between the qubit and the measurement resonator mode, and ω_ r^ m is the mean value between the resonator frequencies for the physical qubit states |0⟩ and |1⟩. Note that in this derivation, fast-oscillating ω_ r(t) is the value averaged over the physical qubit states, and we neglect quantum backaction developing during short time scale Ω_ R^-1.
    
    The sideband drive of the resonator at frequencies ω_ d±Ω_ rf with equal amplitudes (here for simplicity we assume ω_ d=ω_ r^ m and Ω_ rf=Ω_ R), produces the Hamiltonian term
        
        H_ drive/ħ =   εsin (Ω_ R t+φ)  a^† + h.c.,
        
    where ε is the normalized amplitude, φ depends on the initial phase shift between the sideband tones, a^† is the creation operator for the resonator, and we use the rotating frame based on ω_ d. The form of this term follows from the usual trigonometric relation for adding the sideband tones, (ε/2)sin[(ω_ d + Ω_ R)t+φ]-(ε/2)sin[(ω_ d - Ω_ R)t-φ]=εsin(Ω_ Rt+φ)cos(ω_ d t).
    
    
    This leads to the following dynamics of the resonator's coherent state  |α (t)⟩ [or classical field α (t)] in the rotating frame based on ω_ d,
        α̇ = -i χ r_0 cos(Ω_ R t+ϕ_0)  α - i εsin (Ω_ R t+φ ) -κ/2α,
        
    where we also took into account the resonator damping with energy decay rate κ.
    
    Now let us solve the evolution equation (<ref>), assuming κ≪Ω_ R and |χ| ≪Ω_ R. The drive term produces fast oscillation of α with Rabi frequency Ω_ R,
        Δα (t)= i ε/Ω_ Rcos (Ω_ R t +φ ).
        
    Inserting this oscillation into the first term of Eq. (<ref>), using the trigonometric formula cos(Ω_ R t+ϕ_0)cos (Ω_ R t+φ)= 1/2cos(ϕ_0-φ)+1/2cos (2Ω_ R t + ϕ_0 +φ ), and neglecting oscillations with frequency 2Ω_ R, we obtain the equation for the slow evolution,
        
        α̇_ s =   χε/2Ω_ R  r_0 cos(ϕ_0 -φ )  -κ/2α_ s ,
            
              α (t) =α_s (t)+Δα(t).
    
    Note that we can neglect the additional fast oscillations produced by the first term in Eq. (<ref>), Δ_2α (t)=-i(χ/Ω_R) r_0sin(Ω_ R t+ϕ_0) α_ s, in comparison with (<ref>), because we assume χ^2/(Ω_ Rκ)≪ 1.
    
    We see that the evolution (<ref>) of the resonator field α_ s depends on the state of the  effective qubit,
        
      x =r_0 sin(ϕ_0),     y=y_0,       z=r_0 cos (ϕ_0),
        
    which corresponds to the physical qubit (<ref>)–(<ref>) in the rotating frame Ω_ R. Moreover, this dependence is only on the Bloch φ-coordinate of the effective qubit, which is within the xz-plane at an angle φ from the z-axis, we see this since
        
        r_0 cos(ϕ_0 -φ ) =  Tr[σ_φρ (t)],
        
    where ρ(t) is the slowly-varying density matrix of the effective qubit.
    
    At this stage, we make use of the quantum Bayesian approach <cit.> to describe the qubit evolution due to measurement. Since the oscillating part Δα of the resonator field [Eq. (<ref>)] does not depend on the qubit state, it can be neglected in the analysis. In contrast, homodyne measurement of the leaked field α_ s
    gives us information on the value of the φ-coordinate of the effective qubit, which is a two-level system similar to the physical qubit. Inevitably, this information gradually collapses the effective qubit, i.e., changes its state according to the acquired information.
    
    
    The two σ_φ-basis states |1_φ⟩ and |0_φ⟩  of the effective qubit (σ_φ |1_φ⟩ =|1_φ⟩, σ_φ |0_φ⟩ =-|0_φ⟩) produce two steady states of the resonator, respectively (excluding oscillating Δα),
        α_ st,1 =  χε/Ω_ Rκ,    α_ st, 0=-α_ st,1.
        
    This is all what is needed in the Markovian “bad cavity” regime (when the evolution of the effective qubit is much slower than κ), which is assumed in the main text. Since in circuit QED only the difference between α_ st,1 and α_ st,0 is important for the analysis of the qubit evolution due to measurement in the “bad cavity” regime  <cit.>, the situation is  equivalent to the qubit evolution due to measurement in the standard setup <cit.> with the same  α_ st,1-α_ st,0. Correspondingly, the quantum Bayesian formalism in the “bad cavity” regime is exactly the same as for the standard circuit QED setup <cit.>, which coincides with the Bayesian formalism for qubit measurement using a quantum point contact <cit.>. The only difference is that now we discuss the evolution of the effective qubit instead of the physical qubit.
    
    In particular, when a phase-sensitive amplifier is used, the response ΔĨ of the detector to the effective qubit state has the dependence ΔĨ=ΔĨ_ maxcosθ on the phase difference θ between the amplified quadrature (at the microwave frequency ω_ d) and the optimal quadrature [which is real (horizontal), as follows from Eq. (<ref>)]. The informational backaction is proportional to cosθ, while phase backaction is proportional to sinθ, with ensemble dephasing of the effective qubit,
        Γ=κ/2  |α_ st,1-α_ st,0|^2,
        
    not depending on θ<cit.>. Evolution of the effective qubit state due to measurement is given by Eqs. (18) and (25) of Ref. <cit.> in the basis {|0_φ⟩, |1_φ⟩}  (Eqs. (12) and (13) of Ref. <cit.>). In this basis, the diagonal matrix elements ρ_00 and ρ_11 evolve due to the classical Bayes rule, while the off-diagonal elements ρ_01 and ρ_10 evolve due to evolving product ρ_00ρ_11 and also due to phase backaction.
    
    
    In the above derivation we assumed an exactly resonant microwave drive, ω_ d=ω_ r^ m. If this is not the case, |ω_ d-ω_ r^ m| ∼κ, then there will be an extra term -i (ω_ r^ m-ω_ d)α in Eq. (<ref>), which will lead to an extra term -i (ω_ r^ m-ω_ d)α_ s in Eq. (<ref>). Correspondingly, the steady states of the field α_ s for the effective qubit states |1_φ⟩ and |0_φ⟩  are
        α_ st,1 =-α_ st, 0 =  χε/Ω_ R[κ +2i(ω_ r^ m-ω_ d)] 
    instead of Eq. (<ref>), so that the optimal quadrature is no longer horizontal (real). The quantum Bayesian formalism remains the same.
    
    If the effective rotating-frame qubit is measured by only one detector (σ_φ, but no σ_z) and Ω_ R=Ω_ rf, then it is possible to go beyond the “bad cavity” limit and analyze transients within the time scale κ^-1. The derivation of the quantum Bayesian formalism for this case exactly follows the derivation in Ref. <cit.> and uses the field evolution equation (<ref>) instead of the steady-state solution (<ref>).
    
    
    
    
    
     §.§ Derivation of Eqs. (4)–(6) in the main text 
    
    
    In the absence of phase back-action, the quantum Bayesian equations describing continuous measurement of qubit σ_z observable in the Markovian approximation are <cit.>
        ẋ=    - τ_z^-1  xz   I_z(t)    - γ_z x,
            
        ẏ=     - τ_z^-1 yz   I_z(t) - γ_z y ,
            
        ż=     τ_z^-1(1-z^2)  I_z(t) ,
    
    in the Stratonovich form, where
        
        I_z(t)=  Tr[σ_z ρ (t)]+√(τ_z) ξ_z (t)
        
    is the normalized output signal, ξ_z(t) is the normalized white noise, ⟨ξ_z(t)  ξ_z(t')⟩ =δ (t-t'), τ_z is the “measurement” time after which the signal-to-noise ratio reaches 1, the qubit density matrix is ρ = ( +xσ_x+ yσ_y +zσ_z)/2, and qubit dephasing γ_z in individual measurement is related to the ensemble dephasing Γ_z as γ_z=Γ_z-(2τ_z)^-1.
    
    In the Itô form (i.e., using the forward definition of derivatives instead of the symmetric definition) these evolution equations become <cit.>
        ẋ=     - τ_z^-1/2xz  ξ_z  - Γ_z x  ,
            
        ẏ=     - τ_z^-1/2yz  ξ_z  - Γ_z y  ,
            
        ż=    τ_z^-1/2(1-z^2) ξ_z .
    
    
    
    When the observable σ_φ is measured instead of σ_z, these equations remain the same <cit.> in the basis of eigenstates |0_φ⟩ and |1_φ⟩, so that we can simply  change the notation: x→ x_φ, y→ y_φ, z→ z_φ. Rotating back to the usual basis, i.e., using the transformation x=x_φcosφ + z_φsinφ , y=y_φ, z=z_φcosφ -x_φsinφ, we obtain for the Itô form
    
        ẋ=     - τ_φ^-1/2[xzcosφ - (1-x^2)sinφ] ξ_φ ,
            
             - Γ_φcosφ(xcosφ - zsinφ) ,
            
        ẏ=     - τ_φ^-1/2 y [zcosφ + xsinφ] ξ_φ  - Γ_φ  y ,
            
        ż=    τ^-1/2_φ[(1 - z^2)cosφ -xzsinφ]ξ_φ
             + Γ_φsinφ(xcosφ - zsinφ) ,
    
    with
        
        I_φ(t)=  Tr[σ_φρ (t)]+√(τ_φ) ξ_φ (t).
        
    
    
    When both σ_z and σ_φ measurements are performed at the same time, we simply add the terms from Eqs. (<ref>)–(<ref>) and (<ref>)–(<ref>) <cit.> (with uncorrelated noises ξ_z and ξ_φ in the two channels), thus obtaining Eqs. (4)–(6) of the main text.
    
    
    
    
    
     §.§ Correction to measured rotation phase φ
    
    
    As was discussed above, the phase φ in the double-sideband drive εsin (Ω_ R t+φ)cos (ω_ dt) [see Eq. (<ref>)] directly determines the angle for the measured operator σ_φ for the effective qubit. This followed from the approximate solution of Eq. (<ref>). As we will see below, a more accurate solution shows a small correction to the measured direction φ.
    
    Neglecting the first term in Eq. (<ref>) but still keeping the last term, we obtain the (exact) oscillating solution
       Δα (t) = iε/Ω_ R+κ^2/4Ω_ R[ cos(Ω_ Rt+φ) -κ/2Ω_ Rsin (Ω_ Rt+φ) ] ,
        
    which is more accurate than Eq. (<ref>). (Note that the additional term ∝ e^-κ t/2 is naturally included in the slow dynamics.)
    Inserting Δα (t) into the first term of Eq. (<ref>), we obtain the following evolution of the slow part α_ s of the total field α =α_ s+Δα:
        
        α̇_ s =     χε/2( Ω_ R+κ^2/4Ω_ R)  r_0 [ cos(ϕ_0 -φ )-κ/2Ω_ Rsin (φ -ϕ_0) ]
            
            -κα_ s /2 ,
    
    which in the case κ/Ω_ R≪ 1 is approximately
            α̇_ s =   χε/2Ω_ R  r_0 cos(ϕ_0 -φ -κ/2Ω_ R )  -κ/2α_ s .
        
    This equation coincides with Eq. (<ref>), except φ is replaced with φ+κ/2Ω_ R, thus slightly changing the direction of  the measured operator for the effective qubit,
       
        r_0 cos(ϕ_0 -φ -κ/2Ω_ R ) =  Tr[σ_φ+κ/2Ω_ Rρ (t)].
        
    
    In the case when two channels simultaneously measure nominal operators σ_z and σ_φ (with corresponding resonator bandwidths κ_z and κ_φ), the measured directions on the Bloch sphere are actually κ_z/2Ω_ R and φ+κ_φ/2Ω_ R, so that the relative angle is φ+δφ with
        δφ = κ_φ-κ_z/2Ω_ R .
        
    This correction of ∼2^o was used in the main text when we compared the experimental results with the theory results.
    
    
    
    
     §.§ Decoherence of the effective qubit
    
    
    Now let us derive Eq. (7) of the main text, describing the evolution of the effective qubit not related to measurement.
    
    The evolution of the  physical qubit due to Rabi oscillations and environmental decoherence (energy relaxation and pure dephasing) is described by the standard master equation
    
    
        ρ̇_ ph = -i/ħ [H_ ph,ρ_ ph] + 1/2T_ pd ℒ[σ_z] + 1/T_1 ℒ[σ_-],     
              ℒ[A] ≡  Aρ_ phA^† - 1/2(A^†Aρ_ ph - ρ_ ph A^†A),
    
    where ρ_ ph(t) and H_ ph = ħΩ_ Rσ_y/2 are the density matrix and Hamiltonian of the physical qubit, respectively, the Lindblad-operator evolution with A=σ_z describes pure dephasing, and energy relaxation corresponds to  A=σ_- =|0⟩⟨ 1|= (σ_x - iσ_y)/2.
    
    To convert Eq. (<ref>) into the rotating frame (with frequency Ω_ rf), we apply the unitary transformation ρ(t) = U^†(t) ρ_ ph(t) U(t), where ρ(t) is the effective qubit density matrix and U(t) =exp(-iΩ_ rftσ_y/2). This gives
    
    
        ẋ   = Ω̃_ Rz - (1/2T_ pd + 3/4T_1)x + 1/T_1sin(Ω_ rft)
            
             + ( 1/4T_1 - 1/2T_ pd)[xcos(2Ω_ rft) + zsin(2Ω_ rft)],
            
        ẏ   = -(1/T_ pd + 1/2T_1)y,
            
        ż   = -Ω̃_ Rx - ( 1/2T_ pd + 3/4T_1) z - 1/T_1cos(Ω_ rft)
            
             + ( 1/4T_1 - 1/2T_ pd)[xsin(2Ω_ rft) - zcos(2Ω_ rft)],
    
    where Ω̃_ R = Ω_ R - Ω_ rf.
    Since Ω_ rf is much faster than the evolution of the effective qubit,  in Eqs. (<ref>) and (<ref>) we can neglect the oscillating terms. Finally expressing T_ pd via T_2 and T_1, so that T_ pd^-1+ (2T_1)^-1=T_2^-1 and (2T_ pd)^-1+3(4T_1)^-1=
    (T_1^-1+T_2^-1)/2, we obtain
    
    
        ẋ   = Ω̃_ Rz - 1/2(T_1^-1 + T_2^-1) x,
            
        ẏ   = T_2^-1y,
             
        ż   = -Ω̃_ Rx - 1/2(T_1^-1 + T_2^-1)z,
    
    which is Eq. (7) of the main text.
    
    Since Eqs. (<ref>)–(<ref>) describe only the evolution not related to measurement, while Eqs. (4)–(6) of the main text  describe only the evolution due to measurement, we need to simply add terms in these equations to describe the combined evolution of the effective qubit.
    We would like to emphasize that the derivation presented here relies on a significant separation of frequency scales
    (T_1^-1,T_2^-1,|Ω̃_ R|)≪(Γ_z, Γ_φ) ≪(κ_z,κ_φ)≪Ω_ R.
    
    In our experiment these inequalities are well satisfied:
    T_1=60 μs, T_2=30 μs, |Ω̃_ R^-1| 10  μs,
    Γ_z^-1=Γ_φ^-1=1.31 μ s, κ_z^-1= 37 ns,  κ_φ^-1=22.1 ns, and Ω_ R^-1=4 ns. The frequencies of the resonator modes should obviously be much larger than Ω_ R; in our experiment ω_ r, z/2π=7.4  GHz and ω_ r, φ/2π=6.7  GHz.
    
    
    
    
    
    § ANALYTICAL RESULTS FOR CORRELATORS
    
    
    
    We will first derive Eqs. (14)–(16) of the main text for the correlators using the “collapse recipe” and then we discuss the derivation of this recipe from the quantum Bayesian equations and its correspondence to the quantum regression approach.
    
    
    
    
     §.§ Derivation of Eqs. (14)–(16) in the main text using “collapse recipe”
    
    
    The collapse recipe for calculation of the correlators for the output signals (in the absence of phase backaction) was introduced in Ref. <cit.>. It says that in order to calculate the ensemble-averaged correlator,
    
        K_ij(τ) ≡⟨ I_j(t_1+τ)  I_i(t_1)⟩,     τ>0,
    
    we can replace the continuous measurement of σ_i at the earlier time moment t_1 with its projective measurement. It is also possible to replace the continuous measurement of σ_j at the later time moment t_1+τ with its projective measurement, but this is rather obvious and not important, since average values for the continuous and projective measurements coincide. In the following section we will show how this recipe can be derived from the quantum Bayesian equations (essentially repeating the derivation in Ref. <cit.>); here we just use this recipe.
    
    If the qubit state at time t_1 (it would be more accurate to say, right before t_1) is ρ (t_1), then the projective measurement of σ_i would produce the measurement result I_i(t_1)=1 with probability {1+ Tr[σ_i ρ (t_1)] }/2 and the result I_i(t_1)=-1 with probability {1- Tr[σ_i ρ (t_1)] }/2. After this projective measurement, the qubit state is collapsed to the eigenstate |1_i⟩ or |0_i⟩, correspondingly. Ensemble-averaged evolution after that is simple, since for an ensemble a continuous measurement is equivalent to decoherence. If the state was collapsed to |1_i⟩, then the further ensemble-averaged qubit evolution ρ_ av(t |1_i) starts with ρ_ av(t_1|1_i)=|1_i⟩⟨ 1_i|. Then the average result of σ_j measurement at time t=t_1+τ will be Tr[σ_j ρ_ av (t_1+τ|1_i)], which will produce contribution Tr[σ_j ρ_ av (t_1+τ|1_i)] × 1 to the correlator (<ref>) with probability {1+ Tr[σ_i ρ (t_1)] }/2.  Similarly, the contribution corresponding to the state collapse to |0_i⟩ at time t_1, is Tr[σ_j ρ_ av (t_1+τ|0_i)] × (-1) with probability {1- Tr[σ_i ρ (t_1)] }/2. Summing these two cases, we obtain
        
        K_ij(τ) =  Tr[σ_j  ρ_ av(t_1+τ|1_i)] 1+  Tr[σ_i ρ(t_1)]/2
           -   Tr[σ_j  ρ_ av(t_1+τ|0_i)]  1-  Tr[σ_i  ρ(t_1)]/2 ,
    
    which is Eq. (10) of the main text.
    
    Next, we need to find ρ_ av(t_1+τ|1_i) and ρ_ av(t_1+τ|0_i). The ensemble-averaged evolution of the qubit is described by Eqs. (11)–(13) of the main text. It is easy to see that the evolution of the y-coordinate is decoupled and has the simple solution,
        
      y_ av(t_1+τ) = e^-(Γ_z+Γ_φ)τ  y(t_1),
        
    where we use the subscript “av” to indicate ensemble averaging. The evolution equations for the coordinates x_ av and z_ av can be written in the matrix form,
        d/dtx_ avz_ av =  Mx_ avz_ av,
             M =  [ -(Γ_z +cos^2φΓ_φ+γ) sinφcosφ Γ_φ +Ω̃_ R; sinφcosφ Γ_φ -Ω̃_ R    -(sin^2φ Γ_φ +γ) ].
        
    Diagonalizing the matrix M, we find the solution
        
        x_ av(t_1+τ)z_
        av(t_1+τ) = [e^-Γ_-τ +e^-Γ_+τ/2  +
         e^-Γ_-τ - e^-Γ_+τ/2
              ×( Γ_φsin 2φ +2Ω̃_ R/Γ_+-Γ_- σ_x
          - Γ_z + Γ_φcos 2φ/Γ_+-Γ_- σ_z
         )
          ]x(t_1)z(t_1),
    
    where Γ_±  are the  eigenvalues of the matrix - M, given by Eq. (16) of the main text and repeated here,
        
        Γ_± =Γ_z +Γ_φ±[ Γ_z^2 + Γ_φ^2 +2Γ_zΓ_φcos(2φ) -4 Ω̃_ R^2]^1/2/2
             +  γ .
    
    
    Note that our evolution is symmetric under the inversion operation x→ -x,   y→ -y,   z→ -z (it is a unital map), and therefore
         Tr[σ_j  ρ_ av(t_1+τ|0_i)] = - Tr[σ_j  ρ_ av(t_1+τ|1_i)] .
        
    In this case Eq. (<ref>) simplifies to
        
     K_ij(τ) =  Tr[σ_j  ρ_ av(t_1+τ|1_i)] ,
        
    which no longer depends on the initial state ρ (t_1).
    
    
    In particular, if σ_i=σ_z, then to find K_ij(τ)  we need to use initial conditions x(t_1)=0 and z(t_1)=1 in Eq. (<ref>), which gives
       
        x_ av(t_1+τ |1) =
          sin(2φ) Γ_φ+2Ω̃_ R/2(Γ_+-Γ_-)( e^-Γ_-τ - e^-Γ_+τ) ,      
               z_ av(t_1+τ |1 ) = 1/2[ 1+Γ_z+cos(2φ) Γ_φ/Γ_+ -
        Γ_-]    e^-Γ_-τ
                 + 1/2[ 1 - Γ_z+cos(2φ) Γ_φ/Γ_+ -
        Γ_-]    e^-Γ_+τ .
    
    Inserting these results into Eq. (<ref>) and using relations Tr[σ_z  ρ_ av(t_1+τ|1_i)] =z_ av (t_1+τ) and Tr[σ_φ ρ_ av(t_1+τ|1_i)] =z_ av (t_1+τ)cosφ +x_ av (t_1+τ)sinφ, we obtain Eqs. (14) and (15) of the main text for the correlators K_zz(τ) and K_zφ(τ). If σ_i≠σ_z, then we can use rotational symmetry to find K_ij(τ), simply replacing φ with the angle difference between the measured directions and renaming the measurement channels.
    
    Note that the evolution (<ref>) for the y-coordinate was not important for K_zz(τ) and K_zφ(τ). Also note that we were able to use Eq. (<ref>) instead of Eq. (<ref>) because for the effective qubit the states |1⟩ and |0⟩ are equivalent (producing unital ensemble-averaged evolution). In similar calculations for a physical (non-rotating) qubit, energy relaxation would make states |1⟩ and |0⟩ non-equivalent, and then we would need to use Eq. (<ref>). Finally, we emphasize that this recipe is valid only in the absence of phase backaction. It requires a minor modification when phase backaction is present.
    
    
    
    
    
    
     §.§ Derivation via stochastic Bayesian equations
    
    
    Now let us derive Eqs. (14) and (15) of the main text for K_zz(τ) and K_zφ(τ ) using the
    stochastic evolution equations (4)–(6) of the main text instead of the collapse recipe used above. Even though equivalence of these methods was shown in Ref. <cit.>, we will do the derivation
    explicitly, essentially repeating the equivalence proof in
    <cit.>. In the derivation we assume fixed t_1
    and t_1+τ [averaging the correlator (<ref>) over the ensemble of realizations], and for brevity of notations we assume t_1=0. The qubit state
    right before the first measurement is therefore ρ_
    in≡ρ(0). (Note that if we have a distribution of the initial states,
    it is possible to average the correlator over this distribution
    later. However, such averaging is not actually needed because of the linearity of quantum evolution that allows us to use a single initial state, which is equal to the average over the distribution.)
    
    
    We will mainly consider K_zφ(τ ); the derivation for K_zz(τ ) is similar. Using Eqs. (<ref>) and (<ref>), we can write the correlator  K_zφ(τ ) as a sum of two parts, describing a correlation between qubit states at different times and a correlation between the noise
    and the future qubit state [there is no correlation with the past states because of causality, and the noise-noise correlations for τ > 0 are also absent for uncorrelated white noises ξ_z(t) and ξ_φ (t)],
        
        K_zφ(τ)= K^(1)_zφ(τ) +
        K^(2)_zφ(τ),
          
               K^(1)_zφ(τ) = ⟨ Tr[σ_φρ (τ)] ⟩  z(0) ,    
                K^(2)_zφ(τ) = √(τ_z) ⟨ Tr[σ_φρ (τ)]  ξ_z (0)⟩ ,
    
    where averaging is over the noise realizations ξ_z(t) and ξ_φ
    (t), which affect evolution of ρ via Eqs. (4)–(6) of the main text, and the initial state is ρ_ in=[ +x_ inσ_x + y_ inσ_y+z_ inσ_z ]/2 with { x(0),y(0),z(0)}={x_ in,y_ in,z_ in}.
    
    
    The first contribution can also be written as
        
    K^(1)_zφ(τ) =  Tr[σ_φρ_ av (τ|ρ_ in)]   z_ in,
        
    where  ρ_ av (τ|ρ_ in) is the ensemble-averaged density matrix at time τ, which starts with ρ_ in at t=0. Using linearity of the ρ_ av evolution given by Eqs. (<ref>) and (<ref>), we can  formally rewrite it as
        
    K^(1)_zφ(τ) =  Tr[σ_φρ_ av (τ|z_ inρ_ in)],
        
    where the evolution of ρ_ av now starts with state ρ_ av (0|z_ inρ_ in)=[ +z_ inx_ inσ_x + z_ iny_ inσ_y+z_ in^2σ_z ]/2. Note that in the definition of the state z_ inρ_ in we still use physical normalization Tr(z_ inρ_ in)=1, multiplying by z_ in only Bloch-sphere components of ρ_ in.
    
    
    To find the second contribution K^(2)_zφ(τ), we use the
    stochastic equations (4)–(6) of the main text [complemented with Eq. (7) of the main text] and derive the evolution equations for correlators ⟨ x(τ) ξ_z(0)⟩ and ⟨ z(τ) ξ_z(0)⟩:
       
        d/dτ⟨ x(τ) ξ_z(0)⟩⟨
        z(τ) ξ_z(0)⟩ =  M⟨ x(τ) ξ_z(0)⟩⟨ z(τ) ξ_z(0)⟩
               + 1/√(τ_z)-x_
        in z_ in1-z_ in^2 δ (τ),
    
    where M is the evolution matrix (<ref>), and for τ <0 these correlators are zero because of causality. This equation has a simple physical meaning. As follows from Eqs. (4) and (6) of the main text, the noise
    ξ_z(0) slightly changes the initial state after an infinitesimal
    time dt, so that x(0+dt)=x_ in-τ_z^-1/2x_ inz_
    inξ_z(0)   dt and z(0+dt)=z_ in+τ_z^-1/2(1-z_
    in^2) ξ_z(0)   dt. The further evolution starts with this
    slightly different state. Therefore, ⟨ x(dt) ξ_z(0)⟩ = -τ_z^-1/2x_ inz_ in⟨ξ_z(0)^2 dt⟩ =-τ_z^-1/2x_ inz_ in, since ⟨ξ_z(0)^2 dt⟩ = 1, as follows from Eq. (3) of the main text.
    Similarly, ⟨ z(dt) ξ_z(0)⟩ = τ_z^-1/2(1-z_ in^2)⟨ξ_z(0)^2 dt⟩ =τ_z^-1/2(1-z_ in^2). Thus we obtain the last
    term in Eq. (<ref>), while the evolution due to the
    matrix M is rather obvious. Even though the y component is not
    important for our analysis, for generality we can similarly derive ⟨
    y(dt) ξ_z(0)⟩ = -τ_z^-1/2y_ inz_ in⟨ξ_z(0)^2 dt⟩ =-τ_z^-1/2 y_ inz_ in.
    
    Since the evolution of  ⟨ x(τ) ξ_z(0)⟩ and ⟨ z(τ) ξ_z(0)⟩ is governed by the same matrix M as for the components of ρ_ av (similar for y-component), we can write the contribution K_zφ^(2)(τ) as
        
        K^(2)_zφ(τ)=  Tr[σ_φρ_ av(τ|δρ_ in)],
            
              δρ_ in = 1/2 [-x_ inz_ inσ_x
        -y_ in z_ inσ_y+(1-z_ in^2)σ_z ],
    
    where δρ_ in is an unphysical density matrix with zero trace, in which
    the Bloch-sphere components are the shifts discussed above due to the second term in Eq. (<ref>), multiplied by √(τ_z) because of Eq. (<ref>).
    
    It is easy to see that
        
        z_ inρ_ in + δρ_ in = [ +σ_z]/2=|1⟩⟨ 1 |
        
    (recall that z_ inρ_ in is defined with unity trace).
    Therefore, combining Eqs. (<ref>) and (<ref>), we find
        
        K_zφ(τ )=  Tr[σ_φρ_ av(τ|1)],
        
    which coincides with Eq. (<ref>) for σ_i=σ_z and σ_j=σ_φ. Note the slightly different notations for the initial state of ρ_ av, which should not be confusing, for example ρ_ av(τ|1)≡ρ_ av(τ| |1⟩⟨ 1|).
    
    Thus we have shown that the correlator K_zφ(τ ) derived from the stochastic evolution equations coincides with the result of the previous derivation based on the collapse recipe. The derivation for K_zz(τ) from the stochastic equations is similar, we just need to replace φ with z and σ_φ with σ_z in Eqs. (<ref>)–(<ref>), (<ref>), and (<ref>), thus obtaining K_zz(τ )=  Tr[σ_z ρ_ av(τ|1)], which coincides with Eq. (<ref>) for σ_i=σ_j=σ_z.
    
    
    
    
    
    
      §.§.§ Equivalence in a non-unital case
    
    
    We have shown equivalence of the results for correlators
    K_zz(τ) and K_zφ(τ) derived via the stochastic Bayesian equations and via the simple collapse recipe.
    However, in showing the equivalence we implicitly used
    the fact that the ensemble-averaged equations [Eqs. (11)-(13) of the main text] are homogeneous (not only linear). This is the so-called unital evolution (which preserves the center of the Bloch sphere), which originates from full symmetry between the states |0⟩ and |1⟩ of the effective qubit. Let us now prove that even in a non-unital case (for example, when we measure a physical qubit and asymmetry between states |0⟩ and |1⟩ is created by energy relaxation), the two methods for calculation of correlators are still equivalent (in the absence of phase-back-action). We will see that in this general case we can use Eq. (<ref>) originating from the collapse recipe, but cannot use its simplified version (<ref>).
    
    
    Let us use the linearity of the ensemble-averaged quantum evolution E (a trace-preserving positive map) from  t=0 to t=τ,
        ρ_ in→ E(ρ_ in)=ρ̃_c + x_ inΔρ̃_x +  y_ inΔρ̃_y +  z_ inΔρ̃_z,
        
    where ρ̃_c = E(ρ_c) is the state mapped from the Bloch sphere center ρ_c= /2, while Δρ̃_x = E(ρ_x)-ρ̃_c, Δρ̃_y= E(ρ_y)-ρ̃_c, and Δρ̃_z= E(ρ_z) -ρ̃_c describe mapping of the Bloch sphere axes, with density matrices ρ_x, ρ_y, and ρ_z corresponding to pure states (|1⟩ +|0⟩)/√(2), (|1⟩ +i |0⟩)/√(2), and |1⟩, respectively. Following the same logic as above, we can write the first contribution (<ref>) to K_zφ(τ) as
        
        K_zφ^(1)(τ) =
            [ Tr(σ_φρ̃_c) +x_ in Tr(σ_φΔρ̃_x) +y_ in Tr(σ_φΔρ̃_y)
            
                  +z_ in Tr(σ_φΔρ̃_z)]  z_ in.
    
    The second contribution (<ref>) to K_zφ(τ) can be written as
       
        K_zφ^(2)(τ) =
               -x_ inz_ in Tr(σ_φΔρ̃_x)   -y_ inz_ in Tr(σ_φΔρ̃_y)
               
                +(1-z_ in^2)   Tr(σ_φΔρ̃_z).
    
    Combining the two contributions, we find
        
          K_zφ(τ) =
             Tr(σ_φρ̃_c)   z_ in +  Tr(σ_φΔρ̃_z).
        
    
    On the other hand, using Eq. (<ref>) of the collapse recipe with σ_i=σ_z and σ_j=σ_φ , we obtain
        
        K_zφ(τ) = [  Tr(σ_φρ̃_c)   +  Tr(σ_φΔρ̃_z)] 1+z_ in/2
               -  [ Tr(σ_φρ̃_c)  -  Tr(σ_φΔρ̃_z)] 1-z_ in/2,
    
    which coincides with Eq. (<ref>). Thus, equivalence of both methods for K_zφ(τ) is proven in the general (non-unital) case. The proof for the correlator  K_zz(τ) is practically the same, just replacing σ_φ with σ_z. The proof of the equivalence for an arbitrary K_ij(τ) is also similar, but we need to use the basis, corresponding to σ_i.
    
    
    Note that in the non-unital case we should use Eq. (<ref>) of the collapse recipe, which takes into account both scenarios (collapse to the state |0⟩ or to |1⟩) and not the simplified equation (<ref>) (collapse to |1⟩ only), which is valid only for a symmetric (unital) evolution. As seen from Eq. (<ref>) [or Eq. (<ref>)], in the non-unital case the correlator K_ij(τ) depends on the initial state ρ_ in via the term Tr [σ_j ρ_ av (t_1+τ |ρ_ c)]  Tr [σ_i ρ (t_1)], where ρ_ c is the fully mixed state.
    
    
    Also note that in our derivation we assumed absence of the phase backaction terms <cit.> in the Bayesian stochastic equations. These terms would introduce additional contribution to δρ_ in and therefore to correlators. The collapse recipe in this case should be modified accordingly.
    
    
    
    
    
     §.§ Derivation via quantum regression approach
    
    
    Now let us derive Eqs. (14) and (15) of the main text for the
    correlators K_zz(τ) and K_zφ(τ) using the
    standard non-stochastic approach <cit.>, which cannot describe
    individual realizations of the qubit measurement process, but is
    sufficient to calculate correlators. In this section we
    assume t_1=0 and use ρ_ in=ρ(0).
    
    In this approach <cit.> we need to use the Heisenberg picture and associate
    the measurement outcomes I_z(t) and I_φ (t) with the
    operators σ_z(t) and σ_φ (t), which evolve in
    time as σ_z (τ) ≡ e^i H_ totτσ_z
    e^-iH_ totτ and σ_φ (τ) ≡ e^i H_
    totτσ_φ e^-iH_ totτ, where H_ tot
    is the total Hamiltonian describing the qubit, environment, and
    interaction between them (in this approach we consider the measurement
    apparatus as an environment). The correlators of the outcomes then can
    be expressed as symmetrized combinations
        
        K_zz(τ) = 1/2 Tr_tot [
        σ_z(τ) σ_z(0) ρ_ tot(0)
            
              
         + σ_z(0) σ_z(τ)  ρ_ tot(0)]
            
               =  Re{ Tr_tot [ σ_z (τ)  σ_z(0)  ρ_ tot(0) ]},
            
               K_zφ(τ) = 1/2 Tr_tot[ σ_φ(τ) σ_z(0) ρ_ tot(0)
            
               + σ_z(0) σ_φ(τ)
         ρ_ tot(0)]
            
               =  Re{ Tr_tot [ σ_φ
        (τ)  σ_z(0)  ρ_ tot(0) ]},
    
    where ρ_ tot(0)=ρ (0) ⊗ρ_ env(0) is the
    initial density matrix, which includes the environment (“bath”),
    and the trace should be taken over the qubit and environment degrees
    of freedom.
    
    We need to assume that the coupling between the qubit and the
    environment is sufficiently weak, so that the effective decoherence
    rate of the qubit due to its coupling with the environment is much
    smaller than the reciprocal of the typical correlation time for the
    bath degrees of freedom (this includes the “bad cavity” assumption). In this case we can use the standard formula <cit.>
    (related to what is usually called the Quantum Regression Theorem)
         Tr_ tot [ A (τ)  B(0)  ρ_ tot
    (0)] =  Tr _ sys [ A  ρ_ av(τ|Bρ_
    in)] ,
        
    where in the right-hand side the trace is only over the system (qubit),
    operators A and B are system observables, and ρ_ av(τ|Bρ_ in) is the system (reduced) density matrix at
    time τ, which evolves in time according to the
    ensemble-averaged (reduced) evolution equations and starts in the
    state Bρ_ in, i.e., ρ_ av(0|Bρ_ in)=Bρ_ in.
    Note that ρ_ av(τ|Bρ_ in) is unphysical
    because it starts with an unphysical initial state Bρ_ in (it is typically not Hermitian and not normalized). Also note that the
    validity of Eq. (<ref>) requires that the system and the
    environment will be weakly entangled; i.e., ρ_ tot(t)≈ Tr_ env[ρ_ tot(t)]⊗ρ_ env(0).
    This is consistent with the above assumption that the coupling is
    weak.
    
    In our case in Eq. (<ref>) the operator B is σ_z,
    while A is either σ_z or σ_φ. The starting
    state for ρ_ av(τ|Bρ_ in) is σ_zρ_ in, so
        
        K_zz(τ) =  Re{ Tr
            [ σ_z   ρ_ av(τ|σ_z ρ_ in) ]},
            
                K_zφ(τ) =  Re{ Tr
            [ σ_φρ_ av(τ|σ_z ρ_ in) ]}.
    
    Since we have to work with unphysical unnormalized states, we use ρ(t) = [P_ N
    + x(t) σ_x + y(t) σ_y + z(t) σ_z]/2,
    where the normalization is conserved, Ṗ_ N=0, during the ensemble-averaged evolution described by Eqs. (11)–(13) of the main text.
    
    Now let us represent the unphysical initial state σ_zρ_ in as
        
        σ_z [      ρ_11, in     ρ_10,  in; (ρ_10,  in)^*     ρ_00,  in ] =
           σ_z/2+(ρ_11, in-ρ_00, in)/2
              
        +    [             0      ρ_10, in; -(ρ_10, in)^*             0 ] ,
    
    using condition ρ_00, in+ρ_11, in=1. Since the
    ensemble-averaged evolution is linear,
    we can separate ρ_ av(τ|σ_zρ_ in) into
    three terms, corresponding to the three terms in Eq.(<ref>). The first term, σ_z/2, gives the physical evolution ρ_ av(τ |1) starting with the state |1⟩. The second term, (ρ_11, in -ρ_00, in)  /2, does not change in time and gives zero contribution to the correlators (<ref>) and (<ref>). The third term is initially anti-Hermitian, and it will remain anti-Hermitian
    in the evolution, because all coefficients in Eqs. (11)–(13) of the main text are real. The anti-Hermitian term will
    give zero contribution to Eqs. (<ref>) and (<ref>) because the traces will be
    imaginary numbers.
    
    Thus we obtain equations
        
    K_zz(τ) =  Tr
        [ σ_z   ρ_ av(τ|1)],    
        K_zφ(τ) =  Tr
        [ σ_φρ_ av(τ|1)],
        
    which coincide with Eq. (<ref>) for σ_i=σ_z. Therefore, the final result for the correlators is the same as for the derivation based on the collapse recipe.
    
    
    In the general non-unital case (without phase back-action), using representation (<ref>) of a linear quantum map, we can obtain Eq. (<ref>) from Eqs. (<ref>) and (<ref>), thus proving that the derivation via the quantum regression approach is still equivalent to the derivations via the stochastic equations and via the collapse recipe.
    
    
    
    
    
    
    
    
    
    § EXTRACTING CORRELATORS FROM EXPERIMENTAL DATA
    
    
    
    The experimental correlators are calculated as
    
    
        K_ij (τ) = ∫_t_a^t_bdt_1 ⟨(Ĩ_i(t_1)-Ĩ_i^ off)(Ĩ_j(t_1+τ)-Ĩ_j^ off)⟩/(t_b-t_a)  ΔĨ_i  ΔĨ_j,
    
    where Ĩ_i(t) are experimental output signals for σ_i measurement (0≤ t≤ 5 μs), ⟨ ...⟩ denotes ensemble averaging over all selected traces with the same angle difference φ (∼200,000 per angle, the selection includes heralding at the start of the run and checking that the physical transmon qubit is in the subspace of its lowest two energy levels after the run), additional time-averaging is between t_a=1 μs and t_b=1.5 μs, the correlators are normalized by responses ΔĨ_i, and small offsets Ĩ_i^ off are calculated separately for each value of φ (see below). Significantly larger offsets are already removed from Ĩ_i(t) individually for each trace by measuring and averaging the background noise for the non-rotating qubit after each trace.
    
    
    
    
    To find ΔĨ_i and Ĩ_i^ off, for each angle φ we separate the traces (each trace includes outputs for both measurement channels) into two approximately equal groups. These groups correspond to the effective qubit initialized either in the state |1⟩ (z_0=1) or |0⟩ (z_0=-1), which is controlled by the initial state |1⟩ or |0⟩ of the physical qubit before application of 40 MHz Rabi oscillations and stroboscopic sideband measurement. Calibration of the z-axis of the effective qubit is done by maximizing the response (for the lower-κ channel) for zero nominal angle between the two measurement directions (φ=0 neglecting small δφ). For non-zero nominal φ, the stroboscopic measurement directions are -φ/2 (channel 1, ω_ r,1/2π=7.4  GHz, κ_1/2π=4.3  MHz) and φ/2 (channel 2, ω_ r,2/2π=6.7  GHz, κ_2/2π=7.2  MHz) – see Fig. <ref>. Theoretically only the angle difference φ matters for the correlators; however, for calibration we need to use the actual measurement directions ±φ/2 (more accurately, -φ/2 and φ/2+δφ). In the main text the channel 1 is called z-channel, while channel 2 is φ-channel. In this section we will also be using the terminology of channels 1 and 2.
    
    
    
    
    
    
    
    To find the responses ΔĨ_1≡ΔĨ_z and  ΔĨ_2≡ΔĨ_φ for the two channels, we calculate 𝒟_i(t)≡⟨Ĩ_i(t)⟩_z_0=1 - ⟨Ĩ_i(t)⟩_z_0=-1 for each φ, where the subscripts z_0=± 1 denote the group of traces with initial state z_0 of the effective qubit. This quantity can be also obtained theoretically from the ensemble averaged evolution equations (11)–(13) of the main text with the initial condition x(0)=sin(φ/2), y(0)=0, and z(0)=cos(φ/2), see Fig. <ref>. It is equal to 𝒟_1(t)= ΔĨ_1 z_ av(t) and 𝒟_2(t) = ΔĨ_2 [z_ av(t)cos(φ+δφ) + x_ av(t)sin(φ+δφ)] for the first and second channels, respectively. In Fig. <ref>, we plot experimental 𝒟_1(t) and 𝒟_2(t) for 5 values of φ (φ_n=nπ/10, n=0,1,2,3,4), and fit them with theoretical results. We find a good agreement for the responses ΔĨ_1 = 4.0 and ΔĨ_2 = 4.4 in units of the experimental output. In this fitting we disregard any residual Rabi oscillations (Ω̃_ R=0).
    
    
    
    
    
    To estimate the offsets Ĩ_i^ off, we use the symmetric combination 𝒮(t) ≡ [⟨Ĩ_i(t)⟩_z_0=1 + ⟨Ĩ_i(t)⟩_z_0=-1]/2, which is shown in Fig. <ref> for both measurement channels and for 11 angles φ. We see that the 𝒮(t) are approximately independent of time, and therefore we can introduce the offsets Ĩ_i^ off=𝒮 for each channel and each value of φ. The offsets only weakly depend on the angle φ , but are significantly different in the two channels. For the first (z) channel we crudely find Ĩ_1^ off≡Ĩ_z^ off  = 0.15, 0.16, 0.16, 0.16, 0.17, 0.16, 0.16, 0.17, 0.16, 0.17 0.18 for the angles φ in increasing order. For the second (φ) channel we find Ĩ_2^ off≡Ĩ_φ^ off = -0.18, -0.17, -0.17, -0.15, -0.16, -0.17, -0.17, -0.17, -0.18, -0.19 and -0.19 for the angles φ in increasing order.
    
    
    § SELF-CORRELATORS AT SMALL Τ
    
    
    In this section we discuss why the amplified vacuum noise is still white (delta-correlated) for finite damping rate κ of a resonator. We also estimate the self-correlator contribution K_ii(τ) ∝exp (-κ_i τ/2), with small amplitude ∼Γ /κ_i due to qubit evolution.
    
    
    
    It is a somewhat surprising result that the self-correlator K_ii(τ) does not have a significant contribution ∝exp (-κ_i τ/2), originating from the correlation time (κ_i/2)^-1 of vacuum fluctuations inside the resonator. [We may naively expect that this would widen the contribution K_ii (τ)=η_iτ_iδ(τ) from the amplified vacuum noise.] To show why this is not the case, let us consider a resonator with finite κ without a qubit (Fig. <ref>) and calculate the correlator for the amplified vacuum noise, coming from the resonator. The coupling between the resonator and the output transmission line is κ_ out, while the remaining dissipation rate κ-κ_ out is modeled as a coupling to another transmission line.
    
    
    
    Using the standard input-output theory <cit.>, we need to consider the vacuum noise v̂(t) incident to the resonator from the output line, with the operator correlator ⟨v̂(t) v̂^†(t')⟩=δ (t-t'), and write the equations for the annihilation operators in the Heisenberg picture. However, for our purposes it is sufficient to use a simpler approach (e.g., Appendix B of <cit.>), in which we consider the evolution of classical fields (in the usual, i.e., Schrödinger picture) due to “classical” vacuum noise v(t) (complex number) with the correlator
             ⟨ v(t)  v^*(t')⟩ = 1/2 δ (t-t'),    ⟨ v(t)  v(t')⟩ =0,
        
    where real and imaginary parts of v(t) correspond to orthogonal quadratures. As follows from Eq. (<ref>), any quadrature has correlator (1/4) δ (t-t'), and orthogonal quadratures are uncorrelated.
    
    
    
    
    
    The evolution of the resonator field fluctuation α(t) (evolution of the field due to drive is decoupled due to linearity) is
        α̇=-i Δω α -κ/2 α
        +√(κ_ out)  v(t) +√(κ -κ_ out)  v_ a(t),
        
    where Δω=ω_ r-ω_ d is the resonator frequency ω_ r in the rotating frame based on the drive frequency ω_ d (it is important for homodyne detection) and v_ a(t) is the additional vacuum noise from the other transmission line (see Fig. <ref>) with the same correlator (<ref>) and uncorrelated with v(t). In Eq. (<ref>) we use the standard normalizations for the resonator field (based on the average number of photons) and for the propagating fields (based on average number of propagating photons per second). This equation has the simple solution,
        α (t) = ∫_-∞^t [√(κ_ out)  v(t')+
        √(κ -κ_ out)  v_ a(t')]   e^-κ̃(t-t')/2  dt',
        
    where κ̃= κ +2iΔω.
    
    The outgoing field F(t), which is then amplified is
        
        F(t) = -v(t)+√(κ_ out) α (t).
        
    After a phase-sensitive amplification, this field is sent to a homodyne detector, which outputs the signal I(t)∝ Re[F(t) e^-iθ], where θ is the amplified quadrature. Without loss of generality, we can assume θ=0 (by properly defining the quadrature). Therefore, we are interested in the self-correlator of Re[F(t)], which is equal (up to a coefficient) to the output signal correlator.
    
    Our goal is to show that the correlator of Re[F(t)] is the same as for vacuum noise, i.e.,
        
        K_ ReF(τ)≡⟨ Re[F(t)]   Re[F(t+τ)]⟩ = 1/4 δ(τ ).
        
    Actually, it is sufficient to show that K_ ReF(τ )=0 for τ >0, since the coefficient 1/4 in Eq. (<ref>) can be simply obtained from Eqs. (<ref>) and (<ref>). It is also sufficient to choose t=0 in Eq. (<ref>).
    
    As follows from Eq. (<ref>), the correlator K_ ReF(τ) for τ≥ 0 has three contributions, K_ ReF(τ)=K_ ReF^(1)(τ) +K_ ReF^(2)(τ)+K_ ReF^(3)(τ), where
        
        K_ ReF^(1)(τ) = ⟨ Re[v(0)]   Re[v(τ)]⟩  = 1/4 δ(τ),   
               K_ ReF^(2)(τ) = -√(κ_ out) ⟨ Re[v(0)]   Re[α(τ)]⟩  
               K_ ReF^(3)(τ) =κ_ out⟨ Re[α(0)]   Re[α(τ)]⟩ .
    
    Note that for τ≥ 0 there is no contribution due to correlation between α(0) and v (τ) because of causality. We need to show that the second and third contributions exactly cancel each other, K_ ReF^(2)(τ)+K_ ReF^(3)(τ)=0.
    
    Using Eq. (<ref>), correlator ⟨ Re[v(t)]  Re[v(t')]⟩ =(1/4) δ (t-t'), and absence of correlation between Re[v(t)] and Im[v(t')], we easily obtain (τ > 0)
            
    K_ ReF^(2)(τ) =-(κ_ out/4)  Re[e^-κ̃τ/2].
        
    Similarly (with a little more work) we obtain
        
        K_ ReF^(3)(τ) = κ_ outκ/4∫_-∞^0 [  Re( e^κ̃t'/2)   Re( e^-κ̃(τ - t')/2) .
            
              .
        +  Im( e^κ̃t'/2)   Im( e^-κ̃(τ - t')/2) ] dt'
            
              
        = (κ_ out/4)   e^-κτ /2cos(Δω τ) ,
    
    where addition of contributions from v(t) and v_ a(t) gives the coefficient κ=κ_ out+(κ -κ_ out) on the first line. It is easy to see that Eqs. (<ref>) and (<ref>) exactly cancel each other, thus proving Eq. (<ref>). The proof using the standard input-output theory <cit.> is essentially the same as our proof, just using operators instead of complex numbers and associating time-dependence with the Heisenberg picture.
    
    This result explains why we do not see significant exponential contributions ∝exp(-κ_iτ /2) to the self-correlators K_ii(τ ) at small τ in Fig. 2(d) of the main text. However, finite bandwidth of the amplifier leads to widening of the delta-function correlator of the amplified signal, producing crudely exponential time dependence at small τ in Fig. 2(d) of the main text.
    
    While we have shown above that in the ideal case the amplified noise is delta-correlated, small contributions ∝ exp(-κ_iτ /2) to the self-correlators are possible due to various non-idealities. For example, if the temperature of the resonator is significant, then the delta-correlator of the noise v_ a(t) is larger than the vacuum correlator of v(t). Repeating the derivation above, we see that the coefficient in Eq. (<ref>) increases, and the cancellation by Eq. (<ref>) is incomplete. Similarly, if the amplitude of the microwave drive fluctuates in time, these fluctuations are essentially passed through a filter with the resonator bandwidth, creating a contribution  ∝exp(-κ_iτ /2) from the “white-noise” part of the fluctuations.
    
    A similar mechanism is produced by random evolution of the qubit, which is much slower than κ_i, but still has a non-zero spectral weight at frequencies comparable to κ_i. Let us estimate the corresponding contribution to the self-correlator K_ii(τ) in the following way. Assuming φ=π/2 (so that z and x components of the qubit are measured) and assuming Γ_z=Γ_x=Γ, let us consider the uniform diffusion of the qubit state along the x-z great circle on the  Bloch sphere <cit.> with the angular diffusion coefficient 2Γ. Note that we assume an ideal detector by separating a non-ideal detector into an ideal part and extra noise. Even though the Markovian theory <cit.> cannot describe the qubit evolution at the frequency scale κ_i, in this estimate we just assume the same uniform diffusion with coefficient 2Γ.
    
    In the Markovian approximation, the qubit evolution characterized by the angle β(t) from the x-axis, produces the output signal in the σ_z-channel I_z(t)=sinβ (t) (excluding noise).  However, because of the finite bandwidth κ_z of the resonator, there will be a correction -∫_-∞^t e^-κ_z(t-t')/2cosβ β̇  dt' to the output signal due to transient delay. The contribution from this correction to the self-correlator K_zz(τ)=T^-1∫_-T/2^T/2 I_z(t) I_z(t+τ)  dt (with T→∞) is κ_z^-1e^-κ_zτ/2⟨ [ (cosβ) Δβ ]^2/Δ t ⟩, where Δβ=β̇ Δ t is the change during small Δ t. Using ⟨cos^2β⟩ =1/2 and ⟨ (Δβ)^2⟩ =2ΓΔ t, we obtain the contribution (Γ/κ_z) e^-κ_zτ/2 to K_zz(τ). A slightly more accurate calculation, which takes into account qubit diffusion during resonator transients, produces the result (Γ/κ_z) e^-κ_zτ/2e^-Γτ, which is practically the same since in our case Γ≪κ_z. The same crude derivation can be done for the σ_x-channel. Thus, for φ=π/2 we expect the contribution ∼ (Γ /κ_i)  e^-κ_iτ/2 to the self-correlator K_ii(τ ) at small τ. However, these contributions are not visible in Fig. 2(d) of the main text because experimentally Γ /κ_z=0.028 and Γ /κ_φ=0.017, which is almost three orders of magnitude less than the scale of the amplifier-caused effect contributing to the lines in Fig. 2(d) [note that by design of the experiment Γ/κ_i≪ 1– see Eq. (<ref>)].
    
    
    
    
    9S-Shay2016 S. Hacohen-Gourgy, L. S. Martin, E. Flurin, V. V. Ramasesh,
    K. B. Whaley, and I. Siddiqi, Nature (London)  538, 491 (2016).
    
    S-Korotkov-99-01 A. N. Korotkov, A. N. Korotkov, Phys. Rev. B  60, 5737 (1999); Phys. Rev. B  63, 115403 (2001).
    
    S-Korotkov-2011 A. N. Korotkov, arXiv:1111.4016.
    
    S-Korotkov-2016 A. N. Korotkov, Phys. Rev. A  94, 042326 (2016).
    
    S-Wiseman1993 H. M. Wiseman and G. J. Milburn, Phys. Rev. A 47, 642 (1993).
    
    S-Gambetta2008 J. Gambetta, A. Blais, M. Boissonneault, A. A. Houck, D. I.
    Schuster, and S. M. Girvin, Phys. Rev. A  77, 012112 (2008).
    
    S-Blais-2004 A. Blais, R.-S. Huang, A. Wallraff, S. M. Girvin, and R. J.
    Schoelkopf, Phys. Rev. A  69, 062320 (2004).
    
    S-Ruskov2010 R. Ruskov, A. N. Korotkov, and K. Mø lmer, Phys. Rev. Lett.  105, 100506 (2010).
    
    S-Korotkov2001sp A. N. Korotkov, Phys. Rev. B  63, 085312 (2001).
    
    S-Gardiner-book C. W. Gardiner and P. Zoller,  Quantum noise (Springer, Berlin, 2004), Sec. 5.2.
    
    S-Gardiner-1985 C. W. Gardiner and M. J. Collett, Phys. Rev. A  31, 3761 (1985).
    
    S-Clerk-2010 A. A. Clerk, M. H. Devoret, S. M. Girvin, F. Marquardt, and R. J. Schoelkopf, Rev. Mod. Phys.  82, 1155 (2010).
