=0cm

# Probing flux and charge noise with macroscopic resonant tunneling

Alexander M. Whiticar  D-Wave Systems Inc., 3033 Beta Avenue, Burnaby BC
Canada V5G 4M9  Anatoly Y. Smirnov  D-Wave Systems Inc., 3033 Beta Avenue,
Burnaby BC Canada V5G 4M9  Trevor Lanting  D-Wave Systems Inc., 3033 Beta
Avenue, Burnaby BC Canada V5G 4M9  Jed Whittaker  D-Wave Systems Inc., 3033
Beta Avenue, Burnaby BC Canada V5G 4M9  Fabio Altomare  D-Wave Systems Inc.,
3033 Beta Avenue, Burnaby BC Canada V5G 4M9  Teresa Medina  D-Wave Systems
Inc., 3033 Beta Avenue, Burnaby BC Canada V5G 4M9  Rahul Deshpande  D-Wave
Systems Inc., 3033 Beta Avenue, Burnaby BC Canada V5G 4M9  Sara Ejtemaee
D-Wave Systems Inc., 3033 Beta Avenue, Burnaby BC Canada V5G 4M9  Emile
Hoskinson  D-Wave Systems Inc., 3033 Beta Avenue, Burnaby BC Canada V5G 4M9
Michael Babcock  D-Wave Systems Inc., 3033 Beta Avenue, Burnaby BC Canada V5G
4M9  Mohammad H. Amin  D-Wave Systems Inc., 3033 Beta Avenue, Burnaby BC
Canada V5G 4M9  Department of Physics, Simon Fraser University, Burnaby, BC
Canada V5A 1S6

###### Abstract

We report on measurements of flux and charge noise in an rf-SQUID flux qubit
using macroscopic resonant tunneling (MRT). We measure rates of incoherent
tunneling from the lowest energy state in the initial well to the ground and
first excited states in the target well. The result of the measurement
consists of two peaks. The first peak corresponds to tunneling to the ground
state of the target well, and is dominated by flux noise. The second peak is
due to tunneling to the excited state and is wider due to an intrawell
relaxation process dominated by charge noise. We develop a theoretical model
that allows us to extract information about flux and charge noise within one
experimental setup. The model agrees very well with experimental data over a
wide dynamic range and provides parameters that characterize charge and flux
noise.

##  I Introduction

Improving the performance of superconducting quantum computing technologies
relies on reducing the impact of noise sources that lead to decoherence de
Leon et al. (2021). This can be achieved by designing noise-resistant circuits
and developing lower-loss materials Yan et al. (2016a); Nguyen et al. (2019);
Place et al. (2021); Siddiqi (2021). The dominant noise sources affecting
superconducting qubits are flux and charge noise, which are thought to
originate from ensembles of microscopic systems manifesting as materials
defects. For qubits implemented with a superconducting quantum interference
device (SQUID), a ubiquitous $1/f$ flux noise spectrum has been observed
Bialczak et al. (2007); Lanting et al. (2009); Quintana et al. (2017);
Braumüller et al. (2020). Although a concrete microscopic mechanism for flux
noise has yet to be determined, the prevailing models suggest that randomly
oriented electronic spins at the metal-oxide interface lead to inductive
losses Koch et al. (2007); Anton et al. (2013); Lanting et al. (2014, 2020).
Similarly, defects in dielectrics are thought to cause dielectric losses by
coupling to and extracting energy from the qubit’s electric field Martinis et
al. (2005); Müller et al. (2019). To design next generation hardware, it is
crucial to be able to distinguish and quantify the strength of each noise
source in current hardware Braumüller et al. (2020).

Macroscopic resonant tunneling (MRT) uses flux-tunable, multi-well qubits to
measure the noise affecting the flux qubits Harris et al. (2008). An MRT
experiment consists of measuring the incoherent tunneling rate between the
flux states of the left and right wells of the qubit potential as a function
of flux bias $\Phi^{x}$ (see Fig. 1). When the energy levels are aligned, the
observed MRT peak is shaped by details of the noise spectral density. Previous
work on MRT primarily focused on the details of the peak originating from the
tunneling between the lowest energy levels in the initial and target wells,
which is dominated by flux noise Harris et al. (2008); Lanting et al. (2011).
In this article, we report measurements of the lowest energy transition
(zeroth peak) and the first excited transition (first peak) corresponding to
incoherent tunneling to the first excited state within the target well.
Intrawell relaxation from the first excited state to the ground state inside
the target well leads to an additional broadening of the first peak. This
intrawell relaxation is dominated by charge noise and allows us to
characterize the strength of coupling to charge fluctuations. To extract
information on the noise affecting the qubits, we develop a theoretical model
that combines interwell and intrawell relaxation and takes into account both
charge and flux noise that can be fit to experimental data.

##  II Theoretical model

Figure 1:  (a) Schematic diagram of an rf-SQUID flux qubit with two external
fluxes, $\Phi^{x}$ and $\Phi_{\rm CJJ}^{x}$, threading the two loops. (b)
Classical potential energy of the rf-SQUID flux qubit versus the flux $\Phi$
induced in the main loop. Four metastable energy levels are shown, two in each
well.

We consider a compound Josephson junction (CJJ) rf-SQUID flux qubit,
schematically represented in Fig. 1a Harris et al. (2010). The qubit consists
of two loops, main and CJJ, threaded by two external fluxes, $\Phi^{x}$ and
$\Phi^{x}_{\rm CJJ}$. The potential energy of the qubit has double-well shape
as a function of flux $\Phi$ threading the main loop (Fig. 1b). When the
barrier is high, tunneling amplitude between the two wells is small. This
allows us to define metastable states $\ket{n}$ with energies $E_{n}$ in each
well, and introduce tunneling amplitudes $\Delta_{mn}$ between states
$\ket{m}$ and $\ket{n}$ in opposite wells, as described in Appendix A. The
Hamiltonian of the system in this basis is written as

$H_{\rm S}=\sum_{n}E_{n}\ket{n}\bra{n}-{1\over 2}\sum_{m\neq
n}\Delta_{mn}\ket{m}\bra{n}.$ (1)

We numerate states in the left (right) well with even (odd) integers (see Fig.
1b). For simplicity, we assume in the following that we tunnel from the left
initial state $\ket{0}$ into the right well.

The rf-SQUID is dominantly coupled to flux and charge noise via current $I$
through the main loop and voltage $V$ across the junctions, with an
interaction Hamiltonian

$\displaystyle H_{\rm int}=-\sum_{m,n}(I_{mn}\delta\Phi+V_{mn}\delta
q)\ket{m}\bra{n},$ (2)

where $I_{mn}=\bra{m}I\ket{n}$ and $V_{mn}=\bra{m}V\ket{n}$. The flux noise,
$\delta\Phi$, and charge noise, $\delta q$, are characterized by noise
spectral densities $S_{\Phi}(\omega)$ and $S_{q}(\omega)$, respectively. Flux
noise is taken to be a sum of low-frequency and high-frequency components:
$S_{\Phi}(\omega)=S^{\Phi}_{L}(\omega)+S^{\Phi}_{H}(\omega)$. The low-
frequency part is characterized by its r.m.s. value $W_{\Phi}$, and the high
frequency component is assumed to be ohmic parameterized by a dimensionless
parameter $\eta$. Also, charge noise is described by dielectric loss tangent
$\tan\delta_{C}$. Details of the spectral densities and the noise parameters
are provided in Appendix B.

At time $t=0$, the system is initialized in the lowest energy state of the
left well, $\ket{0}$, with probability $P_{0}=1$. The rate of transition out
of this initial state is given by

$\Gamma(\epsilon)=-\left[{dP_{0}\over
dt}\right]_{t=0}=\sum_{n}\Gamma_{0n}(\omega_{01}-\omega_{n1}),$ (3)

where $\omega_{n1}=(E_{n}-E_{1})/\hbar$, and
$\epsilon=E_{0}-E_{1}=\hbar\omega_{01}$ is the energy bias from the degeneracy
point. The functions $\Gamma_{0n}(\omega)$ describe resonant tunneling from
$\ket{0}$ to the state $\ket{n}$ in the target well.

While flux noise directly affects the transition between the wells, charge
noise broadens the transition peak indirectly via intrawell relaxation. When
the state $\ket{n}$ is an excited state in the target well, after tunneling,
the system will quickly relax down to the lowest energy state, $\ket{1}$. The
energy uncertainty due to this intrawell relaxation leads to an additional
broadening of the transition peak, with a width proportional to the rate of
relaxation $\Gamma_{n1}(\omega)$ from $\ket{n}$ to $\ket{1}$. The total
transition rate is described by a convolution of three functions, each
corresponding to one component of noise (see Appendix C):

$\Gamma_{0n}(\omega)={\Delta_{0n}^{2}\over 4\hbar^{2}}\,G_{0n}(\omega),$ (4)

where

$\displaystyle G_{0n}(\omega)=\\!\\!\\!\int\\!{d\omega^{\prime}\over
2\pi}{d\omega^{\prime\prime}\over
2\pi}G_{L}(\omega{-}\omega^{\prime})G_{H}(\omega^{\prime}{-}\omega^{\prime\prime})G_{R}(\omega^{\prime\prime}).\quad$
(5)

Here, we have defined single-peaked functions

$\displaystyle G_{L}(\omega)$ $\displaystyle=$
$\displaystyle{\sqrt{2\pi}\hbar\over
W}\exp\left\\{{-(\hbar\omega-\epsilon_{p})^{2}\over 2W^{2}}\right\\},$ (6)
$\displaystyle G_{H}(\omega)$ $\displaystyle=$
$\displaystyle{2\hbar\gamma\over\hbar^{2}\omega^{2}+\gamma^{2}}{\hbar\omega/k_{B}T\over
1-e^{-\hbar\omega/k_{B}T}},$ (7) $\displaystyle G_{R}(\omega)$
$\displaystyle=$
$\displaystyle{2\Gamma_{n1}(\omega+\omega_{n1})\over\omega^{2}+\Gamma_{n1}^{2}(\omega+\omega_{n1})}.$
(8)

The Gaussian function $G_{L}(\omega)$ represents broadening due to low-
frequency flux noise Amin and Averin (2008); Harris et al. (2008). The width
of the peak, $W$, is proportional to the r.m.s. value of low-frequency flux
noise. The shift, $\epsilon_{p}$, is related to $W$ by the fluctuation-
dissipation theorem:

$W^{2}=2k_{B}T\epsilon_{p}.$ (9)

High frequency flux noise is included via the Lorentzian-like function
$G_{H}(\omega)$, with broadening determined by $\gamma$ Amin and Averin
(2008); Smirnov and Amin (2018). Finally, intrawell relaxation is captured by
$G_{R}(\omega)$. The function $\Gamma_{n1}(\omega)$ represents intrawell
relaxation from state $\ket{n}$ to $\ket{1}$ within the target well. Naturally
for the lowest MRT peak with $n=1$, there is no intrawell relaxation.
Therefore, $\Gamma_{11}=0$ and $G_{R}(\omega)=2\pi\delta(\omega)$, turning Eq.
(4) into a single convolution integral:

$\displaystyle G_{01}(\omega)=\\!\\!\\!\int{d\omega^{\prime}\over
2\pi}G_{L}(\omega{-}\omega^{\prime})G_{H}(\omega^{\prime}).\quad$ (10)

For transition to the first excited state in the target well with $n=3$, we
approximately write (see Appendix C)

$\displaystyle\Gamma_{31}(\omega)$ $\displaystyle=$
$\displaystyle{\zeta\over\hbar}{\tanh(\hbar\omega/k_{B}T)\over
1-e^{-\hbar\omega/k_{B}T}},$ (11)

where $\zeta$ is a charge noise broadening coefficient in units of energy. For
multi-level MRT peaks with $n>3$, there are more than one intrawell relaxation
channels, and $\Gamma_{n1}$ is the sum of all of them. In this paper, however,
we only focus on the first two MRT peaks. This model is an extension to
previous models that relied on a convolution of two noise sources Amin and
Averin (2008); Smirnov and Amin (2018). A more formal derivation of the model
introduced here can be found in [21]. Note that all $G_{\kappa}(\omega)$
functions, with $\kappa=L,H,R,0n$, are approximately normalized (with slight
deviations due to non-ideal Lorentzian form) with normalization condition,

$\int{d\omega\over 2\pi}G_{\kappa}(\omega)=1.$ (12)

Any deviation from a perfect normalization will be absorbed into
$\Delta_{0n}$, when taken as a free parameter.

The parameters $\Delta_{0n}$, $W$, $\gamma$, and $\zeta$ are all in units of
energy and can be calculated using the underlying rf-SQUID Hamiltonian and
noise spectral densities as described in the appendices. The broadening
parameters $W$, $\gamma$, $\zeta$ would then depend on the target state
$\ket{n}$ and the energy bias $\epsilon$. Ignoring these small dependencies
and under some additional assumptions listed in Appendix D, we can treat them
as fitting parameters. This allows us to fit the model to the experimental
data and extract information about noise with no need for diagonalization of
the rf-SQUID Hamiltonian. The energy bias $\epsilon$ at which MRT peaks are
measured is obtained by applying the external flux $\Phi^{x}$. This bias can
be approximated by $\epsilon=2I_{P}\Phi^{x}$, where $I_{P}=(I_{11}-I_{00})/2$
is the persistent current, and $\Phi^{x}$ is measured from the degeneracy
point $\Phi_{0}/2$. One can therefore present the transition rate as a
function of the applied flux, $\Gamma(\Phi^{x})$, and express the noise
parameters in flux units

$W_{\Phi}={W\over 2I_{P}},\quad\epsilon_{\Phi}={\epsilon\over
2I_{P}},\quad\gamma_{\Phi}={\gamma\over 2I_{P}},\quad\zeta_{\Phi}={\zeta\over
2I_{P}}.$ (13)

While $W_{\Phi}$ directly measures the r.m.s. value of low-frequency flux
noise, the dimensionless ohmic coefficient of high-frequency flux noise, and
the loss tangent of charge noise are given by (see Appendix D)

$\eta={4I_{P}\gamma_{\Phi}\over
k_{B}T},\qquad\tan\delta_{C}={\zeta_{\Phi}\over\Phi^{x}_{31}},$ (14)

where $\Phi^{x}_{31}=\hbar\omega_{31}/2I_{P}$ is the distance between the two
peaks in flux units. It is also common to express ohmic flux noise in terms of
shunt resistance $R_{S}$ or inductive loss tangent $\tan\delta_{L}$ Nguyen et
al. (2019):

$R_{S}={8I_{P}^{2}L^{2}\over\hbar\eta}={2I_{P}L^{2}k_{B}T\over\hbar\gamma_{\Phi}},\qquad\tan\delta_{L}={\omega
L\over R_{S}},$ (15)

where $L$ is the inductance of the main loop.

Figure 2:  Macroscopic resonant tunneling rate as a function of external flux
$\Phi^{x}$ (measured from the degeneracy point $\Phi_{0}/2$) for different
values of (a) low-frequency flux noise broadening $W_{\Phi}$, (b) high-
frequency flux noise broadening $\gamma_{\Phi}$, and (c) charge noise
broadening $\zeta_{\Phi}$. Red and blue line-shapes in (a) and (b) correspond
to zeroth peak with left and right well initialization, respectively. The
line-shapes in (c) represent zeroth and first peaks corresponding to the left
initialization at a constant $\Phi^{x}_{\rm CJJ}$. We use $\Delta_{01}=2$ MHz
and $T=5$ mK in all simulations. We also use $\gamma_{\Phi}=10~{}\mu\Phi_{0}$
for (a), $W_{\Phi}=35~{}\mu\Phi_{0}$ for (b), and $\Delta_{03}=20$ MHz,
$W_{\Phi}=35~{}\mu\Phi_{0}$, $\gamma_{\Phi}=3~{}\mu\Phi_{0}$ for (c).

To illustrate the dependence of the MRT peaks on the broadening parameters, we
plot the transition rate $\Gamma(\Phi^{x})$ for different $W_{\Phi}$,
$\gamma_{\Phi}$, and $\zeta_{\Phi}$ in Fig. 2. Figure 2a,b highlights the
zeroth peaks for both left and right well initial state preparations (colored
as red and blue respectively). In Fig. 2a, the peaks are plotted for two
values of the width: $W_{\Phi}=10$ and $30\,\mu\Phi_{0}$. The shift of the
peak position from zero bias is $\epsilon_{\Phi}$, which is also changed
according to (9) assuming constant temperature Harris et al. (2008). Figure 2b
shows MRT peaks for different strengths of high frequency noise
$\gamma_{\Phi}$ while keeping all other parameters constant. For small values
$\gamma_{\Phi}$, low-frequency flux noise dominates, and a Gaussian broadened
line shape of width $W_{\Phi}$ is recovered Amin and Averin (2008); Harris et
al. (2008). When $\gamma_{\Phi}$ increases, an additional broadening develops
with a characteristic asymmetric tail extending to larger flux biases Lanting
et al. (2011). While intrawell relaxation does not contribute to the
broadening of $\Gamma_{01}$ in Eq. (3), it does for higher energy transitions
such as $\Gamma_{03}$ as shown in Fig. 2c. For increasing values of the
broadening parameter $\zeta_{\Phi}$, i.e., higher intrawell relaxation, the
width of the first peak and its contribution to the valley between the two
peaks are increased.

##  III Experimental Results

The MRT measurement protocol involves preparing the qubit in a known initial
state of the rf-SQUID double-well potential (ground state of the left or right
well) and measuring the tunneling rate into the adjacent well as a function of
flux bias as described in [16]. MRT measurements were performed on the
quantum processor of a D-Wave 2000Q™ lower noise system. The qubit has
external lines that apply fluxes $\Phi^{x}$ and $\Phi_{\rm CJJ}^{x}$ to the
main and CJJ loops, respectively. These lines enable time-dependent control
over the qubit potential energy, with $\Phi^{x}(t)$ setting the flux-bias tilt
between the left and right well, and $\Phi_{\rm CJJ}^{x}(t)$ tuning the
tunneling energy $\Delta$. Each qubit is controlled by external lines that
have a 3 and 30 MHz bandwidth, respectively, that enable in-situ MRT
measurements on individual qubits throughout the quantum processor.

We measured 27 parametrically identical qubits across the fabric of the
processor. The measurements were performed in a dilution refrigerator with a
base temperature of $\\!\sim$10 mK on a processor fully calibrated according
to the procedure described in [18]. The qubits had a critical current of
$I_{c}=2.30\pm 0.08~{}\mu{\rm A}$, an inductance of $L=250\pm 7~{}{\rm pH}$
and a capacitance of $C=110\pm 4~{}{\rm fF}$. A constant bias of $\Phi_{\rm
CJJ}^{x}=-0.74~{}\Phi_{0}$ was applied to facilitate measurements of the
tunneling rate that varied over four orders of magnitude as a function of
$\Phi^{x}$. At this value of $\Phi_{\rm CJJ}^{x}$, the qubits had a persistent
current of $I_{p}=1.37\pm 0.01~{}\mu{\rm A}$.

Figure 3 shows a typical dataset of a single qubit’s tunneling rate as a
function of flux bias with an initial state preparation in the left (right)
well shown in red (blue). At the degeneracy point, $\Phi^{x}=0$, the ground
states of the two wells are aligned. The data shows a resonant peak near this
point that corresponds to tunneling between these two states with an offset
depending on the state initialization (see Eq. 9 and [16]). This zeroth MRT
peak has a width dominated by low-frequency flux noise. Away from this peak,
the tunneling rate exhibits an asymmetric tail representative of high-
frequency flux noise, in qualitative agreement with Fig. 2b. Further
increasing $|\Phi^{x}|$ causes a gradual increase in the tunneling rate until
reaching the first peak at $|\Phi^{x}|\sim 2~{}\rm{m}\Phi_{0}$. At this point,
the initial state is aligned with the first excited state in the target well.
An additional broadening is observed on the first peak due to intrawell
relaxation in the target well. The line-shape near this peak and in the valley
between the two peaks provides information about the strength of the charge
noise.

Figure 3:  Measurement of macroscopic resonant tunneling rate $\Gamma$ as a
function of flux bias $\Phi^{x}$ for a single qubit at a fixed $\Phi_{\rm
CJJ}^{x}$ controlling the barrier height. The red (blue) colour shows the
initial state prepared in the left (right) well. The solid line is a fit to
the model described in Eqs. (16). The inset shows the zeroth MRT peaks, which
highlights the asymmetric tail resulting from high-frequency flux noise. The
first MRT peak has additional broadening due to intrawell relaxation, which
allows the extraction of the strength of charge noise.

We fit the model described in Eqs. (3)-(11) to the measured tunneling rates
with,

$\Gamma(\Phi^{x})={1\over
4}\big{(}\Delta_{01}^{2}G_{01}(\Phi^{x})+\Delta_{03}^{2}G_{03}(\Phi^{x})\big{)},$
(16)

where $G_{01}(\Phi^{x})$ and $G_{03}(\Phi^{x})$ are described by Eq. 10 and 5.

A typical best fit to the dataset is shown by the solid black line in Fig. 3
with the tunneling amplitudes $\Delta_{01}/h=2.72\pm 0.01~{}\rm{MHz}$ and
$\Delta_{03}/h=29.8\pm 0.2~{}\rm{MHz}$, and the noise broadening parameters
$W_{\Phi}=37.2\pm 0.1~{}\mu\Phi_{0}$, $\gamma_{\Phi}=0.54\pm
0.05~{}\mu\Phi_{0}$, and $\zeta_{\Phi}=4.53\pm 0.09~{}\mu\Phi_{0}$. Flux
offset $\Phi^{x}_{31}=2153.6\pm 0.5~{}\mu\Phi_{0}$, corresponding to
$\omega_{31}/2\pi=2I_{P}\Phi^{x}_{31}/h=9.17$ GHz, is used to fit the relative
position of the zeroth and first peaks. The fitting temperature of $T=7.3$ mK
matches the thermometry mounted on the mixing chamber plate. Using (14) and
(15) we estimate the noise parameters: $\eta=5.9\pm 0.5\times 10^{\rm{-2}}$,
$R_{\rm S}=147\pm 13~{}\rm{k}\Omega$, $\tan\delta_{C}=2.07\pm 0.04\times
10^{\rm{-3}}$, and $\tan\delta_{L}(1\,\text{GHz})=10.6\pm 0.9\times
10^{\rm{-6}}$.

Using the same measurement procedure we fit all 27 qubits to the hybrid noise
model and find consistent results. The data and fit for each of these qubits
is similar to Fig. 3. A summary of these results is presented in Fig. 4. We
find mean noise parameters of $\eta=5\pm 1\times 10^{\rm{-2}}$, $R_{\rm
S}=180\pm 40~{}\rm{k}\Omega$, $\tan\delta_{C}=2.1\pm 0.2\times 10^{\rm{-3}}$,
and $\tan\delta_{L}(1\,\text{GHz})=8.9\pm 1.8\times 10^{\rm{-6}}$. The
extracted $\tan\delta_{C}$ is consistent with the expected value for amorphous
SiOx and the low and high frequency flux noise is similar to previous
experiments Lanting et al. (2011).

Figure 4:  Variation in noise parameters for the 27 qubits calculated using
(14) and (15) with $\tan\delta_{L}(1\,\text{GHz})$ in (d). The parameters are
calculated from fits to data sets that are similar to Fig. 3.

The approximate model used to fit the data of Fig. 3 is obtained under the
assumptions listed in appendix D. A more accurate model uses the CJJ rf-SQUID
Hamiltonian (17) to calculate the tunneling amplitudes and noise broadening
parameters, as described in appendices A to C and also in [21]. In Fig. 5,
we compare the simplified model to the full model using the same noise
parameters found for Fig. 3 and the qubit parameters reported above. We find
an overall good agreement between the two models, with the simplified model
resulting in a $\sim$5% better $\chi^{2}$ value due to uncertainty in qubit
parameters. The agreement between the models gives us confidence in the noise
parameters extracted from the approximate model, which requires significantly
less computational resources.

##  IV Conclusion

We introduce a hybrid noise model for macroscopic resonant tunneling in rf-
SQUID flux qubits. The model includes contributions of low and high frequency
flux noise as well as charge noise. We fit the experimentally measured MRT
rates to the model and find good agreement over a dynamic range of four orders
of magnitude. Each noise component generates a characteristic line-shape
broadening that is captured by the fit. This allows the noise sources to be
uniquely identified and quantified. The ability to extract information about
different sources of noise in a single experimental setting and in-situ on the
quantum processor is an important step towards understanding the origin of the
measured noise and providing an indication on how to reduce it. This will
ultimately be crucial for the development of quantum computers.

###### Acknowledgements.

We acknowledge fruitful discussions with Richard Harris. We also thank Joel
Pasvolsky for a careful reading of the paper.

## Appendix

In the following appendices we provide details of the theoretical model used
in the main text. While in the main text we have used the full expressions,
for simplicity we will use $\hbar=k_{B}=1$ in the appendices.

##  Appendix A rf-SQUID Hamiltonian

A simplified version of a compound Josephson junction (CJJ) flux qubit Harris
et al. (2010) is sketched in Fig. 1(a). It has two superconducting loops, main
and CJJ loops, with flux degrees of freedom $\Phi$ and $\Phi_{\rm CJJ}$,
subject to external flux biases $\Phi^{x}$ and $\Phi_{\rm CJJ}^{x}$,
respectively. The Hamiltonian of such an rf-SQUID is written as

$\displaystyle H_{S}={q^{2}\over 2C}+{q_{\rm CJJ}^{2}\over 2C_{\rm
CJJ}}+U(\Phi,\Phi_{\rm CJJ})$ (17)

where $C$ and $C_{\rm CJJ}$ are parallel and series combinations of the
junction capacitances, $q$ and $q_{\rm CJJ}$ are the sum and difference of the
charges stored in the capacitors respectively, and

$\displaystyle U(\Phi,\Phi_{\rm
CJJ})=\frac{(\Phi{-}\Phi^{x}{+}\Phi_{0}/2)^{2}}{2L}+\frac{(\Phi_{\rm
CJJ}{-}\Phi_{\rm CJJ}^{x})^{2}}{2L_{\rm CJJ}}$ $\displaystyle-
E_{J}\cos\left(\frac{\pi\Phi_{\rm
CJJ}}{\Phi_{0}}\right)\cos\left(\frac{2\pi\Phi}{\Phi_{0}}\right)\qquad$ (18)

is a 2-dimensional potential with $L$ and $L_{\rm CJJ}$ being the inductances
of the two loops, and $\Phi_{0}{=}\,h/2e=\pi\hbar/e$ is the flux quantum. We
have assumed symmetric Josephson junctions forming the CJJ loop, with a total
critical current $I_{C}=2\pi E_{J}/\Phi_{0}$ through both junctions. Flux and
charge degrees of freedom satisfy commutation relations: $[\Phi,q]=[\Phi_{\rm
CJJ},q_{\rm CJJ}]=i\hbar$.

The CJJ loop typically has a small inductance, $L_{\rm CJJ}\ll L$, making
dynamics of $\Phi_{\rm CJJ}$ much faster than $\Phi$. Therefore, the qubit’s
quantum properties is dominantly determined by tunneling in the $\Phi$
direction. The environment also mostly affects the qubit via $\Phi$ and $q$
degrees of freedom. We therefore write the interaction Hamiltonian as

$H_{\rm int}=-{\Phi-\Phi^{x}+\Phi_{0}/2\over L}\,\delta\Phi-{q\over C}\,\delta
q\equiv-I\delta\Phi-V\delta q,$ (19)

where $\delta\Phi$ and $\delta q$ are flux and charge noise operators,
respectively, and

$I={\Phi-\Phi^{x}+\Phi_{0}/2\over L},\qquad V={q\over C}=-{i\hbar\over
C}{\partial\over\partial\Phi}$ (20)

are loop current and junction voltage operators. The flux and charge noises
are described by noise spectral densities

$\displaystyle S_{\Phi}(\omega)\\!\\!$ $\displaystyle=$
$\displaystyle\\!\\!\int dt\,e^{i\omega
t}\langle\delta\Phi(t)\delta\Phi(0)\rangle,$ (21) $\displaystyle
S_{q}(\omega)\\!\\!$ $\displaystyle=$ $\displaystyle\\!\\!\int dt\,e^{i\omega
t}\langle\delta q(t)\delta q(0)\rangle.$ (22)

where $\langle\dots\rangle$ represents averaging over environmental degrees of
freedom.

Experiments are performed when $U(\Phi,\Phi_{\rm CJJ})$ forms a double-well
potential along the $\Phi$ direction, with a large barrier between the wells.
The lowest energy states in each well are then metastable with small
amplitudes of tunneling to states in the opposite well. At the degeneracy
point, $\Phi^{x}=0$, the minima in the two wells align. We follow the
procedure described in [22] to determine the metastable states $\ket{n}$. We
divide the Hilbert space into two subspaces with
$\Phi{-}\Phi^{x}{+}\Phi_{0}/2<0$ and $\Phi{-}\Phi^{x}{+}\Phi_{0}/2>0$
corresponding to the two wells. We then partially diagonalize the Hamiltonian
in each subspace to determine $\ket{n}$. The system Hamiltonian (17) can now
be written in this new basis as

$H_{\rm S}=\sum_{n}E_{n}\ket{n}\bra{n}-{1\over 2}\sum_{m\neq
n}\Delta_{mn}\ket{m}\bra{n},$ (23)

where

$\displaystyle E_{n}=\langle n|H_{S}|n\rangle,\qquad\Delta_{mn}=-2\langle
m|H_{S}|n\rangle.$ (24)

The interaction Hamiltonian in this representation becomes

$\displaystyle H_{\rm int}=-\sum_{m,n}(I_{mn}\delta\Phi+V_{mn}\delta
q)\ket{m}\bra{n},$ (25)

where

$\displaystyle I_{mn}=\bra{m}I\ket{n},\qquad V_{mn}=\bra{m}V\ket{n}.$ (26)

We numerate states in the left (right) well with even (odd) integers (see Fig.
1b). Due to the construction of the basis, we have $\Delta_{mn}=0$ within each
well (between two even or two odd states) and $I_{mn}=0$ for every pair of
states in opposite wells (for odd $m{+}n$). Also, since $\ket{n}$ is
delocalized in charge, we expect $V_{nn}=0$ for all $n$. We define persistent
current $I_{P}$ as the expectation value of current in the lowest state of
each well when the rf-SQUID is at the degeneracy:

$\displaystyle I_{P}=\frac{I_{11}-I_{00}}{2}.$ (27)

In principle, persistent current is bias-dependent, but the dependence is
expected to be weak.

##  Appendix B Noise spectral density

Both flux noise and charge noise affect the shape of resonant tunneling peaks.
In principle, frequency dependence of flux noise can be different at low and
high frequencies. We therefore write flux noise as a sum of two components,

$S_{\Phi}(\omega)=S^{\Phi}_{L}(\omega)+S^{\Phi}_{H}(\omega),$ (28)

with different frequency dependencies. The low-frequency component typically
has $1/f$ type of spectrum

$S^{\Phi}_{L}(f)={A_{\Phi}^{2}\over|f/f_{0}|^{\alpha}}$ (29)

where $f_{0}=1\,$Hz, $\alpha\lesssim 1$, and $A_{\Phi}$ is typically of order
of a few $\mu\Phi_{0}/\sqrt{Hz}$. In practice, low-frequency noise affects the
MRT line-shape through its rms value

$W_{\Phi}=\sqrt{\int{d\omega\over 2\pi}S^{\Phi}_{L}(\omega)},$ (30)

which can also be expressed in units of energy

$W=2I_{P}W_{\phi}.$ (31)

The relation between $W_{\Phi}$ and $A_{\Phi}$ can be non-trivial. Modeling
this relation requires knowledge of accurate noise frequency dependence at
intermediate frequencies and proper introduction of integration bounds. We
therefore take $W_{\Phi}$ directly as an independent fitting parameter in our
model.

At high frequencies, flux noise is typically ohmic Lanting et al. (2011) with
spectral density

$S^{\Phi}_{H}(\omega)={1\over 4I_{P}^{2}}{\eta\omega
e^{-|\omega|/\omega_{c}}\over 1-e^{-\omega/T}},$ (32)

where $\eta$ is a dimensionless coupling coefficient and $\omega_{c}$ is a
cutoff frequency. We assume $\omega_{c}$ is larger than all relevant
frequencies and ignore it. We expect $S_{\Phi}(\omega)$ to be determined only
by flux noise and not by the qubit’s operation point. This means
$S_{\Phi}(\omega)$, and therefore $W_{\Phi}$, are independent of $I_{P}$,
hence

$W\propto I_{P},\qquad\eta\propto I_{P}^{2}.$ (33)

To have a quantity that is independent of $I_{P}$, we introduce inductive loss
tangent, $\tan\delta_{L}$, via

$S_{\Phi}(\omega)=2L\,{{\rm sgn}(\omega)\tan\delta_{L}(\omega)\over
1-e^{-\omega/T}}.$ (34)

The ${\rm sgn}(\omega)$ is added to make $\tan\delta_{L}$ positive. To agree
with $S_{\Phi}(\omega)$ in both low and high frequency regimes, we should have

$\displaystyle\tan\delta_{L}(\omega)\\!=\\!\left\\{\begin{array}[]{cc}(2\pi
f_{0})^{\alpha}(A_{\Phi}^{2}/2LT)|\omega|^{1-\alpha},&|\omega|\ll T\\\ &\\\
(\eta/8LI_{P}^{2})|\omega|,&|\omega|\gtrsim T\end{array}\right.\\!.\ \ \ $
(38)

Notice that $\tan\delta_{L}$ is independent of $I_{P}$. It is also insensitive
to the rf-SQUID geometry if $S_{\Phi}\propto L$, which is the case if the
length of the qubit wire is changed without changing its width Lanting et al.
(2009).

Similar to flux noise, charge noise also has $1/f$ spectral density at low-
frequencies

$S_{q}(f)={A_{q}^{2}\over|f/f_{0}|^{\alpha_{q}}},\qquad f_{0}=1\,\text{Hz}$
(39)

with $\alpha_{q}\approx 1$ and $A_{q}\sim 10^{-2}$ \-
$10^{-4}\,e/\sqrt{\text{Hz}}$. The $1/f$ spectral density typically crosses
over to a different frequency dependence at higher frequencies, which is
likely to be ohmic Astafiev et al. (2004); Yan et al. (2016b). It is common to
express charge noise in terms of a capacitive loss tangent, $\tan\delta_{C}$,
that characterises the quality of the dielectrics and two-level fluctuators in
the environment, independent of the qubit. We therefore define

$S_{q}(\omega)=2C\,{{\rm sgn}(\omega)\tan\delta_{C}\over 1-e^{-\omega/T}}.$
(40)

The sign function ${\rm sgn}(\omega)$ is needed to make the numerator
antisymmetric while keeping $\tan\delta_{C}$ positive. At low frequencies,
$S_{q}(\omega)\to 2CT\tan\delta_{C}/|\omega|$, which reduces to (39) with
$\alpha_{q}=1$ if

$\tan\delta_{C}=2\pi{E_{C}\over T}\left({A_{q}\over
e/\sqrt{\text{Hz}}}\right)^{2}$ (41)

where $E_{C}=e^{2}/2C$ is the charging energy. At high frequencies, (40) leads
to $S_{q}(\omega)\propto\tan\delta_{C}$, which requires
$\tan\delta_{C}\propto\omega$ for ohmic spectral density. We therefore expect
$\tan\delta_{C}$ to be constant at low frequencies with a crossover to a
different frequency dependence at high frequencies.

##  Appendix C Macroscopic resonant tunneling

Our goal is to calculate the rate of incoherent tunneling between the two
wells. Suppose at time $t=0$ the rf-SQUID is initialized in state $\ket{0}$
with probability $P_{0}=1$. The probability $P_{0}$ will decrease with time as
the system tunnels to states in the opposite well. We define the MRT
transition rate as:

$\Gamma=-\left[{dP_{0}\over dt}\right]_{t=0}.$ (42)

In general transition out of state $\ket{0}$ happens via tunneling to more
than one state in the target well. We therefore write

$\Gamma(\epsilon)=\sum_{n}\Gamma_{0n}(\epsilon-\omega_{n1}),$ (43)

where $\Gamma_{0n}$ is the transition rate from the initial state $\ket{0}$ to
state $\ket{n}$ in the target well, $\epsilon=E_{0}-E_{1}$ is the energy bias
from the degeneracy point, and $\omega_{nm}=E_{n}-E_{m}$. Each $\Gamma_{0n}$
can be calculated independently. In the next two subsections, we describe the
zeroth and first peak $\Gamma_{01}$ and $\Gamma_{03}$.

###  C.1 Tunneling between the lowest energy states

To calculate $\Gamma_{01}$, we need to consider only two states $\ket{0}$ and
$\ket{1}$, corresponding to the ground states in the left and right wells,
respectively. We can therefore represent the system Hamiltonian in terms of
Pauli matrices

$\displaystyle\sigma_{x}$ $\displaystyle=$
$\displaystyle\ket{0}\bra{1}+\ket{1}\bra{0},$ $\displaystyle\sigma_{z}$
$\displaystyle=$ $\displaystyle\ket{1}\bra{1}-\ket{0}\bra{0}.$ (44)

The effective Hamiltonian in this subspace becomes

$\displaystyle
H_{S}=-\frac{\epsilon}{2}\,\sigma_{z}-\frac{\Delta_{01}}{2}\,\sigma_{x}.$ (45)

It can be shown that

$\displaystyle\epsilon\approx 2I_{P}\,\Phi^{x},$ (46)

with the external flux $\Phi^{x}$ measured relative to the degeneracy point.
Substituting the current operator $I=I_{P}\sigma_{z}$ into Eqs. (25), the
interaction Hamiltonian for flux noise becomes

$\displaystyle H_{\rm
int}^{\Phi}=-I_{P}\,\delta\Phi\,\sigma_{z}=-\frac{1}{2}\,Q\,\sigma_{z},$ (47)

with the noise operator

$\displaystyle Q=2\,I_{P}\,\delta\Phi.$ (48)

The effective Hamiltonian of the two-state system describing the rf-SQUID
coupled to flux noise environment is therefore given by

$\displaystyle
H=-\frac{\epsilon}{2}\,\sigma_{z}-\frac{\Delta_{01}}{2}\,\sigma_{x}-\frac{1}{2}\,Q\,\sigma_{z}.$
(49)

We introduce the spectral density corresponding to operator $Q$ as

$\displaystyle S_{Q}(\omega)=\int dt\,e^{i\omega t}\langle
Q(t)Q(0)\rangle=4I_{p}^{2}S_{\Phi}(\omega).$ (50)

Like $S_{\Phi}(\omega)$, we can decompose this spectral density into low and
high frequency components:
$S_{Q}(\omega)=S^{Q}_{L}(\omega)+S^{Q}_{H}(\omega)$. It was shown in [20]
that the MRT transition rate can be expressed by a convolution integral

$\displaystyle\Gamma_{01}(\epsilon)\approx{\Delta_{01}^{2}\over
4}\int{d\omega\over 2\pi}G_{L}(\epsilon{-}\omega)G_{H}(\omega).$ (51)

The effect of low-frequency flux noise is captured by the Gaussian envelope

$G_{L}(\omega)={\sqrt{2\pi}\over
W}\exp\left\\{{-(\omega-\epsilon_{p})^{2}\over 2W^{2}}\right\\},$ (52)

where

$\displaystyle W^{2}\\!\\!$ $\displaystyle=$
$\displaystyle\\!\\!\int{d\omega\over 2\pi}S^{Q}_{L}(\omega),$ (53)
$\displaystyle\epsilon_{p}\\!\\!$ $\displaystyle=$ $\displaystyle\\!\\!{\cal
P}\int{d\omega\over 2\pi}{S^{Q}_{L}(\omega)\over\omega}.$ (54)

Here, ${\cal P}$ represents the principal value integral. Fluctuation-
dissipation theorem requiresAmin and Averin (2008)

$W^{2}=2T\epsilon_{p}.$ (55)

High frequency noise affects the peak through a Lorentzian-like envelope
function

$G_{H}(\omega)\approx{S^{Q}_{H}(\omega)\over\omega^{2}+[{1\over
2}S^{Q}_{H}(\omega)]^{2}},$ (56)

with

$\displaystyle S^{Q}_{H}(\omega)={\eta\omega\over 1-e^{-\omega/T}}.$ (57)

At small frequencies, $\omega\ll T$, we have $S^{Q}_{H}(\omega)=\eta T$,
therefore, the denominator of (56) can be written as $\omega^{2}+\gamma^{2}$,
where $\gamma\equiv\eta T/2$. At large frequencies, $\omega\gg T$, we have
$S^{Q}_{H}(\omega)=\eta\omega$, and the denominator of (56) becomes
$(1+\eta^{2}/4)\omega^{2}\approx\omega^{2}$ for $\eta\ll 1$. Therefore, we can
express the high frequency envelope function in terms of the broadening
parameter $\gamma$ as

$G_{H}(\omega)\approx{2\gamma\omega/T\over(\omega^{2}+\gamma^{2})(1-e^{-\omega/T})}.$
(58)

The two parameters $W$ and $\gamma$ measure the width of the envelope
functions $G_{L}$ and $G_{H}$ in energy units, respectively. Each envelope
function approximately becomes a delta function when the width goes to zero.
In the absence of high frequency noise, $\gamma\to 0$, we have

$\displaystyle\Gamma_{01}(\epsilon)\approx{\Delta_{01}^{2}\over
4}G_{L}(\epsilon),$ (59)

in agreement with [19 and 16]. Similarly, when low-frequency noise is absent
($W,\epsilon_{p}\to 0$), the Gaussian function (52) becomes $\delta$-function,
and we obtain

$\Gamma_{01}(\epsilon)\approx{\Delta_{01}^{2}\over
4}{S^{Q}_{H}(\epsilon)\over\epsilon^{2}+\gamma^{2}}.$ (60)

For $\epsilon>\gamma$ we get

$\Gamma_{01}(\epsilon)\approx{\Delta_{01}^{2}\over
4\epsilon^{2}}S^{Q}_{H}(\epsilon)$ (61)

in agreement with the Bloch-Redfield theory. For small biases, we obtain the
Lorentzian relaxation rate expected for white noise Amin and Averin (2008)

$\Gamma_{01}(\epsilon)\approx{\Delta_{01}^{2}\over
2}{\gamma\over\epsilon^{2}+\gamma^{2}}.$ (62)

Therefore, the convolution form (51), with envelope functions (52) and (58),
gives correct results in all limiting regimes. It can also be numerically
shown that it agrees well with the exact results in other regimes as long as
$\eta\ll 1$. One advantage of the convolution form is that it separates
contributions of low and high frequency noise into two separate envelope
functions. It is therefore possible to study the effect of each noise
separately and calculate the corresponding envelope function. We will use this
property to determine the effect of intrawell relaxation on multi-level MRT
peaks in the next subsection.

###  C.2 Tunneling to a higher energy state

We now consider multi-level MRT peaks when tunneling happens to a higher
energy state in the target well. For simplicity, we consider transition to the
second level in the target well, i.e., $\Gamma_{03}$. As before, we assume
that the system is initialized in state $\ket{0}$. Incoherent tunneling from
$\ket{0}$ to $\ket{3}$ is affected by the flux noise the same way as discussed
in the previous subsection. The broadening due to low-frequency and high-
frequency noise is captured by $G_{L}(\omega)$ and $G_{H}(\omega)$,
respectively, with minor changes to the parameters that we shall mention
below. However, since $\ket{3}$ is an excited state within the target well,
the system will quickly relax to state $\ket{1}$, in a time scale much shorter
than the incoherent tunneling rate. The uncertainty in energy $E_{3}$ due to
the intrawell relaxation creates an additional broadening. Such a broadening
was introduced in [ 19], but the resulting transition peak was symmetric
around its center, violating the detailed balance needed to reach Boltzmann
distribution in thermal equilibrium. Here, we provide a simple calculation of
the broadening effect in a way that satisfies detailed balance. A more
complete derivation is provided in [21].

As we mentioned before, the broadening effect of every component of noise can
be calculated independently and combined together through a convolution
integral. The combined transition rate is

$\Gamma_{03}(\epsilon)={\Delta_{03}^{2}\over
4}\,G_{03}(\epsilon-\omega_{31}),$ (63)

with

$G_{03}(\omega)=\\!\\!\\!\int{d\omega^{\prime}\over
2\pi}{d\omega^{\prime\prime}\over
2\pi}G_{L}(\omega{-}\omega^{\prime})G_{H}(\omega^{\prime}{-}\omega^{\prime\prime})G_{R}(\omega^{\prime\prime}),\quad$
(64)

where $G_{L}(\omega)$ and $G_{H}(\omega)$ are given by (52) and (58),
respectively, and $G_{R}(\omega)$ is a peaked function capturing the
broadening due to the intrawell relaxation. As before, $\epsilon=E_{0}-E_{1}$
is the energy bias measured from the rf-SQUID degeneracy point. We therefore
have

$\omega_{03}=\epsilon-\omega_{31}\approx 2I_{P}^{03}(\Phi^{x}-\Phi^{x}_{31}).$
(65)

Here, we define $\Phi^{x}_{31}$ as the value of the external flux $\Phi^{x}$
when energy states $\ket{0}$ and $\ket{3}$ are in resonance, and introduce a
generalized persistent current

$I_{P}^{0n}=(I_{nn}-I_{00})/2,$ (66)

which captures state dependence of the current matrix element $I_{nn}$. As
pointed out before, the broadening widths due to both low and high frequency
noise are functions of the persistent current: $W\propto I^{0n}_{P},\
\gamma\propto(I_{P}^{0n})^{2}$. One should therefore rescale these parameters
in $G_{L}(\omega)$ and $G_{H}(\omega)$ according to (66). However, since the
state dependence of the persistent current is expected to be weak, we assume
$I_{P}^{0n}\approx I_{P}$ and neglect these small corrections.

We obtain $G_{R}(\omega)$ by calculating the transition rate when the only
broadening effect is due to intrawell relaxation, i.e., other noise
contributions are turned off ($W=\gamma=0$). To simplify the calculation we
focus on a three-state system described by

$H_{S}=\sum_{n=0,1,3}E_{n}\ket{n}\bra{n}-{\Delta_{03}\over
2}(\ket{0}\bra{3}+\ket{3}\bra{0})$ (67)

with interaction Hamiltonian

$H_{\rm int}=-Q_{R}\ket{3}\bra{1}+h.c.$ (68)

where

$Q_{R}=I_{31}\delta\Phi+V_{31}\delta q$ (69)

provides coupling to flux and charge noise. Notice that the interaction
Hamiltonian can only cause transition between states $\ket{3}$ and $\ket{1}$.
The intrawell relaxation rate can be calculated using Bloch-Redfield theory

$\Gamma_{31}(\omega_{31})=S_{R}(\omega_{31}),$ (70)

where

$\displaystyle S_{R}(\omega)=\int dt\,e^{i\omega t}\langle
Q_{R}(t)Q_{R}(0)\rangle$ (71)

is the environment spectral density corresponding to $Q_{R}$ defined in (69).
When this relaxation is strong, it is not possible to separate interwell
tunneling and intrawell relaxation as two independent processes. We therefore
combine them into a single quantum mechanical process that creates transition
from $\ket{0}$ to $\ket{1}$ mediated by state $\ket{3}$ (via virtual
transition).

Using perturbation expansion in $\Delta_{03}/\omega_{03}\ll 1$, we diagonalize
Hamiltonian (67) to obtain

$H_{S}=\sum_{n=0,1,3}\tilde{E}_{n}\ket{\tilde{n}}\bra{\tilde{n}},$ (72)

where $\tilde{E}_{n}$ are perturbed energies and

$\displaystyle\ket{\tilde{0}}\\!\\!$ $\displaystyle=$
$\displaystyle\\!\\!\ket{0}+{\Delta_{03}\over 2\omega_{03}}\ket{3},$ (73)
$\displaystyle\ket{\tilde{1}}\\!\\!$ $\displaystyle=$
$\displaystyle\\!\\!\ket{1},$ (74) $\displaystyle\ket{\tilde{3}}\\!\\!$
$\displaystyle=$ $\displaystyle\\!\\!\ket{3}-{\Delta_{03}\over
2\omega_{03}}\ket{0},$ (75)

are perturbed eigenstates. The interaction Hamiltonian in this basis is

$\displaystyle H_{\rm
int}=-Q_{R}\left(\ket{\tilde{3}}\bra{\tilde{1}}+{\Delta_{03}\over
2\omega_{03}}\ket{\tilde{0}}\bra{\tilde{1}}\right)+h.c.$ (76)

Since there is no off-diagonal term between $\ket{\tilde{0}}$ and
$\ket{\tilde{3}}$ in Hamiltonians (72) and (76), we can now remove state
$\ket{\tilde{3}}$ from consideration. Introducing Pauli matrices

$\displaystyle\sigma_{x}$ $\displaystyle=$
$\displaystyle\ket{\tilde{0}}\bra{\tilde{1}}+\ket{\tilde{1}}\bra{\tilde{0}},$
$\displaystyle\sigma_{z}$ $\displaystyle=$
$\displaystyle\ket{\tilde{1}}\bra{\tilde{1}}-\ket{\tilde{0}}\bra{\tilde{0}},$
(77)

we obtain the familiar two-state Hamiltonian

$\displaystyle H=-{1\over 2}\epsilon\,\sigma^{z}-{\Delta_{03}\over
2\omega_{03}}Q_{R}\,\sigma^{x},$ (78)

where $\epsilon=\omega_{01}=\omega_{03}{+}\omega_{31}$ is the energy bias from
qubit degeneracy. Here, we ignore second order corrections to energies
$\tilde{E}_{n}$. The $\sigma^{x}$ term can be treated using Bloch-Redfield
formalism to obtain

$\Gamma_{03}(\omega_{03})={\Delta_{03}^{2}\over
4\omega_{03}^{2}}S_{R}(\omega_{03}{+}\omega_{31}).$ (79)

Equation (79) is divergent at $\omega_{03}=0$. We can remove the divergence
the same way as in (56) by writing

$\Gamma_{03}(\omega_{03})={\Delta_{03}^{2}\over 4}G_{R}(\omega_{03}),$ (80)

with the envelope function

$G_{R}(\omega)={\Gamma_{31}(\omega{+}\omega_{31})\over\omega^{2}+[\Gamma_{31}(\omega{+}\omega_{31})/2]^{2}},$
(81)

where $\Gamma_{31}(\omega)=S_{R}(\omega)$ is a frequency dependent intrawell
relaxation rate. As we mentioned before, intrawell relaxation can be caused by
both flux and charge noise. Therefore, $S_{R}(\omega)$ could be a sum of two
components

$S_{R}(\omega)=I_{31}^{2}S_{\Phi}(\omega)+V_{31}^{2}S_{q}(\omega),$ (82)

where $I_{n1}$ and $V_{n1}$ are defined in (26). However, with our noise
parameters, contribution of flux noise is negligible compared to charge noise.
We use (40) for spectral density of charge noise, with magnitude of noise
characterized by a constant loss tangent $\tan\delta_{C}$. To avoid
singularity at $\omega=0$, we need to replace ${\rm sgn}(\omega)$ with a
smoother function. The center of the MRT peaks is dominantly broadened by the
low-frequency flux noise, with almost no effect from low-frequency components
of charge noise. We therefore choose

${\rm sgn}(\omega)\to\tanh(\omega/T)$ (83)

which gives maximally flat spectrum $S_{q}(\omega)$ near $\omega=0$ without
additional fitting parameters. We therefore write

$\Gamma_{31}(\omega)=\zeta\,{\tanh(\omega/T)\over 1-e^{-\omega/T}}.$ (84)

where

$\zeta=2CV_{31}^{2}\tan\delta_{C}.$ (85)

is now a frequency independent parameter characterizing the width of
$G_{R}(\omega)$ in (81). $\zeta$ measures the (frequency independent)
intrawell relaxation out of state $\ket{3}$. If we ignore frequency dependence
of (84) we recover the symmetric result of [19]. The relaxation rate then
would not satisfy detailed balance and cannot explain the experimental
results.

For higher energy states in the target well, there are more than one channel
of relaxation. Therefore, the effective broadening becomes larger with higher
energy. A more general and rigorous derivation of (81) is provided in [21].

##  Appendix D Simplified model

Figure 5:  Comparison between the simplified model (black) used to fit the
experimental data (red symbols) and the full model (blue). The simplified
model uses approximations (detailed in appendix D) to reduce computational
time. These approximations are relaxed by using the full model (see main text)
that uses the extracted noise parameters from the fit. The two models agree
very well with each other and with the experimental data.

The formalism described in the previous sections was obtained using
Hamiltonian (23) and interaction Hamiltonian (25). These Hamiltonians were
themselves obtained from the rf-SQUID Hamiltonian (17) after partial
diagonalization. The procedure to extract Hamiltonian parameters is time
consuming and requires accurate knowledge of circuit parameters, such as
inductance, capacitance, and critical current. In this section we introduce an
approximation to this model that allows fitting to experimental data and
extracting noise parameters without diagonalization or knowledge of rf-SQUID
parameters. The assumptions behind this approximation are as follows:

  1. 1.

Energy bias, $\epsilon=E_{0}-E_{1}$, is a linear function of the applied flux
$\Phi^{x}$ (as in (46)) over the experimental range.

  2. 2.

Persistent current $I_{P}$ has weak bias dependence and is measured
independently.

  3. 3.

Tunneling amplitudes $\Delta_{0n}$ have negligible bias dependence.

  4. 4.

Current matrix element $I_{nn}$ ($\approx I_{P}$) is weakly dependent on state
$\ket{n}$, therefore, noise parameters $W$ and $\gamma$ are the same for all
$\Gamma_{0n}(\epsilon)$.

  5. 5.

Capacitive loss tangent $\tan\delta_{C}$ is constant over the range of
frequencies that matter for intrawell relaxation (close to $\omega_{31}$).

  6. 6.

Inter- (intra-) well transitions are dominantly affected by flux (charge)
noise.

With these assumptions, one can fit the model to the experimental data using
six fitting parameters, $\Delta_{01}$, $\Delta_{03}$, $\omega_{31}$, $W$,
$\gamma$, and $\zeta$, with no need for diagonalization. Note that all these
parameters are in energy units. However, the potential tilt is applied to the
rf-SQUID via an external flux bias ($\epsilon=2I_{P}\Phi^{x}$). Therefore, the
MRT peaks are measured as functions of flux (not energy) bias. The distance
between the MRT peaks, $\Phi^{x}_{31}={\omega_{31}/2I_{P}}$, is also directly
measured in flux units. It is therefore convenient to express noise parameters
directly in flux units:

$W_{\Phi}={W\over 2I_{P}},\qquad\gamma_{\Phi}={\gamma\over
2I_{P}},\qquad\zeta_{\Phi}={\zeta\over 2I_{P}}.$ (86)

Each of these broadening parameters characterizes one component of noise.
$W_{\Phi}$ measures the r.m.s. value of the low-frequency flux noise, and
$\gamma_{\Phi}$ measures the magnitude of the high frequency flux noise. From
$\gamma_{\Phi}$, the dimensionless ohmic coefficient and the inductive loss
tangent can be calculated:

$\eta={4I_{P}\gamma_{\Phi}\over
T},\qquad\tan\delta_{L}(\omega)={\gamma_{\Phi}\over 4LI_{P}}\left|{\omega\over
T}\right|.$ (87)

Finally, the broadening parameter $\zeta_{\Phi}$ characterizes the charge
noise. When the potential barrier is high, the bottom of the target well can
be approximated by a parabola. The lowest energy levels inside the well can
therefore be obtained using a harmonic oscillator model. One can then show
that

$V_{31}\approx\sqrt{\omega_{31}\over 2C}.$ (88)

Using (84) and (85), and assuming $\omega_{31}\gg T$, we obtain

$\Gamma_{31}(\omega_{31})\approx\zeta\approx\omega_{31}\tan\delta_{C}.$ (89)

This is what one expects for relaxation in a harmonic oscillator. Converting
to flux units, we obtain

$\tan\delta_{C}\approx{\zeta_{\Phi}\over\Phi^{x}_{31}}.$ (90)

As usual, loss tangent is the ratio of peak broadening and oscillation
frequency, both measured in flux units.

## References

  * de Leon et al. (2021) N. P. de Leon, K. M. Itoh, D. Kim, K. K. Mehta, T. E. Northup, H. Paik, B. Palmer, N. Samarth, S. Sangtawesin, and D. Steuerman, Science 372 (2021).
  * Yan et al. (2016a) F. Yan, S. Gustavsson, A. Kamal, J. Birenbaum, A. P. Sears, D. Hover, T. J. Gudmundsen, D. Rosenberg, G. Samach, S. Weber, et al., Nature communications 7, 1 (2016a).
  * Nguyen et al. (2019) L. B. Nguyen, Y.-H. Lin, A. Somoroff, R. Mencia, N. Grabon, and V. E. Manucharyan, Physical Review X 9, 041041 (2019).
  * Place et al. (2021) A. P. Place, L. V. Rodgers, P. Mundada, B. M. Smitham, M. Fitzpatrick, Z. Leng, A. Premkumar, J. Bryon, A. Vrajitoarea, S. Sussman, et al., Nature communications 12, 1 (2021).
  * Siddiqi (2021) I. Siddiqi, Nat Rev Mater (2021).
  * Bialczak et al. (2007) R. C. Bialczak, R. McDermott, M. Ansmann, M. Hofheinz, N. Katz, E. Lucero, M. Neeley, A. O’connell, H. Wang, A. Cleland, et al., Physical review letters 99, 187006 (2007).
  * Lanting et al. (2009) T. Lanting, A. Berkley, B. Bumble, P. Bunyk, A. Fung, J. Johansson, A. Kaul, A. Kleinsasser, E. Ladizinsky, F. Maibaum, et al., Physical Review B 79, 060509 (2009).
  * Quintana et al. (2017) C. Quintana, Y. Chen, D. Sank, A. Petukhov, T. White, D. Kafri, B. Chiaro, A. Megrant, R. Barends, B. Campbell, et al., Physical review letters 118, 057702 (2017).
  * Braumüller et al. (2020) J. Braumüller, L. Ding, A. P. Vepsäläinen, Y. Sung, M. Kjaergaard, T. Menke, R. Winik, D. Kim, B. M. Niedzielski, A. Melville, et al., Physical Review Applied 13, 054079 (2020).
  * Koch et al. (2007) R. H. Koch, D. P. DiVincenzo, and J. Clarke, Physical review letters 98, 267003 (2007).
  * Anton et al. (2013) S. M. Anton, J. S. Birenbaum, S. R. O’Kelley, V. Bolkhovsky, D. A. Braje, G. Fitch, M. Neeley, G. C. Hilton, H.-M. Cho, K. D. Irwin, et al., Phys. Rev. Lett. 110, 147002 (2013).
  * Lanting et al. (2014) T. Lanting, M. H. Amin, A. J. Berkley, C. Rich, S.-F. Chen, S. LaForest, and R. de Sousa, Phys. Rev. B 89, 014503 (2014).
  * Lanting et al. (2020) T. Lanting, M. Amin, C. Baron, M. Babcock, J. Boschee, S. Boixo, V. Smelyanskiy, M. Foygel, and A. Petukhov, arXiv preprint arXiv:2003.14244 (2020).
  * Martinis et al. (2005) J. M. Martinis, K. B. Cooper, R. McDermott, M. Steffen, M. Ansmann, K. Osborn, K. Cicak, S. Oh, D. P. Pappas, R. W. Simmonds, et al., Physical review letters 95, 210503 (2005).
  * Müller et al. (2019) C. Müller, J. H. Cole, and J. Lisenfeld, 82, 124501 (2019).
  * Harris et al. (2008) R. Harris, M. Johnson, S. Han, A. Berkley, J. Johansson, P. Bunyk, E. Ladizinsky, S. Govorkov, M. Thom, S. Uchaikin, et al., Physical review letters 101, 117003 (2008).
  * Lanting et al. (2011) T. Lanting, M. Amin, M. Johnson, F. Altomare, A. Berkley, S. Gildert, R. Harris, J. Johansson, P. Bunyk, E. Ladizinsky, et al., Physical Review B 83, 180502 (2011).
  * Harris et al. (2010) R. Harris, J. Johansson, A. Berkley, M. Johnson, T. Lanting, S. Han, P. Bunyk, E. Ladizinsky, T. Oh, I. Perminov, et al., Physical Review B 81, 134510 (2010).
  * Amin and Averin (2008) M. H. Amin and D. V. Averin, Physical review letters 100, 197001 (2008).
  * Smirnov and Amin (2018) A. Y. Smirnov and M. H. Amin, New Journal of Physics 20, 103037 (2018).
  * Smirnov et al. (2022) A. Y. Smirnov, A. Whiticar, and M. H. Amin, arXiv preprint arXiv:2209.10605 (2022).
  * Amin et al. (2013) M. H. Amin, N. G. Dickson, and P. Smith, Quantum information processing 12, 1819 (2013).
  * Astafiev et al. (2004) O. Astafiev, Y. A. Pashkin, Y. Nakamura, T. Yamamoto, and J.-S. Tsai, Physical review letters 93, 267007 (2004).
  * Yan et al. (2016b) F. Yan, S. Gustavsson, A. Kamal, J. Birenbaum, A. P. Sears, D. Hover, T. J. Gudmundsen, D. Rosenberg, G. Samach, S. Weber, et al., Nature communications 7, 1 (2016b).