Non-Gaussian noise spectroscopy with a superconducting qubit sensor

TL;DR

This paper demonstrates a quantum control protocol that characterizes non-Gaussian noise affecting superconducting qubits, enabling more accurate noise modeling beyond Gaussian assumptions.

Contribution

It introduces the first experimental validation of non-Gaussian noise spectroscopy using a superconducting qubit sensor, capturing higher-order noise spectra.

Findings

Successfully reconstructs non-Gaussian noise spectra

Validates the protocol with engineered noise sources

Advances spectral estimation tools for quantum devices

Abstract

Accurate characterization of the noise influencing a quantum system of interest has far-reaching implications across quantum science, ranging from microscopic modeling of decoherence dynamics to noise-optimized quantum control. While the assumption that noise obeys Gaussian statistics is commonly employed, noise is generically non-Gaussian in nature. In particular, the Gaussian approximation breaks down whenever a qubit is strongly coupled to discrete noise sources or has a non-linear response to the environmental degrees of freedom. Thus, in order to both scrutinize the applicability of the Gaussian assumption and capture distinctive non-Gaussian signatures, a tool for characterizing non-Gaussian noise is essential. Here, we experimentally validate a quantum control protocol which, in addition to the spectrum, reconstructs the leading higher-order spectrum of engineered non-Gaussian dephasing noise using a superconducting qubit as a sensor. This first experimental demonstration of non-Gaussian noise spectroscopy represents a major step toward demonstrating a complete spectral estimation toolbox for quantum devices.

Tables

Table 1. Table S1: Device parameters.

Parameter	Value
Qubit frequency ${\omega }_0/2\pi$	2.920 GHz
Qubit anharmonicity $A/2\pi$	1.163 GHz
Relaxation time $T_1$	27.0 $\pm$ 2.7 $\mu s$
Spin-echo relaxation time $T_2$	35.9 $\pm$ 4.4 $\mu s$
Free induction decay time $T_2^{\ast }$	12.2 $\pm$ 1.0 $\mu s$
Readout cavity frequency ${\omega }_r/2\pi$	7.348 GHz
Readout cavity linewidth $\kappa /2\pi$	2.548 MHz
Dispersive coupling strength $\chi /2\pi$	0.130 MHz

Table 2. Table S2: Control pulse sequences designed for non-Gaussian spectral estimation.

Seq. Index $p$	Position of $\pi$ pulses [$\mathrm{ns}$]	Repetitions $M$	Filter function at zero frequency $F_p(0,T)$
1	No pulses (free evolution)	1	$\ne$ 0
2	125, 175, 225, 275, 325, 610, 820, 875	10	$\ne$ 0
3	90, 235, 410, 555, 730, 875	10	$\ne$ 0
4	80, 150, 205, 355, 560, 630, 685, 835	10	$\ne$ 0
5	105, 240, 345, 480, 585, 720, 825, 960	10	$\ne$ 0
6	85, 135, 185, 240, 455, 775, 825, 880	10	0
7	130, 180, 285, 335, 475, 765, 870, 960	10	0
8	90, 150, 200, 305, 500, 715, 860, 960	10	0
9	80, 320, 370, 425, 600, 650, 720, 855	10	0
10	205, 310, 360, 545, 645, 725, 850, 960	10	0
11	145, 365, 425, 495, 600, 680, 850, 960	10	0

Equations

$$\chi (t)$$

$$\phi (t)$$

$$B(t)={\beta }_{\Phi }[\Delta \Phi (t){]}^2,\Delta \Phi (t)\equiv \Phi (t)-{\Phi }_0/2,$$

$${\mu }_B$$

$$S(\omega )$$

$$S_2({\omega }_1,{\omega }_2)$$

$${\chi }_p(MT)\approx \frac{M}{T}\sum _{k\in {\mathcal{K}}_1}|F_p(k{\omega }_h,T){|}^{2}S(k{\omega }_h).$$

$${\phi }_p(MT)=-\frac{M}{3!T^2}\sum _{\vec{k}\in {\mathcal{K}}_2}m({\omega }_h\vec{k})\text{Re}[G_p({\omega }_h\vec{k},T)]S_2({\omega }_h\vec{k}),$$

$$\vec{\phi }=\mathbf{A}{\vec{S}}_2,{\mathbf{A}}_{pn}=-\frac{M}{3!T^2}m({\omega }_h{\vec{k}}_n)\text{Re}[G_p({\omega }_h{\vec{k}}_n,T)],$$

$${\vec{S}}_2^{\text{RMLE}}=\underset{{\vec{S}}_2}{\mathrm{argmin}}\left[\frac{1}{2}(\mathbf{A}{\vec{S}}_2-\vec{\phi }{)}^T{\mathbf{\Sigma }}^{-1}(\mathbf{A}{\vec{S}}_2-\vec{\phi })+||\lambda \mathbf{D}{\vec{S}}_2|{|}_2^2\right],$$

$$\Delta \Phi (t)=\sum _{m=1}^{N_h}\left(a_m\cos {\omega }_{m}t+b_m\sin {\omega }_{m}t\right),$$

$$\mathbb{E}[a_m]=\mathbb{E}[b_m]=0,\forall m,\mathbb{E}[a_m{b}_n]=0,\forall m,n,\mathbb{E}[a_m{a}_n]=\mathbb{E}[b_m{b}_n]=0,\forall m\ne n.$$

$$S_{\Phi }(\omega )=\frac{P_0/\pi {\omega }_c}{1+(\omega /{\omega }_c{)}^2}.$$

$$H(t)=\frac{{\omega }_{\mathrm{q}}+B(t)}{2}{\sigma }_z+\epsilon (t)\cos [{\omega }_{\mathrm{d}}t+\theta (t)]{\sigma }_x,$$

$$H_{\mathrm{d}}(t)=\frac{D+B(t)}{2}{\sigma }_z+H_{\mathrm{c}}(t),$$

$$H_{\mathrm{c}}(t)\approx \frac{1}{2}\left[{\eta }_{\mathrm{I}}(t){\sigma }_x+{\eta }_{\mathrm{Q}}(t){\sigma }_y\right],$$

$$H_{\mathrm{T}}(t)=\frac{D+B(t)}{2}{\vec{y}}_p(t)\cdot \vec{\sigma },\vec{\sigma }\equiv ({\sigma }_x,{\sigma }_y,{\sigma }_z),$$

$$⟨{\sigma }_z(t_f)⟩=\mathbb{E}\left\{\mathrm{Tr}\left[R_{\mathrm{T}}^{†}(t_f){\sigma }_z{R}_{\mathrm{T}}(t_f)U_{\mathrm{T}}(t_f){\rho }_0{U}_{\mathrm{T}}^{†}(t_f)\right]\right\},$$

$$⟨{\sigma }_z(t_f)⟩\approx -y_{p,z}(t_f)+(D+{\mu }_B)\left[F_{p,x}(0,t_f)y_{p,y}(t_f)-F_{p,y}(0,t_f)y_{p,x}(t_f)\right],$$

$$⟨{\sigma }_z(t_f)⟩\approx (D+{\mu }_B)T^{\prime },$$

$${\overline{Z}}_j\equiv \frac{1}{N}\sum _{i=1}^N{Z}_{j,i}\sim \mathrm{Normal}\left[(D_j+{\mu }_B)F_{0,y}(0,t_f),\frac{\mathrm{var}({\sigma }_z{)}_j}{N}\right],$$

$$a^{\mathrm{est}}$$

$$b^{\mathrm{est}}$$

$$\mathrm{var}(a^{\mathrm{est}})=\frac{N_D^{-1}{\sum }_j{D}_j^2}{{\sum }_j{\left(D_j-N_D^{-1}{\sum }_k{D}_k\right)}^2}\mathrm{var}(\overline{Z}),\mathrm{var}(b^{\mathrm{est}})=\frac{\mathrm{var}(\overline{Z})}{{\sum }_j{\left(D_j-N_D^{-1}{\sum }_k{D}_k\right)}^2},$$

$$\mathrm{cov}(a^{\mathrm{est}},b^{\mathrm{est}})\equiv \mathbb{E}[(a^{\mathrm{est}}-a)(b^{\mathrm{est}}-b)]=-\frac{N_D^{-1}{\sum }_j{D}_j}{{\sum }_j{\left(D_j-N_D^{-1}{\sum }_k{D}_k\right)}^2}\mathrm{var}(\overline{Z}).$$

$${\tilde{\mu }}_B^{\mathrm{est}}\equiv a^{\mathrm{est}}/b^{\mathrm{est}}.$$

$$\mathrm{var}({\tilde{\mu }}_B^{\mathrm{est}})\approx \frac{b^2\mathrm{var}(a^{\mathrm{est}})+a^2\mathrm{var}(b^{\mathrm{est}})-2ab\mathrm{cov}(a^{\mathrm{est}},b^{\mathrm{est}})}{b^4}.$$

$$\displaystyle\chi_{p}^{\text{est}}(t)=-\frac{1}{2}\text{ln}(\sigma_{x}^{\text{est}\;2}+\sigma_{y}^{\text{est}\;2})=-\frac{1}{2}\text{ln}\big{[}(\tilde{\sigma}_{x}^{\text{est}}+\mu_{x})^{2}+(\tilde{\sigma}_{y}^{\text{est}}+\mu_{y})^{2}\big{]},$$

$${\chi }_p^{\text{est}}(t)\approx -\frac{1}{2}\text{ln}({\mu }_x^2+{\mu }_y^2)-\left(\frac{{\mu }_y}{{\mu }_x^2+{\mu }_y^2}\right){\tilde{\sigma }}_y^{\text{est}}-\left(\frac{{\mu }_x}{{\mu }_x^2+{\mu }_y^2}\right){\tilde{\sigma }}_x^{\text{est}}.$$

$$\mathbb{E}[{\chi }_p^{\text{est}}(t)]=$$

$$\text{var}[{\chi }_p^{\text{est}}(t)]=$$

$$\mathbb{E}[{\chi }_p^{\text{est}}(MT)]\approx \frac{M}{T}\sum _{k=-\infty }^{\infty }|F_p(k{\omega }_h,T){|}^{2}S(k{\omega }_h).$$

$$\displaystyle{\mathbb{E}}[\chi_{p}^{\text{est}}(MT)]\approx\frac{M}{T}\sum_{k=0}^{K-1}\Big{(}\frac{2-\delta_{k,0}}{2}\Big{)}|F_{p}(k\omega_{h},T)|^{2}S(k\omega_{h}).$$

Full text

Non-Gaussian noise spectroscopy with a superconducting qubit sensor

Youngkyu Sung

Research Laboratory of Electronics, Massachusetts Institute of Technology, Cambridge, MA 02139, USA

Department of Electrical Engineering and Computer Science, Massachusetts Institute of Technology, Cambridge, MA 02139, USA

Félix Beaudoin

Department of Physics and Astronomy, Dartmouth College, Hanover, NH 03755, USA

Leigh M. Norris

Department of Physics and Astronomy, Dartmouth College, Hanover, NH 03755, USA

Fei Yan

Research Laboratory of Electronics, Massachusetts Institute of Technology, Cambridge, MA 02139, USA

David K. Kim

MIT Lincoln Laboratory, 244 Wood Street, Lexington, MA 02421, USA

Jack Y. Qiu

Research Laboratory of Electronics, Massachusetts Institute of Technology, Cambridge, MA 02139, USA

Department of Electrical Engineering and Computer Science, Massachusetts Institute of Technology, Cambridge, MA 02139, USA

Uwe von Lüpke

Research Laboratory of Electronics, Massachusetts Institute of Technology, Cambridge, MA 02139, USA

Jonilyn L. Yoder

MIT Lincoln Laboratory, 244 Wood Street, Lexington, MA 02421, USA

Terry P. Orlando

Research Laboratory of Electronics, Massachusetts Institute of Technology, Cambridge, MA 02139, USA

Department of Electrical Engineering and Computer Science, Massachusetts Institute of Technology, Cambridge, MA 02139, USA

Simon Gustavsson

Research Laboratory of Electronics, Massachusetts Institute of Technology, Cambridge, MA 02139, USA

Lorenza Viola lorenza.viola@dartmouth.edu Department of Physics and Astronomy, Dartmouth College, Hanover, NH 03755, USA

William D. Oliver william.oliver@mit.edu Research Laboratory of Electronics, Massachusetts Institute of Technology, Cambridge, MA 02139, USA

Department of Electrical Engineering and Computer Science, Massachusetts Institute of Technology, Cambridge, MA 02139, USA

MIT Lincoln Laboratory, 244 Wood Street, Lexington, MA 02421, USA

Department of Physics, Massachusetts Institute of Technology, Cambridge, MA 02139, USA

**Accurate characterization of the noise influencing a quantum system of interest has far-reaching implications across quantum science, ranging from microscopic modeling of decoherence dynamics to noise-optimized quantum control. While the assumption that noise obeys Gaussian statistics is commonly employed, noise is generically non-Gaussian in nature. In particular, the Gaussian approximation breaks down whenever a qubit is strongly coupled to discrete noise sources or has a non-linear response to the environmental degrees of freedom. Thus, in order to both scrutinize the applicability of the Gaussian assumption and capture distinctive non-Gaussian signatures, a tool for characterizing non-Gaussian noise is essential. Here, we experimentally validate a quantum control protocol which, in addition to the spectrum, reconstructs the leading higher-order spectrum of engineered non-Gaussian dephasing noise using a superconducting qubit as a sensor. This first experimental demonstration of non-Gaussian noise spectroscopy represents a major step toward demonstrating a complete spectral estimation toolbox for quantum devices. **

Introduction

For any dynamical system that evolves in the presence of unwanted disturbances, precise knowledge of the noise spectral features is fundamental for quantitative understanding and prediction of the dynamics under realistic conditions. As a result, spectral estimation techniques have a long tradition and play a central role in classical statistical signal processing [1]. For quantum systems, the importance of precisely characterizing noise effects is further heightened by the challenge of harnessing the practical potential that quantum science and technology applications promise. Such detailed knowledge is key to develop noise-optimized strategies for enhancing quantum coherence and boosting control fidelity in near-term intermediate-scale quantum information processors [2], as well as for overcoming noise effects in quantum metrology [3, 4]. Ultimately, probing the extent and decay of noise correlations will prove crucial in determining the viability of large-scale fault-tolerant quantum computation [5].

Thanks to their exquisite sensitivity to the surrounding environment, qubits driven by external control fields are naturally suited as “spectrometers,” or sensors, of their own noise [6, 7]. Quantum noise spectroscopy (QNS) leverages the fact that open-loop control modulation is akin to shaping the filter function that determines the sensor’s response in frequency space [8, 9, 10, 11, 12] and, in it simplest form, aims to characterize the spectral properties of environmental noise as sensed by a single qubit sensor. By now, QNS protocols employing both pulsed and continuous control modalities have been explored, and experimental implementations have been reported across a wide variety of qubit platforms – including NMR [13], superconducting quantum circuits [14, 15, 16, 17], semiconductor quantum dots [18, 19, 20, 21], diamond nitrogen vacancy centers [22, 23], and trapped ions [24]. Notably, knowledge of the underlying noise spectrum has already enabled unprecedented coherence times to be achieved via tailored error suppression [25].

While the above advances clearly point to the growing significance of spectral estimation in the quantum setting, they all rely on the assumption that the target noise process is Gaussian – that is, one- and two-point correlation functions suffice to fully specify the noise statistical properties. However, the Gaussian assumption need not be justified a priori and it should rather be validated (or falsified) by the QNS protocol itself. A number of realistic scenarios motivate the consideration of non-Gaussian noise regimes. Statistical processes that are responsible for electronic current fluctuations in mesoscopic devices or the 1/$f$ noise ubiquitously encountered in solid-state quantum devices are not Gaussian in general [26]. In superconducting circuits, previous studies have shown that a few two-level defects within Josephson tunnel junctions can interact strongly with the qubit [27, 28, 29, 30, 31], the resulting decoherence dynamics showing marked deviations from Gaussian behavior under both free evolution and dynamical decoupling protocols [7, 32, 33]. More generally, non-Gaussian noise statistics may be expected to arise whenever a qubit is operated outside a linear-response regime, either due to strong coupling to a discrete environment [34] or to a non-linear energy dispersion relationship. The latter feature, which has long been appreciated to influence dephasing behavior at optimal points [35], is common to all state-of-the-art superconducting qubit archetypes [36, 37, 38, 39]. Thus, statistical correlations higher than second order and their corresponding multi-dimensional Fourier transforms must be taken into account for complete characterization. From a signal processing standpoint, this translates into the task of higher-order spectral estimation [40].

In this work, we experimentally demonstrate non-Gaussian QNS by building on the estimation procedure proposed by Norris et al. [41]. While we employ a flux-tunable superconducting qubit as a sensor, our methodology is portable to other physical testbeds in which classical dephasing noise is the dominant decoherence mechanism. We show how non-Gaussianity distinctively modifies the phase evolution of the sensor’s coherence, resulting in an observable signature to which the spectrum (or power spectral density, PSD) is completely insensitive and which is instead encoded in the leading higher-order spectrum, the bispectrum. Unlike the original proposal [41], the QNS protocol we introduce here makes use of a statistically-motivated maximum likelihood approach. This renders the estimation less susceptible to numerical instability, while allowing measurement errors to be incorporated and both the PSD and the bispectrum to be inferred using a single measurement setup. In order to obtain a clean benchmark for our spectral estimation procedure, we engineer a non-Gaussian noise model by injecting Gaussian flux at the sensor’s degeneracy point, resulting in non-Gaussian frequency noise. The noise implementation is validated by verifying the observed power dependence of the leading cumulants against the expected one. Both the reconstructed PSD and the bispectrum are found to be in quantitative agreement with theoretical predictions within error bars.

Results

Non-Gaussian dephasing noise

Before introducing our experimental test bed, we present the general setting to which our analysis is relevant: a qubit sensor evolving under the combined action of non-Gaussian classical dephasing noise and suitably designed sequences of control pulses. By working in an interaction frame with respect to the internal qubit Hamiltonian and the applied control, and letting $\hbar=1$, the controlled open-system Hamiltonian may be written as $H(t)=y_{p}(t)B(t)\sigma_{z}/2,$ where $B(t)$ is a stochastic process describing dephasing noise relative to the qubit’s eigenbasis defined by the Pauli operator $\sigma_{z}$. The control switching function $y_{p}(t)$ accounts for a sequence $p$ of instantaneous $\pi$ rotations about the $x$ or $y$ axis, starting from initial value $y_{p}(0)=+1$ and toggling between $\pm 1$ with every application of a pulse. Under such a pure-dephasing Hamiltonian, the qubit coherence is quantified by the time-dependent expectation value $\langle\sigma_{+}(t)\rangle\equiv\mathrm{e}^{-\chi(t)+i\phi(t)}\langle\sigma_{+}(0)\rangle$, where the influence of the noise is captured by the decay and phase parameters $\chi(t)$ and $\phi(t)$. These parameters may be formally expanded in terms of noise cumulants, $C^{(k)}(t_{1},\ldots,t_{k})$, $k\,\in\,\{1,2,\ldots,\infty\}$, with $\chi(t)$ taking contribution only from even cumulants and $\phi(t)$ only from odd cumulants [41]. Physically, the $k$-th order cumulant is determined by the multi-time correlation functions $\mathbb{E}[B(t_{1}),\ldots,B(t_{j})]$, with $j\leq k$, where $\mathbb{E}[\cdot]$ denotes the ensemble average over noise realizations.

Since the statistical properties of Gaussian noise are entirely determined by one- and two-point correlation functions, cumulants of order $k\geq 3$ vanish identically. By contrast, for non-Gaussian noise, all cumulants can be non-zero in principle. Assuming that noise is stationary, so that the mean of the process, $\mathbb{E}[B(t)]=C^{(1)}(0)\equiv\mu_{B}$ is constant, the phase parameter may be written as $\phi(t)=\mu_{B}F_{p}(0,t)+\varphi(t)$, with the Fourier transform $F_{p}(\omega,t)\equiv\int_{0}^{t}d\mathrm{s}\,\mathrm{e}^{-i\omega s}y_{p}(s)$ being the fundamental filter function (FF) associated to the control [12]. This expression separates the phase due to the noise mean, which arises for both Gaussian and non-Gaussian noise, from a genuinely non-Gaussian phase $\varphi(t)$, which captures the contribution of all odd noise cumulants with $k\geq 3$. For sufficiently small time or noise strength, we can neglect terms of order $k>3$ in the cumulant expansion, leading to

[TABLE]

where $\vec{\omega}\equiv(\omega_{1},\omega_{2})$ and the second and third noise cumulants enter the qubit dynamics through their Fourier transforms: the PSD or spectrum, $S(\omega)\equiv\int_{\mathbb{R}}\mathrm{d}\tau\,\mathrm{e}^{-i\omega\tau}C^{(2)}(0,\tau)$, and the second-order polyspectrum or bispectrum, $S_{2}(\vec{\omega})\equiv\int_{\mathbb{R}^{2}}\mathrm{d}\vec{\tau}\,\mathrm{e}^{-i\vec{\omega}\cdot\vec{\tau}}C^{(3)}(0,\tau_{1},\tau_{2})$, with $\vec{\tau}\equiv(\tau_{1},\tau_{2})$. In the frequency domain, the influence of such spectra is “filtered” by a corresponding generalized FF – in particular, $G_{p}(\vec{\omega},t)\equiv F_{p}(-\omega_{1},t)F_{p}(-\omega_{2},t)F_{p}(\omega_{1}+\omega_{2},t)$ [12]. Since, to leading order, non-Gaussian features arise in our setting from $S_{2}(\vec{\omega})$, non-Gaussianity of a noise process will be detected and characterized through measurements of $\varphi(t)$.

Experimental setup and noise validation

Our circuit QED system [42, 43] contains an engineered flux qubit [44], which is designed to enable fast single-qubit gates with high fidelity at its flux degeneracy point ($F_{g}>99.9\%$; see Supplementary Notes 1 and 2). Single-qubit operations are performed using cosine-shaped microwave pulses, applying an optimal-control technique to suppress leakage to higher levels [45]. Inductive coupling to a local antenna is used to modulate the external flux $\Phi$ threading the qubit loop interrupted by Josephson junctions (Fig. 1a and Fig. 1b). Near the degeneracy (or optimal [35]) point $\Phi=\Phi_{0}/2$, with $\Phi_{0}$ the flux quantum, the $|0\rangle\rightarrow|1\rangle$ transition frequency $\omega_{\mathrm{q}}$ has an approximately quadratic dependence on the external flux $\Phi$ (Fig. 1c). Hence, a sufficiently slow time-dependent external flux $\Phi(t)$ enables adiabatic modulation of the qubit frequency, leading to

[TABLE]

where $\beta_{\Phi}$ is the quadratic coefficient in the dispersion relation between qubit frequency and flux. Crucially, any non-linear function of a Gaussian process leads to non-Gaussian noise. In particular, the quadratic function implemented in Eq. (3) transduces zero-mean Gaussian flux noise into non-Gaussian qubit-frequency noise (Figs. 1c and d). Assuming that the noise is entirely contributed by the applied $\Delta\Phi(t)$, and that $S_{\Phi}(\omega)$ denotes the corresponding PSD, the mean $\mu_{B}$, PSD $S(\omega)$, and bispectrum $S_{2}(\omega_{1},\omega_{2})$ of $B(t)$ are respectively given by

[TABLE]

In the experiment, we choose $S_{\Phi}(\omega)$ to be a zero-mean Lorentzian function, $S_{\Phi}(\omega)=(P_{0}/\pi\omega_{c})/[1+(\omega/\omega_{c})^{2}]$, where $\omega_{c}/2\pi$ $(=0.5~{}\mathrm{MHz})$ and $P_{0}$ denote the cutoff frequency and the power of the applied flux noise, respectively. As is apparent from Eqs. (4)–(6), cumulants of order $k=1$, 2, and 3 are distinguished by their linear, quadratic, and cubic dependence on power, respectively.

We first validate the intended engineered non-Gaussian noise by demonstrating consistency of the measured power dependence of $\chi$ and $\phi$ with the above prediction. The qubit is initialized to the $+y$ axis by applying a $\pi/2$ pulse about $x$ (rotation $R_{x}(\pi/2)$), and Gaussian flux noise is injected while it evolves in the $xy$-plane of the Bloch sphere for time $T$. During this evolution, we apply a Carr-Purcell-Meiboom-Gill (CPMG) sequence consisting of two refocusing $\pi$ pulses about $y$ (Fig. 2a). At the end of this sequence ($t=T$), the effect of the first cumulant of the noise cancels out ($F_{p}(0,T)=0$) and, as a result, the measured phase becomes solely determined by odd cumulants of order $k\geq 3$: $\phi(T)=\varphi(T)$. To estimate both $\phi$ and $\chi$, we measure $\langle\sigma_{x}\rangle$ and $\langle\sigma_{y}\rangle$ by applying appropriate tomography pulses at time $t=T$, before readout in the $\sigma_{z}$-basis.

Figures 2b and 2c show $\chi$ and $\phi$ as a function of injected flux noise power $P_{0}$ for both the experiment (blue triangles) and Monte Carlo simulations accounting for all cumulants of the applied noise (orange squares, see Supplementary Note 5). Substituting Eqs. (5) and (6) into Eqs. (1) and (2), we also plot the resulting ideal weak-power behavior (gray solid) considering only the leading-order cumulants of order two and three for $\chi$ and $\phi$, respectively. For sufficiently small $P_{0}$, these ideal values are in good agreement with data from both experiment and simulation, showing that $\chi$ and $\phi$ obey the quadratic and cubic power dependences that are expected for the square of a Gaussian flux-noise process under the CPMG sequence. In particular, the cubic dependence of $\phi$ at small $P_{0}$ corroborates the presence of a non-zero third-order cumulant, which would not exist for Gaussian noise. Deviations of the simulations and experimental data from the ideal behavior at large $P_{0}$ are attributable to the contribution of cumulants of order $k>3$. The quantitative agreement between theory, experiment, and simulation observed at low power demonstrates our capability to produce and sense engineered noise that dominates over native one over the relevant parameter regime and exhibits well controlled cumulants, a necessary first step in the experimental validation of non-Gaussian QNS.

Non-Gaussian noise spectroscopy

Having established that $\chi$ and $\phi$ follow their expected behavior, we move on to fully characterizing the first three cumulants of our engineered noise source by measuring its mean, PSD, and bispectrum. Since the noise mean, $\mu_{B}$, manifests itself through a qubit-frequency shift, it can be measured from a simple parameter estimation scheme based on Ramsey interferometry. By contrast, we aim to perform a non-parametric estimation of both the PSD and bispectrum, that is, to reconstruct them at a set of discrete points in frequency space without assuming a prior functional form. Figure 3 illustrates our protocol for simultaneous estimation of the PSD and bispectrum, in which filter design – the selection of pulse times in a control sequence so that the corresponding FF has a particular shape – is instrumental. Building on Ref. [13], applying $M\gg 1$ repetitions of a “base” pulse sequence $p\in\{1,2,\cdots,P\}$, with duration $T$, shapes the FF $|F_{p}(\omega,MT)|^{2}$ into a frequency comb with narrow teeth probing $S(\omega)$ at harmonics $k\omega_{h}$, with $k$ an integer and $\omega_{h}\equiv 2\pi/T$ (Figs. 3a and b). This result generalizes to filters relevant to higher-order spectra [12]: under sequence repetition, $G_{p}(\vec{\omega},MT)$ becomes a two-dimensional (2D) “hyper-comb” with teeth probing $S_{2}(\vec{\omega})$ at $\vec{\omega}\,\in\,\{\vec{k}\omega_{h}\}$, where $\vec{k}\equiv(k_{1},k_{2})$ with $k_{1}$ and $k_{2}$ integers (Fig. 3d).

For both the PSD and bispectrum, distinct pulse sequences have the effect of giving different weights to the comb teeth, granting access to complementary information about $S(k\omega_{h})$ and $S_{2}(\vec{k}\omega_{h})$, enabling their reconstruction. More specifically, in both cases the basic steps of our protocol consist of (i) applying a set of sufficiently distinct pulse sequences $p$ (Fig. 3a); (ii) measuring the corresponding decay and phase parameters; and (iii) solving the resulting systems of linear equations, which give $\chi_{p}(MT)$ and $\varphi_{p}(MT)$ as a function of $S(k\omega_{h})$ and $S_{2}(\vec{k}\omega_{h})$. Since classical noise has a spectrum with even symmetry, $S(\omega)=S(-\omega)$, the PSD is specified across all frequency space by its values at positive frequencies. Likewise, the bispectrum is completely specified by its values over a subspace ${\mathcal{D}}_{2}$ known as the principal domain [46, 41], illustrated in Fig. 3d. Reconstructing the bispectrum over ${\cal D}_{2}$ and exploiting the symmetries that $S_{2}(\vec{\omega})$ exhibits (shown in Fig. 3c) thus suffices to retrieve the bispectrum over the whole relevant frequency domain.

Figure 4 presents experimental results for determining the mean and PSD, which suffice to characterize the noise process in the Gaussian approximation. To measure $\mu_{B}$ by Ramsey interferometry, we apply a pair of $\pi/2$ pulses with a drive at frequency $\omega_{\mathrm{d}}$, first about $x$ at time $t=0$ ($R_{x}(\pi/2)$), and then about $y$ at time $t=T$ ($R_{y}(\pi/2)$). We choose a pulse interval $T=50$ ns, which is short enough for cumulants of order higher than one to be negligible, but long enough to avoid pulse overlap. The qubit polarization at time $t_{f}$ after the two pulses is then $\langle\sigma_{z}(t_{f})\rangle\approx(D+\mu_{B})T^{\prime}$, where $D\equiv\omega_{\mathrm{q}}-\omega_{\mathrm{d}}$ is the drive detuning, and $T^{\prime}$ is an effective time interval that accounts for the finite-width pulse shape (see Supplementary Note 6). Thus, plotting $\langle\sigma_{z}(t_{f})\rangle$ as a function of $D$ produces a straight line whose $x$-intercept is $-\mu_{B}$, leading to an estimate that is insensitive to the pulse shape to first order in the cumulant expansion. Figure 4a presents data for measurements of $\langle\sigma_{z}(t_{f})\rangle$, and shows how we isolate the contribution of the engineered noise source by performing the sequence with (blue data set) and without (black data set) applied noise. The mean of the engineered noise is estimated by subtracting the $x$-intercepts of the straight lines that are fitted to each data set. Performing these fits under the conditional normal model of linear regression (see Supplementary Note 6) yields the estimate $\mu_{B}^{\mathrm{est}}/2\pi=127.1\pm 7.56$ kHz, where the uncertainty corresponds to the 95% confidence interval calculated from the asymptotic normal distribution of qubit polarization.

To estimate the PSD by the comb approach outlined above, we use both a period of free evolution ($p=1$) and $M=10$ repetitions of base sequences $p=2,\ldots,11$ illustrated in Fig. 3a (see Supplementary Note 4 for the actual pulse times). For $M\gg 1$, the FF entering the decay constant in Eq. (1) becomes approximately $|F_{p}(\omega,MT)|^{2}\approx\frac{M}{T}|F_{p}(\omega,T)|^{2}\sum_{k=-\infty}^{\infty}\delta(\omega-k\omega_{h})$, which enables us to sample the PSD at the harmonic frequencies in terms of the (known) control FFs,

[TABLE]

Here, we have used the even symmetry of the PSD, and the high-frequency decay of the PSD and FFs to truncate the comb to a finite set of positive harmonics, ${\cal K}_{1}\equiv\{0,\ldots,K-1\}$. Rather than solving the above linear system by matrix inversion as in Ref. [13], we employ a statistically-motivated maximum likelihood estimate (MLE), which takes experimental error into account (see Supplementary Note 7). Using measurements of $\chi_{p}(MT)$ for each of the same $P=11$ control sequences to be used for the bispectrum estimation, we find a well-conditioned system for $K=8$.

Figure 4b compares the experimentally estimated PSD at the $K=8$ harmonics (blue triangles) with the ideal PSD obtained from Eq. (5) for our engineered noise (solid gray line) and Monte Carlo simulations of the QNS protocol (orange squares). The experimental and simulated estimates of the PSD are plotted along with 95$\%$ confidence intervals obtained from the asymptotic normal distribution of the decay constants. Figure 4c shows the experimental and simulated values of $\chi_{p}(MT)$ that were used as input for the reconstructions, along with ideal values obtained by substituting Eq. (5) into Eq. (1) and approximating the FF by the ideal (infinite) comb as given above. The PSD is slightly underestimated at zero frequency in both the experiment and Monte Carlo simulation since the FF of sequence $p=1$ (a 1 $\upmu\mathrm{s}$-long free induction decay) is comparable in bandwidth to the PSD, whereas the reconstruction procedure assumes the PSD is sampled by infinitely narrow FFs. The disagreement of the experimental and simulated $\chi_{p}(MT)$ for $p=1$ with the ideal value is also explained by the non-negligible bandwidth of the FF (Fig. 4c). Apart from these well-understood discrepancies at $\omega=0$, the quantitative agreement of the experimental reconstruction with simulations and ideal values is remarkable, which demonstrates that our protocol is able to reliably characterize Gaussian features of the applied noise.

We are now in a position to present our key result: the reconstruction of the noise bispectrum. As anticipated, this entails a higher-dimensional analogue of the comb-based approach used for the PSD. We estimate the non-Gaussian phase given in Eq. (2) by subtracting the contribution of the noise mean from the total measured phase, $\varphi_{p}(MT)=\phi_{p}(MT)-\mu_{B}F_{p}(0,MT)$, where we replace $\mu_{B}$ by $\mu_{B}^{\mathrm{est}}$ experimentally determined above. After $M\gg 1$ repetitions of sequence $p$, the FF becomes a 2D comb (Fig. 3d), and the non-Gaussian phase becomes a sampling of the bispectrum at the harmonics $\vec{k}\omega_{h}$, that is, $\varphi_{p}(MT)=-\frac{M}{3!T^{2}}\sum_{\vec{k}\in\mathbb{Z}^{2}}G_{p}(\omega_{h}\vec{k},T)S_{2}(\omega_{h}\vec{k}\,).$ Since both the filter and bispectrum decay at high frequencies, we can truncate this sum to a finite number of $\vec{k}=(k_{1},k_{2})$. As the bispectrum is completely specified by its values on the principal domain, we may further restrict our consideration to a subset of harmonics, $\mathcal{K}_{2}\equiv\{\vec{k}_{1},\ldots,\vec{k}_{N}\}\subset{\mathcal{D}}_{2}$ (Fig. 3c). The non-Gaussian phase then becomes

[TABLE]

where the multiplicity $m(\omega_{h}\vec{k})$ accounts for the number of points equivalent to $S_{2}(\omega_{h}\vec{k}\,)$ by the symmetry properties of the bispectrum. Also on account of these symmetries, the imaginary component of $G_{p}(\omega_{h}\vec{k},T)$ cancels when the sum is restricted to $\mathcal{D}_{2}$ (see Supplementary Note 8).

By measuring the non-Gaussian phase for $P\geq N$ different control sequences, we can construct a vector $\vec{\varphi}=[\varphi_{1}(MT),\ldots,\varphi_{P}(MT)]^{T}$ and a linear system of the form

[TABLE]

where $\vec{S}_{2}=[S_{2}(\omega_{h}\vec{k}_{1}),\ldots,S_{2}(\omega_{h}\vec{k}_{N})]^{T}$ contains the bispectrum at the harmonics in $\mathcal{K}_{2}$ and $\mathbf{A}$ is a $P\times N$ reconstruction matrix. The simplest way to estimate the bispectrum from this linear system is the least-squares estimate employed in Ref. [41], involving the (pseudo-)inverse of the reconstruction matrix, $\vec{S}^{\mathrm{est}}_{2}=\mathbf{A}^{-1}\vec{\varphi}$. As in the case of PSD estimation, a potential drawback of this inversion-based approach is numerical instability stemming from an ill-conditioned $\mathbf{A}$, which occurs when the FFs have a high degree of spectral overlap. Since ill-conditioning makes the least-squares estimate sensitive to even small errors in the measured phases, we again utilize a maximum-likelihood approach with optional regularization to further increase stability (see Supplementary Note 8). From the asymptotic Gaussian distribution of the measurement outcomes of $\vec{\varphi}$, the regularized maximum-likelihood estimate (RMLE) is found as

[TABLE]

where $\!|\!|\,\cdot\,\!|\!|_{2}$ denotes the $L_{2}$-norm and $\lambda\geq 0$ parametrizes the strength of the regularization  [47]. Due to its dependence on the covariance matrix $\mathbf{\Sigma}$, the RMLE down-weights phase measurements with larger error. Numerical stability is increased by the regularizer $|\!|\lambda\mathbf{D}\vec{S}_{2}|\!|_{2}^{2}$, which acts as an effective constraint. When the smoothing matrix $\mathbf{D}$ is proportional to $\mathbf{I}$, the regularizer reduces to the well-known Tikhonov (or $L_{2}$) form. Since the numerical stability afforded by regularization comes at the cost of additional bias, choosing the regularization strength is a nontrivial task. In Supplementary Note 8, we detail how we have selected $\lambda$ based on the so-called “L-curve criterion”. Interestingly, since $\mathbf{A}$ is sufficiently well-conditioned for the sequences we have chosen, we find that regularization gives negligible benefit. Accordingly, we use $\lambda=0$ (which recovers standard MLE) in our experimental reconstructions.

Figure 5a compares the results of the non-Gaussian spectral estimation for the harmonics in the principal domain for the experiment (blue triangles) with both the ideal bispectrum obtained from Eq. (6) (gray circles) and from Monte Carlo simulations (orange squares). To estimate the experimental bispectrum, we input the measured data for $\vec{\varphi}$ and $\mathbf{\Sigma}$ shown in Fig. 5b into $\vec{S}_{2}^{\,\text{RMLE}}$ given by Eq. (10). The ideal values of $\varphi_{p}$, also shown in Fig. 5b, are obtained by substituting Eq. (6) into Eq. (2). We further display 3D representations of the full bispectra, obtained by applying relevant symmetries to the data on ${\cal D}_{2}$, for the ideal (Fig. 5c) and experimental (Fig. 5d) cases, respectively. Ignoring error bars, the reconstructed bispectrum appears to be an over-estimate with respect to the ideal one. This error may be attributed to noise during the finite-duration control pulses used in the experiment, leading to effective pulse infidelity. Upon taking the error bars in Fig. 5a into consideration, however, the ideal and simulated values of the bispectrum lie within the 95$\%$ confidence intervals of the experimental reconstruction, suggesting that this estimation error is statistically insignificant and thus successfully extending the validation of our QNS protocol to the leading non-Gaussian noise cumulant.

Although the theoretical bispectrum falls within the 95$\%$ confidence interval of the estimate, reducing the magnitude of uncertainties is clearly necessary to push the application of non-Gaussian QNS to uncontrollable native noise, whose strength may be comparatively weak. We note that the spectral characterization of the non-Gaussian noise process engineered in this experiment requires an extremely precise estimation of $\mu_{B}$. Since reconstructions of the bispectrum are obtained using $\varphi_{p}(MT)=\phi_{p}(MT)-\mu_{B}F_{p}(0,MT)$, the uncertainty in $\mu_{B}^{\mathrm{est}}$ propagates to $\varphi_{p}(MT)$ when $p$ has zero filter order, i.e. $F_{p}(0,MT)\neq 0$. These sequences play a crucial role in estimating the bispectrum at the “zero points”, grid points $(\omega_{1},\omega_{2})$ with $\omega_{1}=0$ or $\omega_{2}=0$. Since $\mu_{B}$ is much larger than the third cumulant for the current noise process, even a small relative uncertainty in $\mu_{B}^{\mathrm{est}}$ can lead to greater error in the bispectrum estimate at the zero points, as the error bars in Fig. 5a attest.

Discussion

In summary, we reported the first experimental demonstration of high-order spectral estimation in a quantum system. By producing and sensing engineered noise with well-controlled cumulants, we were able to successfully validate a spectroscopy protocol that reconstructs both the power spectral density and the bispectrum of non-Gaussian dephasing noise. Our theory and experimental demonstration lay the groundwork for future research aiming at complete spectral characterization of realistic non-Gaussian noise environments in quantum devices and materials. Theoretically, we expect that the regularized maximum-likelihood estimation approach to quantum noise spectroscopy we invoked here will prove crucial to ensure stable spectral reconstructions in more general settings. Devising alternative estimation protocols based on optimally band-limited control modulation and multitaper techniques [48] appears especially compelling, in view of recent advances in the Gaussian regime [24, 49]. We believe that obtaining a complete spectral characterization will ultimately provide deeper insight into the physics and interplay of different microscopic noise mechanisms, including non-classical non-Gaussian noise, as possibly arising from photon-number mediated non-linear couplings [50].

Data availability

The data that support the findings of this study may be made available from the corresponding authors upon request and with the permission of the US Government sponsors who funded the work.

Code availability

The code used for the analyses may be made available from the corresponding authors upon request and with the permission of the US Government sponsors who funded the work.

Acknowledgement

It is a pleasure to thank K. Harrabi, M. Kjaergaard, P. Krantz, G. A. Paz-Silva, J.I.J. Wang, and R. Winik for insightful discussions, and M. Pulido for generous assistance. We thank H. Bethany for the SEM image of the device. This research was funded by the U.S. Army Research Office grant No. W911NF-14-1-0682 (to L.V. and W.D.O.); and by the Department of Defense via MIT Lincoln Laboratory under Air Force Contract No. FA8721-05-C-0002 (to W.D.O.). Y.S. and F. B. acknowledge support from the Korea Foundation for Advanced Studies and from the Fonds de Recherche du Québec – Nature et Technologies, respectively. The views and conclusions contained herein are those of the authors and should not be interpreted as necessarily representing the official policies or endorsements, either expressed or implied, of the U.S. Government.

Author contribution

Y.S., F.Y., and S.G. performed the experiments. F.B. and Y.S. carried out numerical simulations and analyzed the data and L.V., L.M.N., S.G. and W.D.O. provided feedback. L.M.N., F.B., and L.V. designed the pulse sequences and developed the estimation protocol used in the experiment. F.Y. and S.G. designed the device and D.K.K. and J.L.Y. fabricated it. J.Y.Q. and U.L. provided experimental assistance. Y.S., F.B., L.M.N., and L.V. wrote the manuscript with feedback from all authors. L.V., S.G., T.P.O., and W.D.O. supervised the project.

Competing Interests

The authors declare no competing interests.

Supplementary Information

1 Additional experimental and computational detail

1.1 Device parameters and fabrication of the qubit sensor

Device parameters are summarized in Table S1.

The device was fabricated in the same way as in Ref. [37]. It is a generalized version of the capacitively shunted flux qubit [37], comprising a capacitively shunted small-area junction in parallel is a series array of $N$ junctions. In the capacitively shunted flux qubit of Ref. [37], this array comprised $N=2$ junctions. Here, the number of array junctions is $N=8$, far fewer than used in the fluxonium regime of operation [S1]. The area of each Josephson junction forming the array is identical and designed to be $0.2\times 1.2~{}\upmu\textrm{m}^{2}$. The left junction in Fig. 1a is smaller in area by a factor of 8 than the right junction ($\alpha=1/8$). The critical current density $J_{\mathrm{c}}$ is measured to be $0.60\pm 0.01~{}\upmu\mathrm{A}/\upmu\mathrm{m}^{2}$ and the shunt capacitance $C_{\mathrm{sh}}$ is designed to be 20 $\mathrm{fF}$.

1.2 Randomized benchmarking of single-qubit gates

We characterized an average error rate of single-qubit gates by performing Clifford randomized benchmarking [S2] (Fig. S1). As mentioned in the main text, single-qubit operations are performed using cosine-shaped microwave pulses, applying a quadrature correction (DRAG [45]) to minimize unwanted phase evolution and leakage due to the presence of higher levels.

1.3 Measurement setup

1.3.1 Cryogenic setup

The experiments were performed using a Leiden CF-450 dilution refrigerator, capable of cooling to a base temperature of 15 mK. The samples were magnetically shielded with a superconducting can surrounded by a Cryoperm-10 cylinder. The schematic of the cryogenic circuitry is shown in Figure S2. There are two RF lines for the input and the output of the samples; applying microwave readout tone and measuring the transmission of sample respectively. Thermal noise from room temperature on the RF drive lines is attenuated with 20 dB at the 3 K stage, followed by 6 dB at still, and 26 dB at the 20 mK stage. All attenuators in the cryogenic samples are made by XMA. Note that there is one additional RF line for pumping the Josephson traveling wave parametric amplifier (JTWPA) [S3] used as a first-stage pre-amplifier to amplify the readout signal at base temperature. The effective noise temperature is determined primarily by the JTWPA, with a total system noise measured to be about 600 mK. To avoid any back-action of the pump-signal from TWPA, we added a microwave isolator between the samples and the TWPA. On the RF output line, there is a high-electron mobility transistor (HEMT) amplifier (Cryo-1-12 SN508D) at the 3 K stage. Two microwave isolators allow for the signal to pass through to the amplifier without being attenuated, while taking all the reflected noise off of the amplifier and dumping it in a 50 $\Omega$ termination instead of reaching the sample.

There are two additional lines for qubit flux bias: one is for DC flux bias, which is applied globally through a coil installed in the device package, the other is to apply magnetic flux to the qubit thorough a local antenna. The primary requirement of the DC flux bias line is the ability to tune through at least a single flux quantum on the SQUID of the qubit with high precision and low noise. The local flux bias line is attenuated by 20 dB at the 3 K stage, 6 dB at the still, 20 dB at the 20 mK stage to remove excess thermal photons from higher-temperature stages.

1.3.2 Room temperature control

Outside of the cryostat, we have all of the control components which allow us to apply microwave signals that address the cavity and the qubits, as well as the components necessary to resolve the readout signal. All the signals are added using microwave power splitters (Marki PD0R413) used in reverse. Direct digital synthesis of the qubit signals is performed using a high-speed arbitrary waveform generator (AWG Keysight M8195A), which has a 65 GS/sec sampling rate and sufficient bandwidth for this purpose. The output line is further amplified outside of the cryostat with an amplifier (MITEQ AMF-5D-00101200-23-10P) with a quoted noise figure of 2.3 dB, and a preamplifier (Stanford Research SR445A). A detailed schematic is given in Fig S2. We use an IQ demodulation technique to mix down the signal entering the RF port with a reference signal detuned by 40 MHz applied to the LO port. This results in down-converted signals to 40 MHz using a mixer. All components are frequency-locked via a common SRS rubidium frequency standard (10 MHz).

1.3.3 Pulse generation

Qubit control pulse generation is performed via a Keysight M8195A AWG. The pulses are programmed in Labber and then uploaded to the Keysight M8195A.

1.4 Generation of engineered noise

To synthesize a zero-mean flux-noise process $\Delta\Phi(t)\equiv\Phi(t)-\Phi_{0}/2$ with a given PSD, we use an AWG to produce sample waveforms consisting of $N_{h}$ harmonics

[TABLE]

where $\omega_{m}=2\pi m/T_{0}$, and with $2\pi/T_{0}$ the fundamental angular frequency. The Fourier coefficients $a_{m}$ and $b_{m}$ are random variables with

[TABLE]

Further taking $a_{m}$ and $b_{m}$ to have normal (Gaussian) distributions with variance $\sigma_{m}^{2}=2S_{\Phi}(\omega_{m})/T_{0}$, the waveforms $\Delta\Phi(t)$ become a discrete approximation of a Gaussian stochastic process with a frequency-domain PSD $S_{\Phi}(\omega)$.

In all experiments presented in the main text, we consider

[TABLE]

with $P_{0}$ the noise power and $\omega_{c}/2\pi=0.5$ MHz. For the experiment presented in Fig. 2 of the main text, to produce a discrete approximation of a noise process with this spectrum, we take $T_{0}=20$ $\upmu$s and $N_{h}=10^{3}$, corresponding to harmonics separated by the fundamental frequency $1/T_{0}=50$ kHz with a high-frequency cutoff at $\omega_{N_{h}}/2\pi=N_{h}/T_{0}=50$ MHz. For the experiment presented in Figs. 4 and 5 of the main text, we take $T_{0}=200$ $\upmu$s and $N_{h}=10^{4}$, yielding harmonics separated by $1/T_{0}=5$ kHz with the same high-frequency cutoff at $\omega_{N_{h}}/2\pi=N_{h}/T_{0}=50$ MHz. A new waveform is produced by the AWG for each measurement of a Pauli operator performed on the qubit to ensure statistical independence of the samples of the stochastic process, leading to a total number of noise samples of 40,000 for the experiment presented in Fig. 2, and 80,000 for Figs 4 and 5. The waveforms have a duration of $1.25$ $\upmu$s for the experiment presented in Fig. 2 and $20$ $\upmu$s for Figs. 4 and 5.

1.5 Control pulse sequences

The set of control pulse sequences designed for reconstructing the bispectrum are summarized in Table S2 and visualized in Fig. S3. Note that all control pulse sequences start and end with a $\pi/2$ pulse for the purposes of state preparation and tomography.

1.6 Monte Carlo simulations

For the Monte Carlo simulations that are presented in the main text, we consider a single qubit controlled via a microwave drive at angular frequency $\omega_{\mathrm{d}}$, which is used to apply pulses about $\sigma_{x}$ and $\sigma_{y}$. In contrast with the main text, here we do not assume that these pulses are instantaneous. To describe the time-evolution of the qubit under the combined action of these finite-width pulses and classical noise described by the process $B(t)$, we consider the Hamiltonian in the lab frame,

[TABLE]

where $\omega_{\mathrm{q}}$ is the qubit angular frequency, and $\varepsilon(t)$ and $\theta(t)$ are the drive amplitude and phase, respectively.

We next move to the frame that rotates at the drive frequency by applying the unitary transformation $R_{\mathrm{d}}(t)=\exp(-i\omega_{\mathrm{d}}t\,\sigma_{z}/2)$, leading to the Hamiltonian $H_{\mathrm{d}}(t)=R^{\dagger}_{\mathrm{d}}(t)H(t)R_{\mathrm{d}}(t)-iR^{\dagger}_{\mathrm{d}}(t)\dot{R}_{\mathrm{d}}(t)$ in the rotating frame. This gives

[TABLE]

where $D\equiv\omega_{\mathrm{q}}-\omega_{\mathrm{d}}$ is the drive detuning and $H_{\mathrm{c}}(t)\equiv\varepsilon(t)\cos[\omega_{\mathrm{d}}t+\theta(t)][\sigma_{+}\exp(i\omega_{\mathrm{d}}t)+\mathrm{H.c.}]$ is the control Hamiltonian. Assuming that $\varepsilon(t)\ll\omega_{\mathrm{d}}$ and that $\theta(t)$ varies on a timescale much longer than $2\pi/\omega_{\mathrm{d}}$ allows us to invoke the rotating-wave approximation, under which terms oscillating like $\exp\{\pm[2i\omega_{\mathrm{d}}t+\theta(t)]\}$ are neglected. The resulting control Hamiltonian may be simplified as

[TABLE]

where $\eta_{\mathrm{I}}(t)\equiv\varepsilon(t)\cos\theta(t)$ and $\eta_{\mathrm{Q}}(t)\equiv\varepsilon(t)\sin\theta(t)$ are the envelopes of the in-phase and quadrature components of the microwave control signal, respectively. Equation (16) then describes finite-width pulses about $x$ or $y$ axes.

We perform Monte Carlo simulations by solving the time-dependent Schrödinger equation associated with the Hamiltonian $H_{\mathrm{d}}(t)$ given by Eq. (15), under the control Hamiltonian of Eq. (16). The drive detuning $D$ is set to zero in the simulations. The envelope of each control pulse is cosine with total pulse duration 11 $\mathrm{ns}$ (see inset of Fig. S1).

In our Monte Carlo simulations, we account for non-Gaussian noise by setting $B(t)\equiv\beta_{\Phi}\Delta\Phi(t)^{2}$ in Eq. (15), and producing random samples of Gaussian flux-noise $\Delta\Phi(t)$ through the approach described in Section 1.4 of this Supplement. Because this approach relies on an exact solution of the qubit evolution under the noise samples, it is equivalent to accounting for all the terms in the cumulant expansion. The fundamental frequency $1/T_{0}$ and the number of harmonics $N_{h}$ involved in the Fourier-series representation of the noise process are the same as in Section 1.4. To perform the simulations, we generate 100,000 noise samples for the data presented in Fig. 2 of the main text, and 80,000 noise samples for Figs. 4 and 5.

2 Additional theoretical detail on estimation procedure

2.1 Estimation of the noise mean

2.1.1 Ramsey estimation protocol

To measure the noise mean $\mu_{B}$, we use the Ramsey sequence illustrated in Fig. S4a. In this sequence, $\pi/2$ pulses about $\sigma_{x}$ and $\sigma_{y}$ are applied at times $t=0$ and $t=T$, respectively, followed by a measurement of the qubit in the $\sigma_{z}$ eigenbasis at time $t_{f}=T+\Delta T$, where $\Delta T$ is a buffer time much shorter than $T$, but longer than the pulse width. To lay down the theoretical basis of the procedure, we start from the rotating-frame Hamiltonian $H_{\mathrm{d}}(t)$ introduced in Eq. (15), above. To describe the effects of control with finite-width pulses, it is useful to move to the toggling frame using the unitary transformation $R_{\mathrm{T}}(t)=\mathcal{T}\exp[-i\int_{0}^{t}ds\,H_{\mathrm{c}}(s)]$, with $\mathcal{T}$ the time-ordering operator, and where $H_{c}(t)$ is given by Eq. (16). In this toggling frame, the Hamiltonian is

[TABLE]

where $\vec{y}_{p}(t)$ has components $y_{p,a}(t)\equiv\frac{1}{2}\mathrm{Tr}[R_{\mathrm{T}}^{\dagger}(t)\sigma_{z}R_{\mathrm{T}}(t)\sigma_{a}]$, with $a\,\in\,\{x,y,z\}$, for pulse sequence $p$. Remark that a pulse sequence consisting of instantaneous $\pi$ pulses (instead of the Ramsey sequence considered here) would result in $\vec{y}_{p}(t)=[0,0,y_{p}(t)]$, where $y_{p}(t)$ is the switching function used in the main text.

Moving back to the lab frame, the expectation of the $z$ component of the qubit polarization after the pulse sequence is

[TABLE]

where $U_{\mathrm{T}}(t)\equiv\mathcal{T}\exp[-i\int_{0}^{t}ds\,H_{\mathrm{T}}(s)]$ is the time-evolution operator in the toggling frame, and $\rho_{0}$ is the initial qubit density matrix, before application of the pulses. We evaluate $\langle\sigma_{z}(t_{f})\rangle$ perturbatively by performing a Dyson expansion of $U_{\mathrm{T}}(t)$. Assuming that $t_{f}$ is sufficiently short and that $D+B(t)$ is sufficiently small, we truncate this expansion to the first order in $D+B(t)$. Upon substitution into Eq. (18), this approach is equivalent to neglecting any contribution of cumulants of the noise beyond order 1. Taking $\rho_{0}\equiv|0\rangle\langle 0|$, with $|0\rangle$ the eigenstate of $\sigma_{z}$ with eigenvalue $-1$, then yields

[TABLE]

where we have introduced the filter functions $F_{p,a}(\omega,t)\equiv\int_{0}^{t}ds\,\mathrm{e}^{-i\omega s}y_{p,a}(s)$, with $a\,\in\,\{x,y,z\}$. For the pulse sequence illustrated in Fig. S4a (which we label by $p=0$), neglecting any overlap between the pulses, it is straightforward to show that $\vec{y}_{p}(t_{f})\equiv\vec{y}_{0}(t_{f})=(-1,0,0)$ for $t_{f}=T+\Delta t$, leading to

[TABLE]

where $T^{\prime}\equiv F_{0,y}(0,t_{f})=\int_{0}^{t_{f}}ds\,y_{0,y}(s)$ can be viewed as an effective pulse interval accounting for the shape of the pulses. For instantaneous pulses, $T^{\prime}=T$.

According to Eq. (20), setting $D=0$, the first noise cumulant $C^{(1)}(0)\equiv\mu_{B}$ may be estimated simply by measuring $\langle\sigma_{z}(t_{f})\rangle$ and evaluating $\langle\sigma_{z}(t_{f})\rangle/F_{0,y}(0,t_{f})$. However, this approach requires accurate knowledge of $F_{0,y}(0,t_{f})$, and thus of the shape of the control pulses. For the short Ramsey sequences required here to neglect cumulants of order higher than one in Eqs. (19)-(20), the estimate of $\mu_{B}$ then becomes excessively sensitive to distortions of the pulse envelopes that occur in practice. As a result, estimates of $\mu_{B}$ become significantly biased. Since the bispectrum estimation technique that will be discussed in Sec. 2.3 requires precise knowledge of $\mu_{B}$, this bias precludes accurate non-Gaussian QNS.

2.1.2 Robust implementation via linear regression

Crucially, the vulnerability of the above Ramsey scheme to pulse-width effects can be alleviated by estimating $\mu_{B}$ with a linear regression procedure. Indeed, according to Eq. (20), plotting $\langle\sigma_{z}(t_{f})\rangle$ as a function of the detuning $D$ results in a straight line that intersects with the abscissa at $D=-\mu_{B}$ (Fig. S4b), irrespective of $F_{0,y}(0,t_{f})$. Therefore, measuring $\langle\sigma_{z}(t_{f})\rangle$ as a function of $D$ and performing a linear fit of the resulting data leads to an estimate of $\mu_{B}$ that is insensitive to the pulse shape to the first order in $D+B(t)$.

To apply this idea to our experimental data, we now explicitly construct an estimator of $\mu_{B}$ based on linear regression. For each drive detuning $D_{j}$, with $j\,\in\,\{1,2,\ldots,N_{D}\}$, we consider $N$ projective measurements of $\sigma_{z}$ yielding outcomes $Z_{j,i}=+1$ or $Z_{j,i}=-1$, where $i$ labels measurements. In the limit $N\gg 1$, the sample mean $\overline{Z}_{j}$ of the projective measurements for detuning $D_{j}$ becomes Gaussian distributed,

[TABLE]

where $\mathrm{var}(\sigma_{z})_{j}\equiv\langle\sigma_{z}(t_{f})^{2}\rangle_{j}-\langle\sigma_{z}(t_{f})\rangle^{2}_{j}$ is the expected variance of $\sigma_{z}$ averaged over realizations of the noise process, for detuning $D_{j}$. Assuming that $\mathrm{var}(\overline{Z}_{j})=\mathrm{var}(\sigma_{z})_{j}/N$ is the same for all relevant detunings, $\mathrm{var}(\overline{Z})_{j}\equiv\mathrm{var}(\overline{Z})\,\forall\,j$, Eq. (21) then corresponds to the conditional normal model of linear regression [S4], in which deviations of the data points from the expected linear behavior are given by independent and identically distributed (i.i.d.) Gaussian random variables. This assumption of uniform variance is justified in an approximate sense for ideal projective measurements. Indeed, in this situation, $\mathrm{var}(\overline{Z}_{j})=[1-\langle\sigma_{z}(t_{f})\rangle_{j}^{2}]/N$, so that $\mathrm{var}(\overline{Z}_{j})$ is independent of $D_{j}$ to first order in $\langle\sigma_{z}(t_{f})\rangle_{j}\approx(D+\mu_{B})F_{0,y}(0,t_{f})$, with $\mathrm{var}(\overline{Z}_{j})\approx 1/N$. Limiting ourselves to detunings for which $\langle\sigma_{z}(t_{f})\rangle_{j}\lesssim 0.05$ (see Fig. 4a), we find that $\mathrm{var}(\overline{Z}_{j})$ (estimated from the sample mean of measurements of $\sigma_{z}$) varies by less than $5\%$ across values of $D_{j}$.

To define our estimator of $\mu_{B}$ within the conditional normal model of linear regression, we first introduce the quantities $a\equiv\mu_{B}F_{0,y}(0,t_{f})$ and $b\equiv F_{0,y}(0,t_{f})$, corresponding to the $y$-intercept and slope of the linear equation $\langle\sigma_{z}(t_{f})\rangle_{j}=(D_{j}+\mu_{B})F_{0,y}(0,t_{f})\equiv a+b\,D_{j}$, respectively. Maximizing the likelihood of $a$ and $b$ with respect to measurement outcomes $\overline{Z}_{j}$ with the probability distribution given by Eq. (21) then yields the estimators

[TABLE]

The estimators defined by Eqs. (22) and (23) are Gaussian random variables with $E(a^{\mathrm{est}})=a$, $E(b^{\mathrm{est}})=b$, and

[TABLE]

To estimate $\mu_{B}$, we use

[TABLE]

When $\mathrm{var}(a^{\mathrm{est}})^{1/2}$ and $\mathrm{var}(b^{\mathrm{est}})^{1/2}$ are sufficiently small, we expand $\tilde{\mu}_{B}^{\mathrm{est}}$ in powers of $\delta a^{\mathrm{est}}\equiv a^{\mathrm{est}}-a$ and $\delta b^{\mathrm{est}}\equiv b^{\mathrm{est}}-b$. Truncating to the first order in $\delta a^{\mathrm{est}}$ and $\delta b^{\mathrm{est}}$, $\tilde{\mu}_{B}^{\mathrm{est}}$ becomes approximately Gaussian-distributed with ${\mathbb{E}}(\tilde{\mu}_{B}^{\mathrm{est}})\approx\mu_{B}$ and

[TABLE]

For the experimental data presented in Fig. 4a of the main text, $\mathrm{var}(\tilde{\mu}_{B}^{\mathrm{est}})$ is estimated by replacing $a\rightarrow a^{\mathrm{est}}$ and $b\rightarrow b^{\mathrm{est}}$ in Eq. (25), with $a^{\mathrm{est}}$ and $b^{\mathrm{est}}$ given by Eqs. (22) and (23), respectively.

Finally, to isolate the shift in the qubit frequency due to the first cumulant of the engineered source of noise, we apply the above procedure first in the absence of noise, and then in its presence. This yields the two sets of data points shown in Fig. 4a of the main text. Fitting a straight line through each data set and substracting their $x$-intercept then gives our final estimate of $\mu_{B}$, $\mu_{B}^{\mathrm{est}}=\mu_{B}^{\mathrm{on}}-\mu_{B}^{\mathrm{off}},$ where $\mu_{B}^{\mathrm{on(off)}}$ is the estimator defined by Eq. (24) with ($\mu_{B}^{\mathrm{on}}$) or without ($\mu_{B}^{\mathrm{off}}$) engineered noise. The variance of $\mu_{B}^{\mathrm{est}}$ is then simply $\mathrm{var}(\mu_{B}^{\mathrm{est}})=\mathrm{var}(\mu_{B}^{\mathrm{on}})+\mathrm{var}(\mu_{B}^{\mathrm{off}})$. For the experimental data presented in the main text, we find $\mu_{B}^{\mathrm{est}}/2\pi=127.1$ kHz with a standard deviation $\mathrm{var}(\mu_{B}^{\mathrm{est}})^{1/2}=3.86$ kHz, corresponding to the 95% confidence interval $\mu_{B}^{\mathrm{est}}/2\pi=127.1\,\pm\,7.56$ kHz.

2.2 PSD estimation procedure

To estimate the PSD, we build on the frequency-comb approach introduced by Alvarez and Suter in Ref. [13]. As detailed in the main text, treating the FFs as frequency combs generates a system of linear equations, solving which determines the PSD sampled at the harmonic frequencies. Rather than solving this sytem by matrix inversion as in Ref. [13], we employ a statistically-motivated maximum likelihood estimate (MLE), which takes experimental error into account. The likelihood function we use follows from the asymptotic Gaussian distribution of the measured decay constants, which we describe in Sec. 2.2.1. In Sec. 2.2.2, we determine the likelihood, the conditional probability of obtaining a set of decay data conditioned on the actual value of the PSD. The task of maximizing the likelihood can be cast as a linear problem, allowing clear comparison with Ref. [13].

2.2.1 Distribution of the decay constant

The PSD enters the dynamics of the qubit through the decay constant in Eq. (1), which can be obtained from measurements of the transverse Pauli operators, $\sigma_{x}$ and $\sigma_{y}$. Let $\sigma_{i}^{\text{est}}$ denote the estimated expected value of $\sigma_{i}$ for $i\in\{x,y\}$ after the qubit has evolved for a time $t$ under control sequence $p$. In the limit of a large number of measurements, $\sigma_{i}^{\text{est}}$ is approximately Gaussian distributed with mean $\mu_{i}=\mathbb{E}[\langle\sigma_{i}(t)\rangle]$ and variance $\text{var}[\sigma_{i}^{\text{est}}]$. In terms of the estimated expected values, the estimated decay constant is

[TABLE]

where $\tilde{\sigma}_{i}^{\text{est}}=\sigma_{i}^{\text{est}}-\mu_{i}$. Note that ${\sigma}_{i}^{\text{est}}$, $\mu_{i}$ and $\tilde{\sigma}_{i}^{\text{est}}$ on the right side of this expression depend implicitly on the time $t$. When $\text{var}[\sigma_{x}^{\text{est}}],\,\text{var}[\sigma_{y}^{\text{est}}]\ll 1$, we can expand this expression about $\tilde{\sigma}_{x}^{\text{est}},\tilde{\sigma}_{y}^{\text{est}}\approx 0$, yielding

[TABLE]

Since it is a linear combination of Gaussian distributed random variables, the decay constant is also Gaussian distributed with mean and variance,

[TABLE]

2.2.2 Maximum likelihood estimate

Recall from the main text that after $M\gg 1$ repetitions of a control sequence $p$ with cycle time $T$, the FF in Eq. (27) is an approximate frequency comb, enabling us to write the decay constant as

[TABLE]

Using the symmetry $S(\omega)=S(-\omega)$, and the decay of the PSD and FF at high frequencies, we can truncate the sum above to a finite number of harmonic frequencies,

[TABLE]

From the expected value above and the variance in Eq. (28), the conditional probability of measuring $\chi_{p}^{\text{est}}(MT)$ given the actual values of the PSD, $\vec{S}=[S(0),\ldots,S(K\omega_{h})]^{T}$, is

[TABLE]

The estimate of the PSD is based on the likelihood or conditional probability of measuring $\vec{\chi}=[\chi_{1}^{\text{est}}(MT),\ldots,\chi_{P}^{\text{est}}(MT)]^{T}$ for a set of control sequences $p=1,\ldots,P$ with $P\geq K$. Since the measurements of the decay constant are uncorrelated, the likelihood follows from Eq. (29),

[TABLE]

where the $P\times P$ covariance matrix has elements $\mathbf{\Sigma}_{p,q}=\text{var}[\chi_{p}^{\text{est}}(MT)]\,\delta_{p,q}$ and the $P\times K$ reconstruction matrix $\mathbf{B}$ depends on the FFs evaluated at the harmonic frequencies

[TABLE]

Since the likelihood is Gaussian, the maximum likelihood estimate of the PSD is equivalent to the value of $\vec{S}$ that minimizes the argument of the exponential in Eq. (30),

[TABLE]

The least squares estimate of the PSD originally used in Ref. [13], given by $\vec{S}^{\text{\,LS}}=\mathbf{B}^{-1}\vec{\chi}$, is recovered when $\mathbf{\Sigma}=\mathbf{I}$. This implies that all measurements of the decay constant contribute equally to the estimate. In contrast, the dependence of $\vec{S}^{\,\text{MLE}}$ on the actual variances in $\mathbf{\Sigma}$ ensures that measurements with more uncertainty contribute less to the estimate.

In the experiment, the PSD is reconstructed at the harmonics $k=0,\ldots,7$ using the $P=11$ control sequences depicted in Fig. (S3), all with cycle time $T=960\,$ns. This differs from the implementation of Ref. [13], which uses CPMG control sequences of varying cycle times. For each control sequence, the transverse Pauli components are measured to obtain $\sigma_{x}^{\text{est}}$, $\sigma_{y}^{\text{est}}$ and $\chi_{p}^{\text{est}}(MT)$ from Eq. (26). The variances of $\chi_{p}^{\text{est}}(MT)$, which comprise $\mathbf{\Sigma}$, follow from Eq. (28) with $\mu_{x}$ and $\mu_{y}$ replaced by the estimated values $\sigma_{x}^{\text{est}}$ and $\sigma_{y}^{\text{est}}$. In principle, the use of $P=11$ control sequences would enable us to reconstruct the PSD at $K=11$ harmonics. For the particular set of control sequences we used, however, the reconstruction matrix $\mathbf{B}$ becomes ill-conditioned for $K>8$, limiting the number of reconstructable harmonics.

2.3 Bispectrum estimation procedure

While this work is based on the non-Gaussian QNS protocol originally proposed in Ref. [41], the estimation procedure we implemented contains several innovations aimed at generalizing the noise model and improving robustness in the presence of experimental error and numerical instability. First, the zero-mean, non-Gaussian noise model of Ref.  [41] is insufficient to describe the square noise engineered in our quarton-qubit sensor, which is inherently nonzero mean. This complicates the estimation procedure, since both the bispectrum and the mean enter the qubit dynamics through the phase in Eq. (2). Estimating the bispectrum requires that we disambiguate the phase contribution of the bispectrum from that of the mean, which we accomplish by first estimating the noise mean and then isolating the dynamical contribution of the bispectrum in the non-Gaussian phase. A second key difference is the “single-shot” nature of the current estimation procedure. In Ref. [41], control sequences with non-zero filter order were first used to estimate the bispectrum at “non-zero harmonics”, i.e., $(k_{1}\omega_{h},k_{2}\omega_{h})$ for which $k_{1},k_{2}\neq 0$. This estimate was combined with subsequent phase measurements to estimate the bispectrum at “zero harmonics”, i.e. $(k_{1}\omega_{h},k_{2}\omega_{h})$ for which $k_{1}=0$ and/or $k_{2}=0$. In the RMLE estimate of the spectrum presented here, both the zeros and the non-zero harmonics are estimated simultaneously, eliminating any compounding of error that can occur in the two-step procedure. The present work also departs from Ref. [41] significantly in its use of a statistically motivated maximum likelihood estimation procedure. As discussed in the main text, the least-squares estimate of Ref. [41] is susceptible to numerical instability and, additionally, does not take measurement error into account.

In the remainder of this section, we fully detail our bispectrum estimation procedure. We begin by describing the probability distribution of the estimated “non-Gaussian” phase angle ($\varphi_{p}$), which enables us to derive the likelihood function for the probability of obtaining a particular set of phase data conditioned on the actual value of the bispectrum (Subsection 2.3.3). From the likelihood, the task of reconstructing the bispectrum can be mapped into an RMLE problem, as shown Subsection 2.3.4. The RMLE approach increases numerical stability, accounts for experimental error and allows us to deploy prior knowledge of the bispectrum in the estimation procedure. Since regularization can introduce error into the estimate if it is too strong, in Subsection 2.3.5 we determine an appropriate regularization strength for our problem using the L-curve criterion.

2.3.1 Distribution of the non-Gaussian phase

The non-Gaussian phase of the qubit is determined from the estimated expected values of the transverse Pauli operators, $\sigma_{x}^{\text{est}}$ and $\sigma_{y}^{\text{est}}$, when the qubit has evolved under control sequence $p$ for a time $t$. Recall from Sec. 2.2 that in the limit of a large number of measurements, $\sigma_{i}^{\text{est}}$ is approximately Gaussian distributed with mean $\mu_{i}=\mathbb{E}[\langle\sigma_{i}(t)\rangle]$ and variance $\text{var}[\sigma_{i}^{\text{est}}]$. From the estimated expected values, the ordinary phase is determined by

[TABLE]

where ${\sigma}_{i}^{\text{est}}$, ${\mu_{i}}$ and $\tilde{\sigma}_{i}^{\text{est}}=\sigma_{i}^{\text{est}}-{\mu_{i}}$ depend implicitly on the time $t$. When $\text{var}[\sigma_{x}^{\text{est}}]$, $\text{var}[\sigma_{y}^{\text{est}}]\ll 1$, we can expand $\phi^{\text{est}}_{p}$ about $\tilde{\sigma}_{x}^{\text{est}}$, $\tilde{\sigma}_{y}^{\text{est}}$ $\approx 0$, yielding

[TABLE]

As a linear combination of Gaussian distributed random variables, the phase is also Gaussian distributed with mean and variance

[TABLE]

The second equality on the right-hand side of $\mathbb{E}[\phi_{p}^{\text{est}}]$ follows from Eq. (2). Subtracting out the contribution of the noise mean from the phase produces the non-Gaussian phase,

[TABLE]

where $\mu_{B}^{\text{est}}$ is the estimated noise mean described in Sec. 2.1. Using the asymptotic Gaussian distribution of $\mu_{B}^{\text{est}}$ with mean $\mu_{B}$ and variance $\text{var}[\mu_{B}^{\text{est}}]$, the non-Gaussian phase is similarly Gaussian with mean and variance

[TABLE]

Note that $\mathbb{E}[\varphi_{p}^{\text{est}}(t)]$ depends to leading order on the bispectrum, unlike $\mathbb{E}[\phi_{p}^{\text{est}}(t)]$ above.

2.3.2 Restriction to the principal domain

For any real, classical process, the bispectrum has three general symmetries: (1) $S_{2}(\omega_{1},\omega_{2})=S_{2}(\omega_{2},\omega_{1})$ (permutation symmetry); (2) $S_{2}(\omega_{1},\omega_{2})=S_{2}(-\omega_{1},-\omega_{2})$ (invariance under complex conjugation); (3) $S_{2}(\omega_{1},\omega_{2})=S_{2}(-\omega_{1}-\omega_{2},\omega_{2})$ (stationarity). These symmetries define the 12 regions of the frequency plane depicted in Fig. 3c. If $(\omega_{1},\omega_{2})\in\text{int}(\mathcal{D}_{2})$ is contained in the interior of the principal domain, the symmetries imply

[TABLE]

In other words, the bispectrum takes a value equivalent to $S_{2}(\omega_{1},\omega_{2})$ in each of the 12 regions. This is summarized by the multiplicity, $m(\omega_{1},\omega_{2})=12$. For $(\omega,\omega)$, which lies on the boundary of $\mathcal{D}_{2}$,

[TABLE]

implying $m(\omega,\omega)=6$. For $(\omega,0)$, which also lies on the boundary of $\mathcal{D}_{2}$,

[TABLE]

similarly implying $m(\omega,0)=6$. Note that $(\omega_{1},\omega_{2})=(0,0)$ is invariant under all of the symmetries, implying that $m(0,0)=1$.

The symmetries and multiplicities simplify the expected value of the phase substantially. Recall that after $M\gg 1$ repetitions of control sequence $p$ with cycle time $T$, the frequency comb approximation enables us to write the expected phase as a discrete sum depending on the bispectrum and the FF evaluated at the harmonic frequencies,

[TABLE]

In terms of the multiplicities, we can rewrite the sum as

[TABLE]

Using $S_{2}(\omega_{h}\vec{k})=S_{2}(-\omega_{h}\vec{k})$ and $G_{p}(\omega_{h}\vec{k},T)^{*}=G_{p}(-\omega_{h}\vec{k},T)$, and truncating the sum to a finite subset $\mathcal{K}_{2}$, we obtain

[TABLE]

2.3.3 Likelihood function $P(\vec{\varphi}|\vec{S}_{2})$

Given the actual bispectrum in $\mathcal{K}$, $\vec{S}_{2}=[S_{2}(\omega_{h}\vec{k}_{1}),\ldots,S_{2}(\omega_{h}\vec{k}_{N})]^{T}$, the conditional probability of measuring $\varphi_{p}^{\text{est}}(MT)$ follows from Eqs. (35)-(37),

[TABLE]

Reconstructing the bispectrum requires measurements the non-Gaussian phase for a set of control sequences $p=1,\ldots,P$ with $P\geq N$, which we gather into the column vector $\vec{\varphi}=[\varphi_{1}^{\text{est}}(MT),\ldots,\varphi_{P}^{\text{est}}(MT)]^{T}$. Because the non-Gaussian phase measurements are uncorrelated, the likelihood or probability of obtaining $\vec{\varphi}$ given $\vec{S}_{2}$ is a product of the conditional probabilities for the complete set of control sequences,

[TABLE]

Here, the $P\times P$ covariance matrix $\mathbf{\Sigma}$ is diagonal with elements $\mathbf{\Sigma}_{p,q}=\text{var}[\varphi_{p}^{\text{est}}(MT)]\,\delta_{p,q}$, and the $P\times N$ reconstruction matrix $\mathbf{A}$ depends on the filter functions evaluated at the harmonic frequencies,

[TABLE]

In the experiment, the likelihood in Eq. (39) is determined by measuring the non-Gaussian phase for each of the $P=11$ control sequences in Fig. (S3). Since we also rely on these sequences to estimate the PSD, both the bispectrum and the PSD can be estimated with the the same set of transverse Pauli measurements. For each of the control sequences, $\varphi_{p}^{\text{est}}(MT)$ was determined from Eq. (32) using measurements of $\sigma_{x}^{\text{est}}$, $\sigma_{y}^{\text{est}}$ and $\mu_{B}^{\text{est}}$. The variances of the $\varphi_{p}^{\text{est}}(MT)$, which constitute the covariance matrix $\mathbf{\Sigma}$, are given by Eq. (36) with $\mu_{x}$ and $\mu_{y}$ replaced by the estimated values $\sigma_{x}^{\text{est}}$, $\sigma_{y}^{\text{est}}$. The reconstruction matrix $\mathbf{A}$ is determined from Eq. (40), with each FF evaluated on the set of harmonics $\mathcal{K}_{1}$ depicted in Fig. 3(b).

2.3.4 Regularized maximum likelihood estimation

For the Gaussian likelihood derived in the previous section, the maximum likelihood estimate (MLE) of the bispectrum is equivalent to the value of $\vec{S}_{2}$ that minimizes the exponent in Eq. (39),

[TABLE]

In the special case where $\mathbf{\Sigma}\propto\mathbf{I}$, we recover the least-squares estimate used in Ref. [41] with solution $\vec{S}_{2}^{\text{\,LS}}=\mathbf{A}^{-1}\vec{\varphi}$. Even with a nonuniform covariance matrix, Eq. (41) is a simple convex optimization problem admitting an analytic solution for $\vec{S}_{2}^{\,\text{MLE}}$. When $\mathbf{A}$ is ill-conditioned, however, the MLE suffers from numerical instability, which can introduce significant error into the estimate of the bispectrum, despite the existence of an analytic solution. The problem can be made more stable by introducing a regularization term or “regularizer” $Q(\vec{S}_{2})$, producing the regularized maximum likelihood estimate (RMLE) of the bispectrum,

[TABLE]

The regularizer imposes additional structure on the solution, making it more robust to numerical instability arising from $\mathbf{A}$. It also prevents overfitting, in which the estimated bispectrum is unduly influenced by errors in $\vec{\varphi}$ and is, thus, a poor predictor of the qubit dynamics under more general control settings.

There are numerous methods of regularization for ill-conditioned and/or ill-posed problems. An approach particularly ammenable to maximum likelihood estimation is Tikhonov regularization, which employs an L2 regularizer $Q(\vec{S}_{2})=|\!|\lambda\vec{S}_{2}|\!|_{2}^{2}$ with strength controlled by the regularization parameter $\lambda\geq 0$ [S5]. In Eq. (42), this regularizer has the effect of penalizing $\vec{S}_{2}$ with larger $L_{2}$-norms. To estimate the bispectrum, we consider a variation of Tikhonov regularization in which

[TABLE]

where $\mathbf{D}=\mathrm{diag}~{}(d_{1},\cdots,d_{N})$ is the diagonal “smoothing matrix”. Note that the Tikhonov regularizer is recovered when $d_{1}=\ldots=d_{N}=1$ (Fig. S5a) and the standard maximum likelihood estimate is recovered when $\lambda=0$. Using non-uniform values for the diagonals enables us to incorporate prior information about the bispectrum. For example, if the magnitude of the bispectrum is known to decay at the high-frequency border of $\mathcal{K}$, we can make the corresponding harmonics in $\mathbf{D}$ large compared to those of the interior (Fig. S5b). Such a smoothing matrix favors a solution with small magnitude at the border. The connection between the smoothing matrix and prior knowledge of the bispectrum is more explicit in a Bayesian formulation of the estimation problem in which the RMLE estimate in Eq. (42) with the regularizer in Eq. (43) is equivalent to a posterior mean estimate in which the prior distribution of the bispectrum is Gaussian and zero-mean with covariance matrix $(2\lambda^{2}\mathbf{D}^{2})^{-1}$, provided $\mathbf{D}$ is full-rank. For any smoothing matrix, the regularized maximum likelihood estimate has the simple analytic solution,

[TABLE]

2.3.5 The L-curve criterion

Although it guards against numerical instability and overfitting to errors in the measured data, the regularizer can introduce its own error into the estimate if $\lambda$ too large. A fundamental challenge in regularization is selecting a value of $\lambda$ that balances these sources of error. While this problem is still an active area of research, one of the most widely used strategies for selecting $\lambda$ is the L-curve criterion [S6]. A graphical technique, the L-curve criterion enables one to visualize the magnitude of the regularization error in proportion to other errors in the estimate and choose $\lambda$ accordingly.

For a given $\lambda$, the sources of error that contribute to the regularized maximum likelihood estimate in Eq. (44) are described by the residual norm

[TABLE]

and the solution norm

[TABLE]

The regularization parameter $\lambda$ enters both the residual norm and solution norm implicitly through $\vec{S}_{2}^{\,\text{RMLE}}$. The residual norm increases as $\lambda$ grows, while the solution norm decreases. When $E(\lambda)$ is too large relative to $R(\lambda)$, the estimate does not account for the measured data due to error introduced by the regularization. Conversely, when $E(\lambda)$ is too small relative to $R(\lambda)$, the estimate is fits the measured data too closely, making it susceptible to overfitting and numerical error. The influence of $\lambda$ on the error conributions is captured by the L-curve, a parametric plot of $\text{log}R(\lambda)$ vs. $\text{log}E(\lambda)$ as a function of $\lambda$. A typical L-curve, with its characteristic “L” shape, is shown in Fig. S6. Note that as $\lambda$ increases from left to right, $\text{log}R(\lambda)$ sharply decreases and then plateaus, while $\text{log}E(\lambda)$ is initially stable followed by a rapid increase. The corner of the L-curve, marks a point at which the solution norm and residual norm are small simultaneously. The corner, thus, signifies the optimal value of $\lambda$ according to the L-curve criterion.

Figure  S7 shows L-curves generated by our experimental data for the two different smoothing matrices illustrated in Fig.  S5. For both the uniform and non-uniform smoothing matrices, the L-curves lack corners. Unlike the typical L-curve in Fig. S6, $\text{log}R(\lambda)$ does not exhibit a sharp increase as $\lambda\rightarrow 0$. This indicates that, for the control sequences we have selected, the reconstruction matrix $\mathbf{A}$ is sufficiently well conditioned to make regularization the dominant source of error [S7]. Consequently, it is not optimal to utilize regularization in this setting and the reconstruction presented in the main text uses $\lambda=0$.

Note that this finding is contingent on both $\mathbf{A}$ and the particular regularizer we employ. Estimating the bispectrum at a greater number of harmonics demands a larger $\mathbf{A}$, which is more likely to be near singular and/or poorly conditioned. This scenario will likely require some form of regularization. Additionally, the error introduced by regularization is reduced when prior knowledge of the bispectrum (if available) is used to select the regularizer. For example, suppose a noise model or previous experiment indicates that $\vec{S}_{2}$ takes a value in the vicinity of $\vec{S}_{\mu}$. This information is captured by the regularizer

[TABLE]

This corresponds to a Gaussian prior distribution of $\vec{S}_{2}$ with mean $\vec{S}_{\mu}$ and covariance matrix $(2\lambda^{2}\mathbf{D}^{2})^{-1}$. In contrast, naively employing Tikhonov regularization amounts to assuming a zero-mean prior distribution of $\vec{S}_{2}$.

Additional references

• S1.

V. E. Manucharyan, J. Koch, L. I. Glazman, and M. H. Devoret, “Fluxonium: Single Cooper-pair circuit free of charge offsets,” Science 326, 113–116 (2009).

• S2.

E. Magesan, J. M. Gambetta, and J. Emerson, “Scalable and robust randomized benchmarking of quantum processes,” Phys. Rev. Lett. 106, 180504 (2011).

• S3.

C. Macklin, K. O’Brien, D. Hover, M. E. Schwartz, V. Bolkhovsky, X. Zhang, W. D. Oliver, and I. Siddiqi, “A near-quantum-limited Josephson traveling-wave parametric amplifier,” Science 350, 307–310 (2015).

• S4.

G. Casella and R. L. Berger, Statistical Inference (Duxbury Pacific Grove, CA, 2002).

• S5.

A. Tikhonov, “Solution of incorrectly formulated problems and the regularization method,” Dokl. Akad. Nauk 151, 1035–1038 (1963).

• S6.

P. C. Hansen, “The L-curve and its use in the numerical treatment of inverse problems,” in Computational Inverse Problems in Electrocardiology (WIT Press, 2000) pp. 119–142.

• S7.

T. Regińska, “A regularization parameter in discrete ill- posed problems,” SIAM J. Sci. Comput. 17, 740–749 (1996).