Quantum information science is poised to deliver unprecedented opportunities in terms of both fundamental physics and device technologies, by pushing existing boundaries in areas as diverse as quantum computation and simulation, secure communication, and quantumenhanced sensing and metrology. Notably, advances in quantum control and systems engineering are enabling access to intermediate-scale quantum processors whose capabilities are beyond what may be tractable classically [1], with impressive achievements having recently been reported [2][3][4][5]. Ultimately, however, realizing the full potential of these technologies will crucially depend on sustained progress in characterizing and overcoming the effects of noise that limit qubit performance.
In the context of entanglement-assisted quantum metrology, spatial and temporal correlations of the noise dictate the extent by which the standard quantum limit on precision may be overcome in the estimation of physical parameters [6,7], in particular when the quantummechanical nature of the environment must be explicitly accounted for [8]. Once characterized, spatial noise correlations may be exploited for augmenting the performance of quantum sensors via tailored quantum encoding [9] or error correction [10]. Likewise, quantitative knowledge about noise properties and their correlations may prove instrumental in designing optimized quantum-control or error-mitigation strategies that can be substantially more efficient than general-purpose schemes [11,12], as well as in determining optimal parameter regimes for "spectator qubits" to be useful in improving control performance of proximal data qubits [13].
Ultimately, noise correlations will play a key role in determining the feasibility of large-scale fault-tolerant quantum computation: establishing that noise correlations decay sufficiently rapidly in space and time is central for validating the locality assumptions under which the existence of an accuracy threshold may be rigorously derived beyond the paradigm of independent errors [14,15]. In turn, the structure of the noise, including the relevance of correlated error processes, is expected to strongly influence the value of the error threshold itself and inform ways in which resource-optimized architectures may be designed [16,17]. As a result, developing viable methodologies to detect and simultaneously characterize both spatial and temporal correlations present in realistic multiqubit noise environments is an imperative next step.
In a single-qubit setting, temporal correlations of dephasing noise that may be assumed to be stationary and Gaussian are characterized in the frequency domain by a single noise spectrum -namely, the Fourier transform of the two-point correlation function of the noise operator with respect to the time lag [18]. Estimation of the spectrum from experimental data may be achieved through various "quantum noise spectroscopy" protocols, which employ either pulsed or continuous control modulation of the qubit sensor to suitably shape its spec-
FIG. 1. Concept of two-qubit quantum noise spectroscopy experiment. Two qubits are coupled to a common bath, leading to spatiotemporally correlated noise. The reduced dynamics of the qubit system is modeled through a master equation which contains all the noise spectra. Using spin-locking control sequences, we measure decay curves of different observables. Comparison between the evolution of these observables predicted by the master equation against measured data yields the target spectra as fit parameters. While our methodology is device-independent, we use superconducting qubits for experimental validation.
tral response. To date, measurements of the noise spectrum have been reported across a wide variety of experimental qubit platforms, including nuclear magnetic resonance [19,20], superconducting quantum circuits [21][22][23][24][25], nitrogen-vacancy centers [26,27], spin donors in semiconductors [28], and trapped ions [29,30].
In a multi-qubit setting, complete characterization of Gaussian dephasing noise necessitates estimation of the full set of spectra {S jk (ω)}, defined by the Fourier transform of the correlation functions of noise operators acting on each possible combination of qubits j and k. While temporal noise correlations that affect qubits individually are now described in terms of self-spectra S jj (ω), coexisting spatial and temporal correlations are captured by the two-qubit cross-spectra {S jk (ω)}, with j = k. Spatial noise correlations have been probed and their strength upper-bounded in recent experiments using superconducting fluxonium qubits [31], nitrogen-vacancy centers in diamond [32] and spin qubits in semiconductors [33]. However, all the protocols implemented thus far lack the frequency sensitivity needed for full-fledged multi-qubit spectroscopy of noise that may be in general spatiotemporally correlated and non-classical. Despite promising theoretical proposals [34][35][36][37] as well as an experimental approach using a specific correlation measure [38], measurements of a two-qubit cross-spectrum remain yet to be reported.
In this paper, we theoretically develop a quantumcontrol protocol for two-qubit spectral estimation, and validate it experimentally in a circuit quantum electrodynamics (cQED) system using engineered photon shot noise. Building on continuous-control modula-tion spin-locking techniques previously implemented in a single-qubit setting [23,25], we formulate the problem in the framework of robust estimation theory [39] and demonstrate the simultaneous reconstruction of all the two-qubit self-spectra and cross-spectra that characterize our engineered noise source within a Gaussian approximation. In contrast to existing proposals employing dynamical-decoupling comb-based noise spectroscopy [34][35][36], our approach does not require the design and application of long sequences of nearlyinstantaneous pulses. As an additional advantage, for instance with respect to methods exploiting decoherencefree subspaces [33], no entangled states nor two-qubit gates are needed. Instead, our protocol relies on continuous driving of the individual qubits followed by simultaneous single-qubit readout. We thus anticipate this protocol to be compatible with any multiqubit architecture in which these resources are available.
Figure 1 illustrates the methodology underlying our two-qubit spectroscopy approach. To appreciate the underlying physical principles, it is useful to contrast our method with single-qubit noise spectroscopy via spinlocking relaxometry [23,25]: there, a microwave tone is applied on the sensor qubit to effectively create a "dressed qubit," whose level splitting equals the frequency of the Rabi oscillations induced by the drive. Since this dressed qubit is predominantly sensitive to the noise spectrum at its transition frequency, the (self-) spectrum may then be sampled, point by point, by measuring the decoherence rate the dressed qubit experiences as a function of the applied Rabi frequency. In our generalized two-qubit spin-locking protocol, we simultaneously drive two sensor qubits with amplitudes set to give an identical Rabi frequency, under ideal conditions. This effectively produces two dressed qubits, which are now sensing the same frequency component of the noise through the self-and cross-spectra that enter the master equation governing their dynamics. These spectra are then reconstructed at each Rabi frequency of interest by fitting the numerical solution of the master equation to experimentally measured decay curves.
The paper is organized as follows. In Sec. II, we provide the theoretical foundation of the reconstruction protocol we qualitatively described above. While Sec. II A presents the derivation of a master equation (ME) for the two-qubit reduced dynamics, Sec. II B is devoted to the robust statistical procedure we employ to extract the spectra from fits of time-dependent expectation values computed from the ME to experimental data. In order to set the stage for the experimental validation of our protocol, we devote Sec. III to the demonstration of engineered correlated photon shot noise in the cQED test-bed we employ (see Fig. 3), and to the experimental determination of the parameters entering the shot-noise spectra. We first perform a parametric analysis of this correlated noise in our experimental setup (Sec. III A), and then demonstrate our capability to selectively probe a specific frequency region of the cross-spectra (Sec. III B). After dis-cussing in Sec. IV A some technical modifications needed to adapt our protocol to the non-idealities specific to our circuit-QED platform, our main results are reported in Sec. IV B: namely, the non-parametric experimental reconstruction of the two-qubit self-and cross-spectra that characterize our engineered noise source. Full technical detail about the derivation of the MEs used to model the reduced dynamics in both the ideal and experimentally relevant parameter regimes is included in Appendix A, whereas Appendix B and Appendix C discuss various aspects pertinent to the robust statistical estimation procedure we employ.
A. Noise model and reduced master equation for spin-locking dynamics
We summarize here the basic ideas of spin-locking relaxometry [23,25] and generalize them to the two-qubit setting. To describe the driven evolution of two qubits under a noisy environment, we consider the Hamiltonian
where H S (t) is the time-dependent Hamiltonian of the two-qubit system S, H B the Hamiltonian of the bath B, and H SB describes system-bath coupling. We assume the two qubits to be characterized by angular-frequency splittings ω qj , and each to be coherently driven at frequency ω dj with a drive strength Ω j , j ∈ {1, 2}. Setting ≡ 1, we thus consider the system Hamiltonian
where σ x j and σ z j are the Pauli matrices for qubit j. Throughout this paper, the +1 and -1 eigenstates of σ z j will be denoted by |1 j and |0 j , respectively. While we leave H B unspecified, we specialize to single-axis systembath couplings of the form
where B j is the bath operator that couples to qubit j.
In the absence of coherent drives (Ω j = 0, j = 1, 2), H SB generates pure-dephasing evolution. However, in an appropriate rotating frame, the coherent drive "tilts" the quantization axis of each qubit by a π/2 angle, thus turning dephasing noise into a source of energy absorption and emission at rates that probe the two-qubit spectra at frequencies Ω j . To describe this phenomenon and exploit it for noise-spectroscopy applications, we apply the unitary tranformation R(t) = exp[-i j ω dj tσ z j /2] to move to a reference frame that rotates at the drive frequencies, and in which the Hamiltonian is H R (t) ≡ R † (t)H(t)R(t)-iR(t) † Ṙ(t). By effecting the frame transformation and invoking the rotating-wave approximation (RWA) to drop counter-rotating terms oscillating at frequency 2ω dj , the rotating-frame Hamiltonian is then
where H S ≡ j (∆ qj σ z j + Ω j σ x j )/2 and ∆ qj ≡ ω qj -ω dj is the detuning of drive j from ω qj .
The time-independent Hamiltonian H S can be diagonalized by transforming to a rotated spin-locking basis, spanned by |+x j ≡ cos(ϑ j /2)|1 j + sin(ϑ j /2)|0 j and |-x j ≡ -sin(ϑ j /2)|1 j + cos(ϑ j /2)|0 j , where ϑ j ≡ arctan(Ω j /∆ qj ) is the angle by which the qubit quantization axis is rotated under the drives. In this basis,
where we have introduced the Pauli matrices τ x j ≡ |+x -x| j +|-x +x| j and τ z j ≡ |+x +x| j -|-x -x| j . Setting ∆ qj = 0 ∀ j leads to ϑ j = π/2 ∀ j, which further simplifies the above Hamiltonians to
In this spin-locking basis, the coherent drives then define two effective "dressed" qubits quantized along τ z j , with angular-frequency splittings equal to the Rabi frequencies Ω j , and subject to purely transverse noise along τ x j . We may describe the evolution of the dressed qubits in the spin-locking basis by tracing out the bath and deriving a reduced ME for the density operator ρ(t) ≡ Tr B [ρ tot (t)]. As detailed in Appendix A 1, to do so we employ a standard time-convolutionless (TCL) ME approach [40]. We assume an initially separable initial state ρ tot (0) ≡ ρ(0) ⊗ ρ B , where ρ(0) and ρ B are the initial density operator of S and B, respectively. In addition, we consider stationary noise with zero mean, so that B j (t) = 0 and
, where B j (t) ≡ e iHBt B j e -iHBt is the time-dependent noise operator for qubit j in the interaction picture associated to the free bath Hamiltonian, the time lag τ ≡ ts, and • denotes expectation with respect to the initial bath state ρ B . We also assume that the coupling between the system and the bath is weak enough to truncate the TCL generator at second order, and employ a secular approximation to drop terms oscillating with frequency Ω 1 + Ω 2 . Setting Ω 1 = Ω 2 ≡ Ω, the two dressed qubits are most sensitive to the noise spectra in a frequency window of width ∼ 1/t around their splitting ±Ω. This enables us to simplify the reduced ME in the limit of a sufficiently long evolution time. For spectra S jk (ω) that vary sufficiently slowly with frequency about ω = ±Ω, we consider only the contribution of the spectra at ω = ±Ω, and finally arrive at
where the superoperators L jk are defined by
with {A, B} ≡ AB -BA denoting the anticommutator of A and B. The superoperators introduced in Eq. ( 9) describe correlated decay and absorption processes with strength proportional to the two-qubit spectra evaluated at ω = ±Ω. These two-qubit spectra are given by
The reduced ME defined by Eqs. ( 8)-( 9) involves all the spectra that are needed to characterize stationary noise that may in general be non-classical (in the sense that the commutator [B j (t), B k (s)] = 0, t = s), and display arbitrary temporal and spatial correlations in a twoqubit system -provided that noise can be assumed to be Gaussian and acting only along the σ z j axes (purely dephasing) in the laboratory frame. The Gaussian assumption may be satisfied exactly, for example for bosonic baths at equilibrium [41], or in an approximate sense when either a large number of independent bath degrees of freedom are involved [42] or the coupling is sufficiently weak for any higher-order cumulants to be negligible. It is worth noting that, from a rigorous spectral estimation standpoint, the assumptions of single-axis noise and Gaussianity cannot be expected to be valid a priori, and should always be verified experimentally through spectroscopic means [43][44][45].
Since Eqs. ( 8)-( 9) contain all the spectra of interest for the present two-qubit problem, an appropriate choice of experimental observations can enable us to infer S jk (Ω) for arbitrary j, k ∈ {1, 2}. For notational convenience, for a given Rabi frequency Ω, we collect the two-qubit spectra that we aim to estimate into a spectrum vector,
To devise a protocol for estimation of S, we will adopt an inverse-problems perspective and infer S by performing a non-linear regression that fits numerical solutions of the reduced ME to experimental data. Figure 2(a) illustrates the control and measurement cycle employed to gather experimental data from which 2. Protocol for spectroscopy of spatiotemporally correlated noise. (a) Control and measurement cycle. The qubit system S is initialized in state ρ(0) ∈ {ρs}, where s ∈ {1, 2, . . . , Nstates} (step 1). Continuous drives (shown as orange waves) with equal Rabi frequency Ω1 = Ω2 = Ω are then applied on the two qubits for time t ∈ {tq}, with q ∈ {1, 2, . . . , Ntimes}, during which S evolves under the influence of the noise spectra evaluated at ±Ω, {S jk (±Ω)} with j, k ∈ {1, 2} (step 2). After this evolution time, a projective measurement of a two-qubit observable O ∈ {Or}, with r ∈ {1, 2, . . . , N obs } is performed (step 3). For a given combination α of initial state ρs, evolution time tq, and observable Or, this cycle is repeated M times, and a sample mean Oα of all outcomes O (m) α is obtained [Eq. (11)], where m labels the outcome of the projective measurement for each cycle. The Rabi frequency Ω is swept to gather experimental observations for all frequencies at which S jk (±Ω) will be sampled. (b) Schematic of the reconstruction procedure of the twoqubit spectra from the data produced in (a). An initial value of Ω is chosen, for which an initial guess of the spectrum vector S(Ω) is assumed. This guess is fed into a non-linear regression algorithm to find the value of S(Ω) that globally minimizes the discrepancy between the measured sample means Oα (blue dots with error bars; note, in the first plot, the occurrence of an "outlier") and the corresponding expectation values Oα S obtained by numerically solving Eq. ( 8) for all chosen combinations α of initial state, evolution time, and observable. An estimate Ŝ(Ω n ) of S(Ω n ) at the n-th Rabi frequency Ω n is obtained following Eq. ( 13) (solid red lines). The latter is then used as the initial guess for the reconstruction at the next frequency Ω n+1 . The procedure is repeated until S(Ω) is reconstructed over all frequencies of interest.
S(Ω) is reconstructed at a given frequency Ω. At time t = 0, the system is prepared in ρ(0) = ρ s , where s ∈ {1, 2, . . . , N states } labels elements of an arbitrary set {ρ s } of two-qubit initial states.
Continuous drives resonant with each qubit are then applied with Rabi frequency Ω, so that the evolution of the two-qubit system is approximately given by the solution to Eq. (8). Multiple evolution times t ∈ {t q }, with q ∈ {1, 2, . . . , N times }, are considered, after which projective measurements of a system's observable O ∈ {O r } are performed, with r ∈ {1, 2, . . . , N obs }. Though we need not specify the initial states and observables at this stage, in the experiment presented here we will consider initial product states in the spin-locking basis, namely, ρ s = |ψ s ψ s |, with |ψ s ∈ {|+x, +x , |+x, -x , |-x, +x , |-x, -x }, and measurements of products of Pauli operators, for which O r = τ 1 1 ⊗ τ 2 2 , with 1 , 2 ∈ {0, x, y, z}, and where τ 0 j ≡ I j is the identity operator for qubit j. These initial states and observables are accessible through simultaneous preparation and measurement of each qubit, and thus using purely local resources.
Figure 2(b) illustrates the procedure by which experimental observations resulting from the control and measurement cycle described above are used to reconstruct the spectrum vector S(Ω). To simplify the notation, we label all combinations of initial states ρ s , evolution times t q , and observables O r using a single collective index α ∈ {1, 2, . . . , d}, where d = N states × N times × N obs . In addition, for each α, we consider sample means O α of all outcomes O (m) α of projective measurements m ∈ {1, 2, . . . , M }. These sample means are defined by
where M is the total number of projective measurements, which we take to be the same for each α for simplicity. It will also be convenient to collect all such sample means measured experimentally for a given Rabi frequency into a single observation vector, O ≡ (O 1 , O 2 , . . . , O d ) T . In the asymptotic limit of a large number of projective measurements, M → ∞, the sample means O α converge to their expectation values O α S by the weak law of large numbers, O α → O α S . These expectation values are determined by S through
where ρ(t qα )| S,sα is the solution of Eq. ( 8), for noise spectra S and initial state ρ sα at time t qα . Hence, for finite M , we estimate the spectra by finding the value of S that minimizes the deviation between the data O α and the predictions of the model O α S . Formally, we define our estimator of S for a particular Rabi frequency Ω as
Here, λ(z) is a loss function that penalizes deviations between the model and the data, which are quantified by the normalized residuals z S,α ≡ (O α -O α S )/σ α , where σ 2 α ≡ var(O α ). Throughout this text, we use hats to denote estimators. Estimators like Eq. ( 13), which minimize a total cost function, are called M-estimators [39].
The most natural choice of loss function is arguably the quadratic function λ(z) = z 2 /2, since Eq. ( 13) then reduces to a simple weighted least-squares estimation. In addition to being well-suited for numerical optimization, the weighted least-squares estimate is statistically well motivated when the probability distribution of O is Gaussian. Indeed, in this case, weighted least-squares optimization can be derived from maximum-likelihood estimation of S, and can be shown to be asymptotically efficient, that is, the variance of Ŝ achieves the Cramér-Rao bound on precision [46].
In practice, however, the statistics of sample means O may not be perfectly described by the expected Gaussian probability distribution, compromising the asymptotic efficiency of Ŝ under weighted least-squares estimation. In particular, experimental data is often contaminated by outliers: data points that do not follow the probability distribution of the majority, for example because of isolated experimental errors. Because of its quadratic dependence on the residuals, weighted least-squares estimation notoriously gives excessive weight to outliers that are distant from normal observations, namely, for which z S,α
1. This can make weighted least-squares estimators inefficient, causing estimates to wander very far away from their expected behavior (as will be seen experimentally in Sec. IV and Appendix C), or, in the worst case, leading to catastrophic divergent behavior of the variance of the estimator. Overcoming these limitations motivates the use of robust estimators, whose performance is not significantly impaired by the presence of outliers [39,46]. A prevalent way to achieve robust estimation is to employ the Huber loss function, a mixture of weighted least-squares estimation with mean absolute error minimization defined by
where δ 0 ∈ [0, ∞[ is a tuning parameter that controls the mixing. For |z| ≤ δ 0 , the Huber loss function provides the statistical and numerical efficiency of weighted least squares, while residuals with |z| > δ 0 only contribute to the total cost through their absolute value. This avoids overweighting outliers in the observations and favors robustness of the estimator, as desired.
In the experiment we will present, we employ the superconducting circuit shown in Fig. 3 two-qubit spectroscopy based on the spin-locking technique described in Sec. II. The two qubits are encoded by the lowest energy levels of a pair of transmons dispersively coupled to a common bath formed by a mode of a common resonator, which is brought to a steady-state population under the competing action of a constant applied microwave drive and photon emission into an external environment. When the coupling between the qubits and this bath is sufficiently weak, and the evolution time is sufficiently long, the common bath may be traced out following the derivation presented in Appendix A 2, and the spin-locking experiment may be described by Eq. ( 8) for evolution of the two-qubit system under the photon shot-noise spectra [18,25]:
with j, k ∈ {1, 2}. Here, χ j is the strength of the dispersive coupling for qubit j, κ quantifies the resonator damping rate, ∆ c ≡ ω c -ω d is the detuning between the drive frequency ω d and the bare frequency of the common resonator ω c , and n is the average photon number in the steady state of the resonator. While the parameters κ, χ 1 , and χ 2 are set by design during device fabrication, both the amplitude and asymmetry of the spectra around zero frequency may be tuned in situ via n and ∆ c by varying the strength and detuning of the drive applied onto the common resonator, granting us the capability to produce spatiotemporally correlated noise with an engineered spectrum.
To validate this capability, in this section we demonstrate the presence of noise correlations consistent with photon shot noise. We first use a Ramsey interferometry technique to witness spatial correlations by measuring a two-qubit correlation function in a free-evolution setting. Fitting the results to the solution of the quantum-optical ME describing the joint evolution of the qubits and the resonator enables us to measure the parameters χ 1 , χ 2 , and κ entering Eq. ( 15). We then use the spin-locking technique to measure noise correlations in a frequencysensitive fashion and demonstrate experimental control over the amplitude and asymmetry of the engineered noise spectra through both n and ∆ c .
As mentioned, we create a source of correlated photon shot noise by applying a coherent signal to the common resonator. We then perform simultaneous Ramsey sequences by applying a pair of π 2 -pulses on each qubit [Fig. 4(a)]. The qubit drives are detuned from the qubit frequencies by ∆ qj = ω qj -ω dj , with ∆ q1 /2π = -1.26 MHz and ∆ q2 /2π = 0.3 MHz. To observe Ramsey fringes, we vary the wait time t between the pulses applied on each qubit, as well as the injected photon number n in the common resonator, which was calibrated in advance by varying the power of the noise drive and observing the frequency shift of the qubit due to the increasing resonator population [47]. After each Ramsey sequence, we perform simultaneous single-shot dispersive readout of both qubits in the σ z j -eigenbasis through their individual resonators, and reinitialize the qubits in the ground state by letting them relax for 500 µs, a time much longer than the observed qubit relaxation time T 1 . During the readout we turn off the otherwise continuous noise drive to improve the single-shot readout fidelity. Taking the sample means σ z 1 and σ z 2 of the resulting readout outcomes for each qubit [Eq. ( 11)], we obtain the Ramsey fringes shown by the color schemes in Fig. 4(b) and Fig. 4(c), respectively. To access the correlation induced by the shared noise source, we then calculate the sample mean of measurement outcomes of the two-qubit observable (a) Ramsey sequence applied simultaneously to the two qubits, with π 2 -pulses around the y-axis (blue boxes). A coherent drive (the "noise" tone, shown in green) is turned on well before the first pair of π 2 -pulses, to bring the common resonator into a steady state with finite population, thus implementing a stationary source of engineered correlated noise. We read out both qubits simultaneously through their individual resonators, using drives at frequencies ωr1 and ωr2 (red box), thus acquiring a series of ±1 outcomes for each qubit. (b)-(c) Sample mean of measured single-qubit observables σ z 1,2 for different wait times t and different injected photon numbers n in the common cavity. The frequency of the observed Ramsey fringes varies with increasing n due to the dispersive shift 2χjn of the qubit frequencies (tilt of the vertical red and blue lines). For higher n, we see a rapid dephasing of the qubit states due to the added photon shot noise (blurring of the features in the top right corners).
(e) Solution of the ME describing the qubits coupled to the common resonator [Eq. ( 16)], in which parameters κ, χj, ∆qj, and γ φj (see Table I the quantum-optical ME [48]:
Here, ρ QC (t) denotes the joint density matrix for the two qubits plus the resonator mode, which evolves under a dispersive Hamiltonian of the form [49,50]
with ε being the amplitude of the drive applied on the resonator. In Eq. ( 16), the Lindblad superoperators L R
x and L R z are defined by (17) and describe the effect of noise coupling to qubit j through σ x j and σ z j , leading to relaxation and dephasing at rates Γ 1,j = 1/T (j) 1 and γ φj , respectively. In Eq. ( 16), L R z phenomenologically describes any uncorrelated source of noise that may couple to σ z j in addition to photon shot noise.
Solving Eq. ( 16) using standard numerical packages [51], we calculate the time-dependent expectation The qubit and resonator frequencies are measured spectroscopically, and the reported relaxation and coherence times are averages of repeated T1decay and echo measurements over 24h, with T2 < 2T1 indicating native sources of dephasing noise. The coupling constants κ and χj, the qubit detunings ∆qj and the native dephasing rates γ φ j , are results of the fits to the Ramsey measurements (Fig. 4). These values are found to agree within error bars with those obtained from independent measurements.
after two instantaneous Ramsey pulses, and estimate relevant physical parameters by fitting to the measurements of C zz displayed in Fig. 4(d). Through this procedure, we obtain the values of κ, χ 1 , χ 2 , ∆ q1 , ∆ q2 , γ φ1 , and γ φ2 collected in Table I. Throughout the simulations, we take ∆ c = 0, and use Γ 1,j measured independently from relaxation experiments in the absence of injected photons. The results of the simulations using the fitted parameters are displayed in Fig. 4(e). The remarkable quantitative agreement between simulation and experiment obtained here [see also the linecuts in Fig. 4(f) for representative n values] provides strong evidence that the measured correlations arise from our engineered source of photon shot noise.
Several important features of this Ramsey interferometry experiment may be understood qualitatively. First, C zz vanishes when the noise drive is off (n = 0). That is, having (approximately) no photons entering or leaving the common microwave cavity turns off the interaction responsible for a non-zero correlation between the two qubits, as desired. Second, C zz becomes "blurry" in the top right corner of Fig. 4(d) for high n and t. As we increase the amplitude of the photon shot noise, we induce correlated dephasing, implying that more photons enter and leave the common resonator, each one causing the qubits to simultaneously dephase as it takes with it information about the phase of the qubits. At a high enough photon number, the dephasing rate is the dominant source of decoherence, causing the signal to decrease significantly within the wait time of 10 µs between the Ramsey pulses. Furthermore, we observe oscillations of C zz at frequencies |∆ q1 -∆ q2 |/2π = 1.56 MHz and |∆ q1 + ∆ q2 |/2π = 0.96 MHz. These two overlapping signals cause a beating pattern at their difference frequency 0.6 MHz, which corresponds to the blurry white lines that appear to enter Figs. 4(d)-(e) from the bottom left (indicated by dashed arrows).
In the previous section, we observed the signature of a spatially correlated noise source on two qubits. However, spectral estimation requires the ability to select the frequency at which the noise spectrum is probed. To this end, we use the generalized spin-locking technique we described in Sec. II. Specifically, our protocol is illustrated in Figs. 5(a)-(b). Two sets of simultaneous π 2 -pulses are applied along the y-axis and separated by a time t, during which the qubits are driven by constant, resonant drives. The first pair of pulses initializes both qubits in state |-x . As explained in Sec. II A, the drives define dressed qubits whose eigenstates |±x j are split by the Rabi frequencies Ω j . During the time t, these dressed qubits undergo a non-unitary evolution due (predominantly) to noise coupling to σ z j at angular frequencies Ω j . The second pair of π 2 -pulses returns the qubits back to the z-axis, where they are immediately measured by dispersive readout. We use this sequence to probe correlations in qubit dephasing by sweeping the Rabi frequency Ω 1 across the fixed Rabi frequency Ω 2 . The peak frequency of our engineered noise being at -∆ c [Eq. ( 15)], we set Ω 2 = |∆ c | to maximize the contribution of shot noise in qubit dynamics. We expect to see a non-zero correlation,
, only if both dressed qubits are sensing the same frequency component of a non-vanishing spectrum (Ω 1 ≈ Ω 2 ). This insight is confirmed by the experimental data shown in Fig. 5(c), in which K zz is only significant around Ω 1 = Ω 2 = |∆ c |, in a frequency region that narrows with the duration t over which the spin-locking drive is applied. This follows from a key feature of the qubit evolution during the spinlocking drive: namely, the terms in the ME that contain the cross-spectra oscillate at ±|Ω 1 -Ω 2 | [see Eq. (A11)] and average out for times 1/|Ω 1 -Ω 2 |. Thus, setting Ω 1 = Ω 2 lets us isolate the influence of the cross-spectra in a specific frequency region.
To verify compatibility of the observed correlation with photon shot noise, we derive the quantum-optical ME for the two qubits and the resonator in the spin-locking frame. This leads to an expression formally identical to Eq. ( 16), but in which we replace
Above, L SL z phenomenologically describes uncorrelated sources of noise coupling to σ z j in addition to photon shot noise, in the spin-locking frame. In the ME simulations,
5. Spin-locking experiments. (a) Spin-locking sequence applied simultaneously on two qubits with a shared source of engineered noise. (b) Bloch sphere sketch of the sequence in (a) for qubit 1 and qubit 2. The initial π 2 -pulse rotates the qubits from the ground state to |-x . The spin-locking drive effectively creates a dressed two-qubit system, with each level splitting being equal to Ωj. Due to dephasing noise (from both the injected photon shot noise and from weaker native sources), the dressed qubits decay to their steady-state values along the x-axis. The final π 2 -pulse turns the qubits back to the initial (σ z j ) quantization axis to be measured.
, where the Pauli matrix τ z j is diagonal in the spin-locking basis {|+x j , |-x j }. The coherent drive creating photon shot noise is detuned by ∆c/2π = -2.03 MHz away from the common resonator. The Rabi frequency Ω2 of the spin-locking drive of qubit 2 is held constant at Ω2 = |∆c|, while Ω1 is swept (y axis). The correlation is significant in the region indicated by the solid black line, with width ∼ 2/t. (d) Numerical solution of the quantum-optical ME [see Eqs. (18) to (20)]. Inset: Qualitative shape of the cross-spectrum S12 [Eq. ( 15 we take values of χ j and Γ 1,j = 1/T (j) 1
given in Table I. In addition, we take the bare cavity drive detuning to be ∆ c = ∆ (00) c + χ 1 + χ 2 , where ∆ (00) c /2π = -1.95 MHz is the cavity drive detuning obtained from experimental measurements of the Stark-shifted resonator transmission peak when both qubits are in their ground state. This yields ∆ c /2π ≈ -2.03 MHz; the remaining parameters are then obtained by fitting the solution of the ME to the experimental data shown in Fig. 5(c), leading to
2 ≈ 9 × 10 3 rad/s, and γ ↓ 2 ≈ 14 × 10 3 rad/s. The resulting correlation, shown in Fig. 5(d), displays strong quantitative agreement with experimental observations. Though the phenomenological decay rates γ ↑ j and γ ↓ j are not entirely negligible, they lead to decay on a timescale 100 µs, while the dynamics due to photon shot noise occur on a shorter timescale, 50 µs. This confirms that our engineered noise source predominantly drives the dynamics in this experiment, as intended.
We finally discuss an additional intriguing feature in the results of the experiment. While the correlation starts out at zero [see Fig. 5(e)] when we initiate the spin-locking drive, and then quickly rises to a maximum of about 0.13, K zz ultimately decays over several tens of µs to a negative value of about -0.03. Numerical simulations predict that, despite the presence of intrinsic qubit decay sources (non-zero Γ 1,j , γ ↑ j and γ ↓ j ), the system subsequently reaches a steady state with K zz ≈ -0.025. While this negative correlation value may seem puzzling at first, it may be readily understood from the steady state obtained numerically, which contains significant coherences between |+x, -x and |-x, +x : in turn, this is due to the exchange of dressed qubit excitations mediated by the common source of photon noise.
IV. VALIDATION OF TWO-QUBIT QUANTUM NOISE SPECTROSCOPY
In Sec. III, we used fits to experimental data to estimate the parameters entering the two-qubit spectra of photon shot noise from a coherently-driven common resonator mode, Eq. (15). We now use this information to achieve our key objective: validating the spectroscopy protocol presented in Sec. II, by adapting it to the nonidealities specific to our circuit-QED platform and comparing the resulting experimental spectrum reconstructions with ideal spectra.
To produce engineered noise with known spectra, we use a microwave tone to give n = 0.127. Setting the Stark-shifted resonator-drive detuning with both qubits in their ground state to ∆ (00) c /2π = 2.05 MHz results in a Lorentzian spectrum peaked at angular frequency
To produce the experimental data needed to reconstruct these spectra, we apply the two-qubit spin-locking sequence illustrated in Fig. 2(a) and Fig. 5(a)-(b), by letting Ω 1 ≈ Ω 2 ≈ Ω ≡ (Ω 1 + Ω 2 )/2 to have both qubits sample the spectrum at ω = Ω, thereby maximizing the sensitivity to the noise spatial correlations. We perform a total of 26 spin-locking experiments, between which the Rabi frequency Ω/2π is swept through 26 values uniformly distributed to probe the Lorentzian peak from -2.2 to -1.8 MHz, along with the corresponding positive frequencies from 1.8 to 2.2 MHz. In the spin-locking experiments, we use the 4 initial states, 11 observables, and 26 evolution times given in Table II. To obtain the sample means of the observables, we average over M = 10 4 simultaneous projective measurements of all 9 combinations of Pauli matrices τ 1 1 and τ 2 2 , 1 , 2 ∈ {x, y, z}, thus performing two-qubit state tomography for each data point [52]. Separate numerical simulations indicate that the spectra may be accurately reconstructed with fewer initial states and observables for sufficiently highquality data; nevertheless, we find that full tomography with the 4 initial states of Table II is more useful in practice, enabling us to diagnose and fix experimental issues such as imperfect calibration of the measurements.
After completing all the spin-locking experiments, we condense the tomography data into sample means of projective measurements of single-qubit observables τ z j , j ∈ {1, 2}, and two-qubit observables of the form
, 1 , 2 ∈ {x, y, z}, along with corresponding standard deviations. To accurately reconstruct the engineered noise spectra, we find that the ME derived in Sec. II A under pure-dephasing noise (in the lab frame) and used for non-linear regression in the procedure described in Sec. II B must be adapted to two types of nonidealities in our cQED setting:
, where {τ j }, ∈ {x, y, z}, is the set of Pauli matrices for the dressed qubit j and the bar indicates a sample mean over M = 10 4 projective measurements.
(i) Finite Rabi-frequency difference.-In Sec. II A, we assumed Ω 1 = Ω 2 to arrive at the reduced ME, Eq. ( 8). Experimentally, however, we observe that the amplitude of the spin-locking drives can drift over a timescale of a few hours, most plausibly due to drifts in electronics, thus making δΩ ≡ Ω 1 -Ω 2 = 0. As a consequence, in the interaction picture with respect to H S [explicitly given below Eq. ( 4)], terms involving cross-spectra in Eq. ( 8) oscillate at frequency δΩ, which significantly suppresses their influence over times 1/δΩ (see Appendix A 3). This biases any estimate of the spectra based on Eq. ( 8).
(ii) Relaxation noise.-Although the engineered dephasing mechanism used in this experiment is made dominant by applying a sufficiently strong resonator drive so that n is large, superconducting qubits always suffer from significant intrinsic noise coupling to σ x j . This leads to T 1 relaxation (e.g., from Purcell decay [53]) in the lab frame. Such noise has distinct dynamical effects that are not captured by Eq. ( 8), and thus biases any estimates of S jk (ω) based on Eq. ( 8) with respect to their true physical value.
In order to simultaneously account for both types of non-idealities, we derive a modified ME for the reduced two-qubit dynamics. Under the approximations described in Appendix A 3, we arrive at
where Ω is now defined as the average Rabi frequency, Ω ≡ (Ω 1 + Ω 2 )/2. In addition, in Eq. ( 21), {L jk } is the usual set of Lindblad superoperators for correlated dephasing noise given by Eq. ( 9), and Γ 1,j ≡ 1/T (j)
1 , with T (j) 1 the longitudinal relaxation time of qubit j in the lab frame (without spin-locking drives). We then apply the spectral reconstruction procedure described in Sec. II B, with the exception that Eq. ( 21) is used in lieu of Eq. ( 8) in the calculation of expectation values involved in the non-linear regression procedure, Eq. ( 12).
In comparison with Eq. ( 8), the modified ME in Eq. ( 21) involves three additional parameters:
1 , and δΩ. While the longitudinal relaxation times are known from prior characterization of the qubits in the lab frame (Table I), the Rabi frequency differences are due to slow random drifts, and cannot be known in advance. However, δΩ can be accounted for as an additional unknown parameter in the estimation scheme; formally, we simply replace S → θ ≡ [S, δΩ] in Eq. ( 13) and simultaneously estimate S and δΩ. In this approach, we thus assume that δΩ is a constant parameter for a given reconstruction at target Rabi frequency Ω, but allow δΩ to vary between spin-locking experiments aiming to reconstruct spectra about distinct target Rabi frequencies, thus modelling slow, quasi-static Rabi-frequency drifts.
Figure 6 summarizes the results of the spectral reconstruction procedure described above, using a spectrum vector that is uniform across its components, S = 1 kHz ∀ , and δΩ = 0, respectively, as initial guesses for the non-linear regression technique we employ. In Fig. 6(a)-(c), we show examples of experimentally measured decay curves with the dressed qubits initialized in state |+x, +x for Ω/2π = 2.184 MHz (blue dots) and Ω/2π = 1.976 MHz (yellow dots) for three of the eleven observables given in Table II. The fitted decay curves are shown by solid lines of corresponding color, and display reasonable agreement with data. The spectra that follow from this robust estimation procedure are shown in Fig. 6(d), in which the gray error bars indicate the 95% confidence intervals derived from the asymptotic statistics of M-estimators (Appendix B). These reconstructions are compared with shaded orange areas representing 95% confidence intervals for the theoretical shot-noise spectra [Eq. ( 15)], using estimates for χ 1 , χ 2 , n, and κ obtained as explained in Sec. III. The four experimental reconstructions clearly capture the asymmetry of the two-qubit spectra arising from the non-commuting nature of the noise operator, and qualitatively reproduce the Lorentzian shape associated with photon shot noise. Remarkably, this includes a reconstruction of the crossspectrum S 12 (ω) that characterizes spatial correlations of the noise.
To produce the reconstructions shown in Fig. 6, we considered all initial states, evolution times and observables given in Table II in a global fit defined by the Mestimator given by Eq. ( 13) with S replaced by θ. To minimize the total cost function in Eq. ( 13) numerically, we employed the least-squares optimization routine of the SciPy package, which implements a trust-region reflective algorithm that allows the use of arbitrary loss functions [54]. While the quadratic loss function is the most standard in nonlinear regression, we find that due to the significant presence of outliers in the experimental data, quadratic loss leads to fitted decay curves that can deviate from the large majority of experimental observations. To mitigate this adverse behavior, we thus perform robust estimation using the Huber loss function, Eq. ( 14), in which we set the tuning parameter δ 0 = 1, leading to less noisy reconstructions and significantly better agreement between fitted decay curves and experimental data (see Appendix C for further discussion). This result showcases the advantage of robust estimation strategies over traditional weighted-least squares estimates in a quantum noise spectroscopy context. Indeed, for each of the 26 Rabi frequencies considered here, the current approach involves, as mentioned, 4 initial states, 11 observables, and 26 evolution times, for a total of 29, 744 data points. For such a large dataset, systematically identifying, explaining, and eliminating outliers in an experimental setting would represent, at best, an extremely tedious and impractical task.
Though the reconstructed and predicted spectra agree within error bars for several frequency values, statistically significant deviations are also observed. In particular, the reconstructed spectra are larger than predicted by an amount 10 kHz for wide ranges of frequencies. To investigate the physical origin of this discrepancy, we first simulate the spectroscopy procedure by exactly solving the coupled evolution of the two qubits and the driven resonator mode using the full quantum-optical ME defined by Eqs. ( 18)-( 20), setting γ ↑ j = γ ↓ j = 0 ∀ j and Ω 1 = Ω 2 = Ω. We then employ the resulting density matrix to calculate the probabilities of all relevant measurement outcomes, and use these probabilities to produce 10, 000 simulated projective measurements for each data point measured experimentally. The resulting decay curves are then fitted using the approach presented in Sec. II, producing the spectrum estimates illustrated by the blue dots in Fig. 6. A small discrepancy with the ideal spectra for the same shot-noise parameters then arises because the assumptions used to arrive at the reduced ME for two-qubit evolution are only valid in an approximate sense. In particular, noise may not be sufficiently weak, and the filters for finite evolution time may not be sufficiently narrow for Eq. ( 21) to hold exactly. Nevertheless, this discrepancy remains too small to explain the excess noise measured experimentally.
To attempt to explain the observed discrepancies, we invoke non-idealities that, taken together, may help to better understand our experimental results. First, native dephasing noise sources (in addition to engineered shot noise) may couple to the qubits via σ z j . We expect the spectrum of this excess noise to add up with engineered noise in the spectra S jk (ω) appearing in Eq. ( 9); such noise would thus also be measured by our protocol. The presence of intrinsic noise would be consistent with the numerical fits performed in Sec. III B, which yielded decay rates between 2 × 10 3 rad/s and 14 × 10 3 rad/s, phenomenologically accounting for spatially uncorrelated noise along σ z j . In addition, characterization of the qubits without engineered noise led to T 2 values significantly shorter than 2T 1 (see Table I). Native dephasing noise may arise from a combination of residual thermal photons in the readout resonators [25] with sources of noise intrinsic to transmon qubits, such as two-level fluctuators. In 15), with parameters taken from fits described in Sec. III: κ/2π = 198 kHz, χ1/2π = -29.1 kHz, χ2/2π = -59.5 kHz, ∆c/2π = 1.961 MHz, n = 0.127. Blue dots: numerical simulation of the spectroscopy protocol using the quantum-optical ME in Eqs. ( 18)- (20). Inset of the bottom-right panel: Rabi-frequency difference δΩ.
particular, the E J /E C ratios in our sample, which determine the sensitivity of the transmon frequency to charge noise, were 28 for qubit 1 and 45 for qubit 2. These are well below the value of 50 normally associated with transmon devices [55,56]. We also observe that Rabi frequencies significantly drift during the 24 hours needed to complete all the spinlocking experiments. Indeed, the inset of the bottomright panel of Fig. 6 shows that a Rabi frequency difference δΩ/2π ∼ 50 kHz builds up as the target value of Ω/2π is swept from 1.8 MHz to 2.2 MHz, most plausibly due to drifts in electronics. It is thus very likely that the average Rabi frequency Ω/2π ≡ (Ω 1 + Ω 2 )/2 defining the frequency at which the spectrum is reconstructed drifts by a similar quantity. In addition, both T 1 and T 2 are known to fluctuate significantly under noise processes that naturally occur in superconducting qubits [57,58].
These combined effects may thus explain the apparent shift of some reconstructed spectrum values away from the predicted Lorentzians.
Finally, the two reconstructed values of S 22 (ω) at frequencies nearest to ω/2π = -1.8 MHz deviate significantly from predictions. For these Rabi frequencies, we find that the nonlinear regression fails to produce a good fit to the experimental decay curves. We attribute this effect to the weakness of the engineered noise source at these frequencies, leading to a poor signal-to-noise ratio in the measured dynamics that prevents convergence of the nonlinear regression to the correct spectrum values. By artificially down-weighting the residuals associated with two-qubit observables in the loss function, we were able to achieve a more accurate reconstruction of the self-spectra for these frequencies; however, this came at the expense of an inaccurate reconstruction of the cross-spectrum. To achieve simultaneously accurate reconstructions of all spectra far from the peak frequency, it is possible in principle to perform more projective measurements to increase signal-to-noise and thus ease the convergence of the non-linear regression procedure into a physically meaningful global minimum.
We have proposed a protocol for two-qubit noise spectroscopy and validated it with an engineered source of photon shot noise in a superconducting-qubit architecture. Despite the complexity of the two-qubit dynamics and the resulting estimation problem in the presence of both non-classical self-and cross-correlation spectra, we were able to successfully extend the spin-locking technique previously used on a single qubit to the two-qubit setting. This enabled us to demonstrate what constitutes, to the best of our knowledge, the first experimental reconstruction of a two-qubit noise cross-spectrum.
Our approach offers several advantages over available proposals for two-qubit spectroscopy. Indeed, our continuous-wave protocol avoids the experimental issues arising from long tailored sequences of nearly instantaneous pulses that are typically required in comb-based dynamical-decoupling spectroscopy [34]; it does not require two-qubit gates or entangled initial qubit states; and it allows one to probe noise in a MHz frequency scale that is difficult to access with pulsed techniques. In addition, the numerical approach taken here, in which the numerical solution of a ME containing the spectra of interest is fitted to experimental data, offers substantial flexibility in the choice of theoretical model. This enabled us to quantitatively account for non-idealities of our experimental platform in a straightforward manner, by simply modifying the ME to include T 1 effects and Rabi-frequency drifts. We expect the simplicity and flexibility of our approach to be instrumental in its application to study correlated noise sources in other state-ofthe-art quantum devices -including nuclear-spin, phonon or charge noise affecting spin qubits in semiconductors or NV centers [33,59,60].
Our work also paves the way to further advances of quantum noise spectroscopy protocols and to the characterization of more general sources of noise. Indeed, applying the same non-linear regression techniques to a different master equation should enable spectral estimation to be extended to multi-axis noise, in a simpler continuouswave approach than in existing pulsed schemes [45] -with further extensions possible to a multi-qubit scenario. In addition, as photon shot noise is genuinely non-Gaussian, increasing the strength of the engineered noise studied here should enable one to use non-Gaussian noise spectroscopy methods [43,44] to investigate a quantum non-Gaussian environment with non-commuting degrees of freedom. As noise spectroscopy methods will continue to grow in complexity and generality, more elaborate mea-surements and larger experimental data sets will also be more likely to contain outliers. We thus envision that robust estimation strategies based on M-estimators, as we demonstrated here, will become increasingly needed and prove useful for quantum sensing applications.
The TCL formalism enables the systematic derivation of a ME for a system of interest S coupled to a bath B by tracing out the bath degrees of freedom. Though the method presented here is more general, to apply the TCL formalism to the cQED sytem studied experimentally in Sec. III and Sec. IV, we will assume that the bath B may be separated into two components: a central bath labeled by C, corresponding to a mode of the microwave resonator shared by the two qubits in the main text, and a larger external bath labeled by E, corresponding to the environment of the resonator mode in the main text. While the central bath couples to the system, the external bath only couples to the central bath; its effect is to allow non-unitary evolution of the central bath even in the absence of coupling to the system. In our cQED context, this allows us to describe damping of the microwave resonator. Beyond the current setting, such "structured environments" can naturally arise, for instance, for qubits coupled to two-level charge fluctuators undergoing incoherent transitions due to an electronic or phononic reservoir [42,61].
More formally, let the system-bath Hamiltonian be
where H SB ≡ H SC ≡ j B j Q j describes coupling between qubit j and C through the qubit operator Q j and the central-bath operator B j . Moving to the interaction picture with respect to H S + H B , we have
where the transformed operators Bj (t) ≡ e iHBt B j e -iHBt and Qj (t) ≡ e iHSt Q j e -iHSt . To derive a TCL ME describing the reduced evolution of the qubit system only, we introduce the projection superoperator P that projects any density matrix ρ onto the relevant (system) part of the Hilbert space: that is,
where the trace is taken over both central and external baths and ρ B is the joint initial state of C and E.
A complementary projection superoperator Q on the irrelevant part of the density matrix is then also defined by Q ≡ I -P, where I is the identity superoperator, Iρ ≡ ρ. As customary, we assume that the density matrix ρ tot (0) describing the joint initial state of the system and the bath is of the form ρ tot (0) = ρ(0) ⊗ ρ B , where ρ(0) is the initial density matrix of S. Further assuming that coupling between the system and the central bath is sufficiently weak, the TCL ME may be truncated at second order, leading to
Here, ρtot (t) denotes the joint interaction-picture density matrix of the system, central and external baths at time t, and K(t) is the second-order TCL generator given by
with L(t)ρ ≡ -i[ HSB (t), ρ] defining the Liouvillian superoperator associated with HSB (t). Importantly, Eqs. (A3) and (A4) are valid when PL(t)P = 0 ∀t, a property that is satisfied when Tr B [ Bj (t)ρ B ] = 0, i.e., for noise with vanishing mean in state ρ B . Under the above assumption, substituting Eq. (A2) into Eq. (A4) and tracing over both central and external baths produces the integro-differential equation
which describes evolution of the reduced density matrix of S, ρ(t) ≡ Tr B ρtot (t), in terms of the two-point correlation functions
(A7) Equality between Eq. (A6) and Eq. (A7) is only respected when noise arising from the central and external baths is stationary, so that correlation functions are invariant under time translations. Noting that Bj (t) = e iHBt B j e -iHBt is formally equivalent to the Heisenbergpicture evolution of B j under the Hamiltonian H B of the bath only, C jk (t, s) is independent of the evolution of the system. Stationarity then arises in the following two situations, the second one being the most relevant to this paper: (i) The initial density matrix of the bath commutes with H B , [H B , ρ B ] = 0. (ii) Evolution of the reduced density matrix of the central bath in the absence of the system, namely,
is accurately described by a ME of the form ρC (t) = L C ρ C (t), where L C is a Markovian (Lindblad) superoperator acting on C only. This is the case, for example, when the central and external baths are initially in a product state, ρ B = ρ C (0) ⊗ ρ E (0), and C is sufficiently weakly coupled to an external bath containing enough degrees of freedom for the Born-Markov approximation to hold. The desired correlation function C jk (t, s) is then given by the following expression of multitime averages for evolution under a Markovian ME [62]
Further assuming that the initial state of the central bath is a steady state of the ME, L C ρ C (0) = 0, then directly leads to a correlation function respecting Eq. (A7), and thus to stationary noise.
For stationary noise, Eq. (A5) can be rewritten in the frequency domain as
where the integration kernel Υjk (ω, t) is given by
To describe spin-locking, we replace H S → H S and Q j → -τ x j in Qj (t) = e iHSt Q j e -iHSt , with H S given by Eq. ( 7), leading to Qj (t) = -(e iΩj t τ + j + H.c.). Substituting this into Eq. (A10) above, we obtain terms oscillating at angular frequencies |Ω j -Ω k | and Ω j + Ω k . Assuming that (Ω j + Ω k )t D 1, ∀j, k, where t D is the typical timescale over which qubit observables decay, we neglect terms oscillating with Ω j + Ω k by invoking a secular approximation [40]. This leads to the simplified expression Υjk (ω, t)ρ ≈ (A11)
where ∆ jk ≡ Ω j -Ω k and F (ω) ≡ t 0 ds e iωs is the firstorder fundamental filter function for free-induction decay [63]. Crucially, for finite time, F (ω ± Ω j ) is peaked around ω = ∓Ω j with a width ∼ 1/t and thus acts as a bandpass filter for the noise spectra in the integral over frequencies, Eq. (A9). Assuming that all spectra vary negligibly over this passband, we replace F (ω) by its infinite-time limit in the sense of distributions, lim t→∞ F (ω) = πδ(ω), which gives
Moving back to the spin-locking reference frame, which rotates at the qubit drive frequencies [see the discussion above Eq. ( 4)], and taking Ω 1 = Ω 2 then results in Eq. ( 8) in the main text, which provides the theoretical basis of the spectroscopy method we presented.
We now apply the theory described in Appendix A 1 to the cQED experimental setting considered in the main text (see Fig. 3). In this setting, S comprises a pair of qubits encoded by the two lowest energy levels of transmons. The two qubits are coupled to a central bath consisting of a microwave resonator mode, which is itself subject to damping due to an external environment. We manipulate the qubits and the central bath by irradiating the input port of the common resonator with microwave drives. A continuous microwave drive with strength ε is first applied at frequency ω d near the fundamental resonator mode frequency ω c , in order to bring the mode into a coherent steady state. An additional pair of drives is then employed to initialize the qubits and measure their final state using the dispersive readout. Most importantly for this work, between initialization and readout, continuous spin-locking drives are also applied on each qubit j with constant strength Ω j at frequency ω dj near the bare qubit frequency ω 0 qj . To describe the above experiment theoretically, we consider the sample parameters given in Table I, enabling us to make several simplifying assumptions. First, we assume that the qubit-resonator detunings ∆ j ≡ ω 0 qj -ω c are large compared with: (i) ε and Ω j , enabling us to neglect any cross-talk effect of the resonator drive on the qubits and of the qubit drives on the resonator; and (ii) the qubit-resonator coupling strengths g j , such that the two-qubit dispersive Hamiltonian can be employed [49]. Second, we also assume that the qubit-qubit detuning ω q1 -ω q2 is much larger than the strength of any interaction between the qubits, for example mediated by virtual transitions with the resonator [64]. We then model the two qubits, the resonator mode, and its environment with the Hamiltonian
where χ j ≡ g 2 j /∆ j is the dispersive coupling strength between qubit j and the resonator mode with annihilation and creation operators a and a † , and ω qj = ω 0 qj + χ j is the Lamb-shifted qubit frequency. In Eq. (A14), H CE is the Hamiltonian that describes coupling of the resonator mode with an external environment with free Hamiltonian H E . Typically, this environment is modeled by an ensemble of harmonic oscillators corresponding to the modes of the free electromagnetic field or phonons by taking
, where ω k and v k are the frequency and coupling strength of the environmental mode k, whose excitations are annihilated by the bosonic operator b k [48].
We next move to the frame that rotates at the drive frequencies using the unitary tranformation
This transformation produces terms oscillating at frequencies 2ω d and 2ω dj , which we neglect under the RWA, assuming that ε ω d and Ω j ω dj . In the rotating frame, the total Hamiltonian for the qubits, resonator, and environment is then given by Eq. (A1), in which
and where we have taken
In Eq. (A18), ∆ qj ≡ ω qj -ω dj and ∆ c ≡ ω c -ω d are the detunings of the qubit and resonator drives, respectively, and
. is the interaction-picture coupling between the resonator and its external environment. In addition, to arrive at a noise operator B j that has zero mean in the initial state, we have added and subtracted the term j χ j nσ z j in H S and H SB , respectively, where n ≡ a † a(0) is the average photon number in the initial state of the resonator. In this approach, the mean of the noise is captured by the Stark shift 2χ j n of each qubit j in Eq. (A19) for ∆ qj . Assuming that this shift is measured with sufficient accuracy, we set ∆ qj = -2χ j n, leading to ∆ qj = 0 ∀ j in Eq. (A16). In the spin-locking basis, Eqs. (A16) and (A17) then take the form of Eq. ( 7) we have taken in the derivations of Appendix A 1.
We now follow the steps laid down in Appendix A 1 to obtain the reduced ME for the qubits. In Eq. (A2), we consider the interaction-picture bath operator Bj
] is the evolution operator under H B (t) given in Eq. (A18), with T denoting time ordering. For weak coupling between the resonator mode and an environment consisting of a large number of degrees of freedom, the Born-Markov approximation may be invoked, following which reduced evolution of any central bath (resonator) operator under R B (t) is approximated by solving the Lindblad ME [48,65],
Assuming that the resonator couples to a continuum of environmental modes, the resonator damping rate appearing in Eq. (A21) is given by κ
where D(ω) and v(ω) are the frequency-dependent density of modes and coupling strength of the external environment, respectively. We also assume that the resonator drive is applied since time t 0 -1/κ so that, at t = 0, the resonator mode has reached a steady state defined by ρC (t) = 0, in which noise is stationary. In this limit, Eq. (A20) is easily solved to give n = a † a(0) = ε 2 /[(κ/2) 2 + ∆ 2 c ]. Substituting Eq. (A20) into Eq. (A8), we may evaluate the central bath correlation functions C jk (t, s). In practice, to evaluate these correlators, we find that it is easier to use the equivalent quantum regression theorem [62]. This results in [18,49]
Fourier-transforming this correlation function then yields the shot-noise spectra given by Eq. ( 15).
In Sec. II A and Appendix A 1, we considered an ideal setting in which noise only couples to σ z j , j ∈ {1, 2}. In addition, we took the Rabi frequencies for the two spinlocking drives to be exactly equal, Ω 1 = Ω 2 ≡ Ω. Here, we derive Eq. ( 21), which accounts for both a finite Rabi frequency difference (Ω 1 = Ω 2 ) and noise coupling to the qubits off-axis, via σ x j , j ∈ {1, 2}. Throughout this Appendix, we will assume that the Rabi frequencies Ω 1 and Ω 2 are approximately time-independent within any measurement of decay curves for a given target Rabi frequency Ω, even though Ω 1 and Ω 2 are allowed to undergo small random fluctuations between distinct sets of measurements at different target values of Ω.
To derive the ME, we follow the same steps as in Appendix A 1, but replace
] in the definition of the system-bath Hamiltonian, below Eq. (A1). In contrast with H SB , H x SB includes an additional central-bath operator B x j that couples to the system through Q x j (t). To describe T 1 -effects, we take Q x j (t) to be σ x j in the frame co-rotating with the qubit drives, Q
t) arise. Among these terms, some oscillate with angular frequencies Ω j and ω qj ± Ω j . Since we must have Ω j ω qj , ∀j (which we assumed in Sec. II A to arrive at Eq. ( 5) within the RWA), we may neglect these fast oscillating terms within a secular approximation, as done below Eq. (A10), without further loss of generality compared with Eq. ( 8). In addition, analogously to Eq. (A11), in the additional terms due to noise along x, noise spectra S xx,j (ω) ≡ ∞ -∞ dτ e -iωτ Bx j (τ ) Bx j (0) are probed by filter functions in passbands of width ∼ 1/t about frequencies ±ω qj , ω qj ± Ω j , and -ω qj ± Ω j . Assuming that the noise spectra vary negligibly within these passbands, we take the infinite-time limit as in Appendix A 1, and move back to the spin-locking reference frame rotating at the qubit drive frequencies. We then arrive at the ME
where L ± jk are introduced in Eq. (A13) and
Allowing for Ω 1 and Ω 2 to be distinct, we introduce the Rabi frequency difference δΩ ≡ Ω 1 -Ω 2 and the average Rabi frequency Ω ≡ (Ω 1 + Ω 2 )/2, so that Ω 1 = Ω + δΩ/2 and Ω 2 = Ω -δΩ/2. We then assume that S jk (ω) varies negligibly around ω = Ω over the range of values taken by δΩ, allowing us to approximate S jk (±Ω 1 ) ≈ S jk (±Ω 2 ) ≈ S jk (±Ω), ∀j, k. For the shotnoise spectra produced experimentally, this amounts to assuming that δΩ is negligible in comparison with the width κ of the peak, δΩ κ. This enables us to approximate L - jk + L + jk ≈ L jk in Eq. (A25), with L jk given by Eq. ( 9). Likewise, as in Ref. [23], we assume that S xx,j (ω) varies negligibly over the small (typically, ∼ 1-100 MHz) frequency ranges between ±ω qj -Ω j and ±ω qj + Ω j for all j, and thus take S xx,j (ω qj ± Ω j ) ≈ S xx,j (ω qj ) and S xx,j (-ω qj ± Ω j ) ≈ S xx,j (-ω qj ) in Eq. (A27). Under these approximations, Eq. (A25) then reduces to Eq. ( 21), in which the damping rates due to noise along x are given by Γ 1,j ≡ S xx,j (ω qj ) + S xx,j (-ω qj ) = 1/T
with T (j) 1
the longitudinal relaxation time of qubit j in a free-evolution experiment. The effects of noise along x are then accounted for phenomenologically in the nonlinear regression approach, by using the average T (j) 1 values measured in independent free-evolution experiments.
Though a simultaneous characterization of all noise spectra S jk (ω) and S xx,j (ω) for all qubit axes would be much preferable in principle [45], the above approach, in which the effect of S xx,j (ω) is captured by a single parameter T (j) 1 , has been successfully employed in a singlequbit context [25]. For sufficiently strong photon shot noise, longitudinal qubit relaxation effects are a significant, but weak correction that need only be taken into account at Rabi frequencies for which S jk (ω) is the weakest, i.e., at the tails of the Lorentzian spectra. Though not accounting for T (j) 1 would make the reconstructed S jk (ω) to be significantly off at the tails, we do not expect the details of the additional noise spectra S xx,j (ω) to lead to significant effects beyond corrections obtained with Eq. ( 21).
M-estimators were introduced in the 1960s for robust estimation of a parameter using data whose distribution function is only approximately known: for example, the observations may follow a normal distribution, except for a fraction of them which is affected by experimental error [66]. Following the steps of Ref. [67], in this Appendix we discuss the asymptotic normality of these estimators when the number of observations tends to infinity. This will lead us to the approximate confidence intervals for the spectrum values presented in the main text. No excessive emphasis is put on mathematical rigor; the following steps are valid under suitable regularity conditions that an interested reader may find in Ref. [67].
Given a loss function λ(z) that penalizes deviation between a model and experimental observations, the general form of an M-estimator θ for a p-dimensional parameter vector θ is
where z α is the α-th realization of a random variable Z associated with experimental observations, and whose probability distribution is parameterized by θ. In Sec. II B, the parameter θ is the spectrum vector, θ = S, while in Sec. IV, θ is the spectrum vector supplemented by the Rabi frequency difference δΩ between the qubit drives. Throughout the main text, z α quantifies the relative deviation of observations from their expected value O α θ via z α = (O α -O α θ )/σ α , where O α is the sample mean of M projective measurements performed on the two-qubit system and σ 2 α = var(O α ). Throughout this Appendix, we will assume that all realizations z α are independent and identically distributed (i.i.d.). For example, under ideal conditions (in the absence of any experimental error), and in the asymptotic limit in which each O α is obtained from M → ∞ projective measurements, the central limit theorem implies that each z α is sampled from the standard normal distribution N (0, 1).
Since the value of θ that minimizes the right-hand side of Eq. (B1) may be found by setting derivatives with respect to θ to zero, M-estimators are often alternatively defined by the solution of a set of p estimating equations
To describe the minimization problem of Eq. (B1), above, we set ψ (z α ) ≡ ∂λ(z α )/∂θ . It is also convenient to introduce the d-dimensional column vectors ψ(z α ) and Ψ(θ), whose components are defined as ψ (z α ) and Ψ (θ) ≡ d α=1 ψ (z α ), respectively. With this notation, the estimating equations for the M-estimator of θ may be written succinctly as
To discuss the asymptotic behavior of θ, we Taylorexpand Eq. (B3) about θ = θ * , where θ * is the true value of θ, and truncate to the first order. This gives
where Ψ(θ * ) is the p × p matrix of derivatives of Ψ(θ) at θ = θ * , namely,
By definition of M-estimators, Eq. (B3), the left-hand side of Eq. (B4) must vanish, Ψ( θ) = 0. Assuming that Ψ(θ * ) is invertible then gives
Recall that Ψ(θ * ) and Ψ(θ * ) are ultimately functions of the realizations z α of the random variable describing the deviation of an experimental observation O α from its expected value. We thus investigate the asymptotic behavior of Ψ(θ * ) and Ψ(θ * ) as the number of observations tend to infinity, d → ∞, assuming that the z α 's are i.i.d. Since, by definition, Ψ (θ * ) = d α=1 ψ (z α ), the central limit theorem implies the following convergence in distribution:
where N p (µ, Σ) is the p-dimensional multivariate normal distribution with mean µ and covariance matrix Σ. To arrive at Eq. (B7), we have used the fact that E[ψ] = 0 at θ = θ * , where E(ψ) is the expectation value of ψ over realizations of the random variable Z. This property must be satisfied for any M-estimator since, by the weak law of large numbers, E[ψ] = lim d→∞ 1 d d α=1 ψ(z α ), where d α=1 ψ(z α ) = 0 at θ = θ by Eq. (B3), and lim d→∞ θ = θ * due to asymptotic consistency of Mestimators [67]. Because E[ψ] = 0, the covariance matrix in Eq. (B7) is simply given by cov(ψ) = E(ψψ T ).
In addition, applying the weak law of large numbers to the last expression in Eq. (B5) leads to the following convergence in probability:
where we have introduced the matrix ψ(z α ), whose components are ψk (z α ) ≡ ∂ψ k (z α )/∂θ | θ=θ * . Substituting Eqs. (B7) and (B8) into Eq. (B6), applying Slutsky's theorem [46], and using the linear transformation property cov(Bx) = Bcov(x)B T , where x is a random column vector and B a constant matrix, then implies that θθ * converges in distribution to θθ * → N p 0, Σ θ , (B9)
with the covariance matrix
where θ is evaluated at θ * .
We may now use the definition ψ (z α ) ≡ ∂λ(z α )/∂θ to approximate the covariance matrix of θθ * in terms of derivatives of loss functions for d 1. Evaluating the derivatives using the chain rule and approximating expectation values of functions of observations f (Z) by
where we have introduced the Jacobian matrix with elements J α ≡ ∂z α /∂θ | θ=θ * , and where ∆ and D are diagonal matrices whose non-zero elements are D αα = ∂λ/∂z| z=zα and ∆ αα = ∂ 2 λ/∂z 2 z=zα , respectively. These derivatives are readily obtained for simple loss functions. In particular, for quadratic loss, λ(z) = z 2 /2, we have D αα = z α and ∆ αα = 1, whereas for Huber loss, the following hold: Substituting Eq. (B11) and (B12) into Eq. (B10) gives an approximation of the covariance matrix Σ θ of θθ * for d 1. We then estimate 95% confidence intervals presented in the main text from θ ± 1.96 Σ θ .
For a sufficiently weakly non-linear relationship between z α and the parameters, the term proportional to second-order derivatives in Eq. (B12) may be neglected, leading to the compact expression Σ θ ≈ (J T ∆J) -1 (J T D 2 J)(J T ∆J) -1 .
(B14)
Though we do account for the term containing secondorder derivatives in all confidence-interval calculations presented in the main text for completeness, in all cases studied here, this term is found to be negligible. Since it avoids the computation of second-order derivatives of z α with respect to all possible pairs of θ k and θ , evaluation of the confidence intervals through Eq. (B14) is much more efficient numerically, leading to significant savings in computation time.
we plot an example of decay curve -here, τ z 2 -as a function of evolution time for Ω/2π = 1.848 MHz as an inset in Fig. 7(b). This inset reveals that while the fit using Huber loss closely follows the bulk of the data points, the fit using weighted least squares goes astray, most plausibly due to the influence of outliers in the many decay curves involved in the global regression procedure. As indicated by the dashed vertical lines in Fig. 7(a)-(b), the reconstructed spectra at the corresponding frequencies ω/2π = ±1.848 MHz lie closer to the theoretical value using Huber loss than with weighted least squares, in particular for S 22 (ω).
To verify that such a failure of weighted least squares can indeed arise from outliers, we reproduce the effect with simulated data in Fig. 7(c)-(d). To produce the simulated data used in the spectral estimation, we first calculate the two-qubit density matrix ρ(t) by substituting the spectra given by Eq. ( 15) into the ME defined by Eqs. ( 8)-( 9), which we solve numerically for the shot-noise parameters measured in the experiment (caption of Fig. 6). From the resulting ρ(t), we evaluate probabilities associated with binary outcomes (±1) of projective measurements, and generate M = 2000 simulated outcomes for each observable O r by sampling a Bernoulli distribution. We then evaluate the sample means O α of projective measurement outcomes for each combination of initial state, observable, and evolution time given in Table II, and their corresponding standard deviations. Finally, to emulate experimental error, we use a δ-contaminated model [46], in which each value of O α is assigned a probability 0.1 to be replaced by an outlier, which we model by sampling a uniform probability distribution between -1 and +1.
After constructing this simulated data set, we perform the reconstructions using the procedure explained in Sec. II B. Despite the presence of numerous outliers, Fig. 7(c) shows that our procedure successfully reconstructs all the spectra when using Huber loss. However, in Fig. 7(d), the equivalent reconstructions using weighted least squares are remarkably more noisy, and even involve a spurious positive-frequency component. The discrepancy between outcomes of the two loss functions is significantly more pronounced than in the experiment, as may be expected for the rather pessimistic model of outliers employed here. The inset of Fig. 7(d) displays the outcome of the fitting procedure for a pair of frequencies at which weighted least squares failed to produce accurate spectrum estimates [Ω/2π = 2.120 MHz, dashed vertical lines in Fig. 7(c)-(d)]. Again, while the decay curve fitted using weighted least squares wanders away, Huber loss produces a decay curve that closely follows the bulk of the data.
Master equation for noise along σ z j only