Correlators in simultaneous measurement of non-commuting qubit observables

We consider the simultaneous and continuous measurement of qubit observables $\sigma_z$ and $\sigma_z\cos\varphi + \sigma_x\sin\varphi$, 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.

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 [27][28][29][30]. For the measurement of non-commuting observables of a qubit, statistics of timeintegrated detector outputs and fidelity of state monitoring directly via time-integrated outputs has been analyzed in Ref. [31]. The evolution of the qubit state due to simultaneous measurement of incompatible variables has been described theoretically in Ref. [32], and has been recently demonstrated experimentally in Ref. [33].
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 σ ϕ ≡ σ z cos ϕ + σ x sin ϕ, where σ x and σ z are the Pauli matrices and ϕ is an angle between the two measurement directions on the Bloch sphere (Fig. 1). The experimental setup is described in detail in Ref. [33]; it is based on a Rabi-rotated physical qubit, which is measured stroboscopically [34] 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 We consider the simultaneous continuous measurement of qubit observables σz and σϕ, which differ by an angle ϕ on the Bloch sphere, and calculate time-correlators for the output signals Iz(t) and Iϕ(t) resulting from this measurement.
quantum trajectories [8,9,35,36] has been developed in Ref. [33]. In this letter we will use a simpler approach based on quantum Bayesian theory [10,37]. The quantum Bayesian description of the rotating-frame experiment [33] is developed in the Supplemental Material [38]. The goal of this letter is calculation of the timecorrelators 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 [39]. We note that the analyzed output signal correlators are different from qubit-state correlators [40].
Quantum Bayesian theory.-A simultaneous continuous measurement of the qubit observables σ z and σ ϕ by two linear detectors ( Fig. 1) produces noisy output signals I z (t) and I ϕ (t), respectively [32,37], 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 and ξ z (t) ξ ϕ (t ) = 0. The qubit state is characterized in the Bloch-sphere representation as ρ(t) ≡ [1 1 + x(t) σ x + y(t) σ y + z(t) σ z ] /2. We assume phasesensitive 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 [35][36][37]. Then there is only the quantum informational backaction, which for measurement of σ z and σ ϕ is described [32,37,38] by the evolution equations (in the Itô interpretation) Here Γ z and Γ ϕ are the ensemble dephasing rates due to measurement, so that the quantum efficiencies [37] for the two channels are η z = 1/(2τ z Γ z ) and η ϕ = 1/(2τ ϕ Γ ϕ ). In the experiment η z = 0.49 and η ϕ = 0.41. Equations (4)-(6) 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 frequencỹ Ω R . In the experiment,Ω R = Ω R − Ω rf is a small (kHzrange) undesired mismatch between the physical qubit Rabi frequency Ω R and rotating frame frequency Ω rf defined by detuning of sideband pumps [33,38]. 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 Evolution of the effective qubit is described by adding terms from Eqs. (4)-(6) and (7).
Correlators.-Our next goal is to calculate the twotime correlators, K ij (τ ), for the output signals, Self-and cross-correlators correspond to i = j and i = j, respectively. The averaging in Eq. (9) 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. (9) can also be understood as averaging over time t 1 . We assume that the parameters in Eqs. (4)-(7) 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 [38], calculation of the correlators from Eqs. (1)-(7) is equivalent to the following recipe [41]: 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 eigen- This gives the correlator where ρ av (t 1 +τ |1 i ) is the ensemble-averaged density matrix at time t 1 + τ with the initial condition ρ av ( The evolution of ρ av is given by Eqs. (4)-(7) without noise, ξ z = ξ ϕ = 0 (because of the Itô form), so thaṫ These equations have an analytical solution presented in [38] (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 [38]): Because of the rotational symmetry, the results for the correlators K ϕϕ (τ ) and K ϕz (τ ) can be obtained from Eqs. (14) and (15) 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 [37], which only affects the zero-time self-correlators K ii (0).
Comparison with experimental results.-Experimental data have been taken in the same way as in Ref. [33] (see also [38]). Experimental parameters correspond to wellseparated frequency scales, as needed for the theoretical results, (T −1 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 [33], the effective ϕ includes a small correction δϕ = (κ ϕ − κ z )/2Ω R ≈ 0.036 (see [38]), 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 twolevel subspace after the run. The recorded signalsĨ i (t) are linearly related to the normalized signals I i (t) in Eqs.   (1) and (2) asĨ i (t) = (∆Ĩ i /2) I i (t)+Ĩ off i , where responses ∆Ĩ i have been calibrated using ensemble-averaged Ĩ i (t) (see details in [38]), giving in arbitrary units ∆Ĩ z = 4.0 and ∆Ĩ ϕ = 4.4. The offsetsĨ off i are approximately zeroed individually for each trace by measuring the non-rotating physical qubit after each run. Additional offset removal, |Ĩ off i | ≈ 0.15 − 0.20, for all traces with the same ϕ is done using Ĩ i (t) [38]. For calculating the correlators, we average over the ensemble of ∼200,000 traces and additionally average over time t 1 in Eq. (9) 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 opera-tor 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 2(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. (15). 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. (16). Note that because of the Markovian assumption, our theory is formally valid only for τ κ −1 i ∼ 30 ns; however, the experimental results agree with the theory even at τ < κ −1 i . Figure 2(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. 2(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 i , which vary from trace to trace. Note that the lines in 2(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. 2(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 . 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 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 fre-quencyΩ R in the experiment. From Eq. (15) we find Since in the case |Ω R | Γ z,ϕ we can neglectΩ R in Eq. (16) for Γ ± , Eq. (19) gives a direct way to findΩ R from the experimental antisymmetrized cross-correlator. The solid line in Fig. 3 shows K zϕ (τ ) − K ϕz (τ ) from the experimental data for ϕ = π/2. Fitting this dependence on τ with Eq. (19) (dashed line), we find the valuẽ Ω 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 Founda-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 [S1], 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 lumpedelement 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 Rabirotated frame, used in the experiment [S1] 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 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, where χ is the (small) dispersive coupling between the qubit and the measurement resonator mode, and ω m r 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 Ω −1 R . The sideband drive of the resonator at frequencies ω d ± Ω rf with equal amplitudes (here for simplicity we assume ω d = ω m r and Ω rf = Ω R ), produces the Hamiltonian term 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, This leads to the following dynamics of the resonator's coherent state |α(t) [or classical field α(t)] in the rotating frame based on ω d , where we also took into account the resonator damping with energy decay rate κ. Now let us solve the evolution equation (S6), assuming κ Ω R and |χ| Ω R . The drive term produces fast oscillation of α with Rabi frequency Ω R , Inserting this oscillation into the first term of Eq. (S6), using the trigonometric formula cos( , and neglecting oscillations with frequency 2Ω R , we obtain the equation for the slow evolution, Note that we can neglect the additional fast oscillations produced by the first term in Eq. (S6), ∆ 2 α(t) = −i(χ/Ω R )r 0 sin(Ω R t + φ 0 ) α s , in comparison with (S7), because we assume χ 2 /(Ω R κ) 1. We see that the evolution (S8) of the resonator field α s depends on the state of the effective qubit, which corresponds to the physical qubit (S1)-(S3) 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 where ρ(t) is the slowly-varying density matrix of the effective qubit. At this stage, we make use of the quantum Bayesian approach [S2-S4] to describe the qubit evolution due to measurement. Since the oscillating part ∆α of the resonator field [Eq. (S7)] 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 ∆α), 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 [S3-S6], the situation is equivalent to the qubit evolution due to measurement in the standard setup [S7] 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 [S3], which coincides with the Bayesian formalism for qubit measurement using a quantum point contact [S2]. 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 ∆Ĩ = ∆Ĩ max cos θ 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. (S12)]. The informational backaction is proportional to cos θ, while phase backaction is proportional to sin θ, with ensemble dephasing of the effective qubit, not depending on θ [S3, S4, S6]. Evolution of the effective qubit state due to measurement is given by Eqs. (18) and (25) of Ref. [S3] in the basis {|0 ϕ , |1 ϕ } (Eqs. (12) and (13) of Ref. [S4]). 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 = ω m r . If this is not the case, |ω d − ω m r | ∼ κ, then there will be an extra term −i(ω m r − ω d )α in Eq. (S6), which will lead to an extra term −i(ω m r − ω d )α s in Eq. (S8). Correspondingly, the steady states of the field α s for the effective qubit states |1 ϕ and |0 ϕ are (S14) instead of Eq. (S12), 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. [S4] and uses the field evolution equation (S8) instead of the steady-state solution (S12).

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 [S2, S4] in the Stratonovich form, where 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 ρ = (1 1+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 [S2, S4] When the observable σ ϕ is measured instead of σ z , these equations remain the same [S8] 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ô forṁ with When both σ z and σ ϕ measurements are performed at the same time, we simply add the terms from Eqs. (S19)-(S21) and (S22)-(S24) [S8] (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 doublesideband drive ε sin(Ω R t + ϕ) cos(ω d t) [see Eq. (S5)] directly determines the angle for the measured operator σ ϕ for the effective qubit. This followed from the approximate solution of Eq. (S6). As we will see below, a more accurate solution shows a small correction to the measured direction ϕ.
Neglecting the first term in Eq. (S6) but still keeping the last term, we obtain the (exact) oscillating solution which is more accurate than Eq. (S7). (Note that the additional term ∝ e −κt/2 is naturally included in the slow dynamics.) Inserting ∆α(t) into the first term of Eq. (S6), we obtain the following evolution of the slow part α s of the total field α = α s + ∆α: which in the case κ/Ω R 1 is approximatelẏ This equation coincides with Eq. (S8), except ϕ is replaced with ϕ + κ/2Ω R , thus slightly changing the direction of the measured operator for the effective qubit, 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 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 equatioṅ 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. (S31) 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Ω rf tσ y /2). This giveṡ Since Ω rf is much faster than the evolution of the effective qubit, in Eqs. (S33) and (S35) we can neglect the oscillating terms. Finally expressing T pd via T 2 and T 1 , so that T −1 which is Eq. (7) of the main text. Since Eqs. (S33)-(S35) 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 In our experiment these inequalities are well satisfied: T 1 = 60 µs, T 2 = 30 µs, |Ω −1 R | 10 µs, Γ −1 z = Γ −1 ϕ = 1.31 µs, κ −1 z = 37 ns, κ −1 ϕ = 22.1 ns, and Ω −1 R = 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. [S9]. It says that in order to calculate the ensemble-averaged correlator, 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. [S9]); 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 (S40) 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 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, 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 dt (S44) Diagonalizing the matrix M, we find the solution where Γ ± are the eigenvalues of the matrix −M, given by Eq. (16) of the main text and repeated here, Note that our evolution is symmetric under the inversion operation x → −x, y → −y, z → −z (it is a unital map), and therefore In this case Eq. (S41) simplifies to 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. (S45), which gives Inserting these results into Eq. (S48) 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 (S42) for the y-coordinate was not important for K zz (τ ) and K zϕ (τ ). Also note that we were able to use Eq. (S48) instead of Eq. (S41) 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 nonequivalent, and then we would need to use Eq. (S41). 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. [S9], we will do the derivation explicitly, essentially repeating the equivalence proof in [S9]. In the derivation we assume fixed t 1 and t 1 + τ [averaging the correlator (S40) 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. (S18) and (S25), 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)], 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 = [1 1 + 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 The second contribution (S57) to K zϕ (τ ) can be written as Combining the two contributions, we find K zϕ (τ ) = Tr(σ ϕρc ) z in + Tr(σ ϕ ∆ρ z ).
On the other hand, using Eq. (S41) of the collapse recipe with σ i = σ z and σ j = σ ϕ , we obtain which coincides with Eq. (S64). Thus, equivalence of both methods for K zϕ (τ ) is proven in the general (nonunital) 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. (S41) of the collapse recipe, which takes into account both scenarios (collapse to the state |0 or to |1 ) and not the simplified equation (S48) (collapse to |1 only), which is valid only for a symmetric (unital) evolution. As seen from Eq. (S41) [or Eq. (S64)], 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 [S2-S4, S6] 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 [S10], 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 [S10] 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 iHtotτ σ z e −iHtotτ and σ ϕ (τ ) ≡ e iHtotτ σ ϕ e −iHtotτ , 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 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 [S10] (related to what is usually called the Quantum Regression Theorem) 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. (S68) 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. (S68) the operator B is σ z , while A is either σ z or σ ϕ . The starting state for ρ av (τ |Bρ in ) is σ z ρ in , so Since we have to work with unphysical unnormalized states, we use ρ(t) = [P N 1 1 + 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 ) 1 1 2 using condition ρ 00,in + ρ 11,in = 1. Since the ensembleaveraged evolution is linear, we can separate ρ av (τ |σ z ρ in ) into three terms, corresponding to the three terms in Eq. (S71). The first term, σ z /2, gives the physical evolution ρ av (τ |1) starting with the state |1 . The second term, (ρ 11,in − ρ 00,in ) 1 1/2, does not change in time and gives zero contribution to the correlators (S69) and (S70). 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. (S69) and (S70) because the traces will be imaginary numbers.
Thus we obtain equations which coincide with Eq. (S48) 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 backaction), using representation (S61) of a linear quantum map, we can obtain Eq. (S64) from Eqs. (S70) and (S71), 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 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Ĩ off i 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Ĩ off i , 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. S1. 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 estimate the offsetsĨ off i , we use the symmetric combination S(t) ≡ [ Ĩ i (t) z0=1 + Ĩ i (t) z0=−1 ]/2, which is shown in Fig. S3 for both measurement channels and for 11 angles ϕ. We see that the S(t) are approximately independent of time, and therefore we can introduce the offsetsĨ off

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 selfcorrelator 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. S4) 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 [S10-S12], we need to consider the vacuum noisev(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 [S4]), 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. (S74), 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) iṡ 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. S4) with the same correlator (S74) and uncorrelated with v(t). In Eq. (S75) 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, whereκ = κ + 2i∆ω.
The outgoing field F (t), which is then amplified is 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 [S8] 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 [S8] 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 −