NASA Technical Reports Server (NTRS) 20160001868: Computational Role of Tunneling in a Programmable Quantum Annealer

4 MB·English

by  NASA Technical Reports Server (NTRS)

Computational Role of Tunneling in a Programmable Quantum Annealer
Sergio Boixo1, Vadim N. Smelyanskiy2, Alireza Shabani1, Sergei V. Isakov1, Mark
Dykman3, Mohammad Amin4, Masoud Mohseni1, Vasil S. Denchev1, Hartmut Neven1
1Google
2NASA Ames
3Physics and Astronomy Department, Michigan State University
4D-Wave Systems
Quantum tunneling is a phenomenon in which a quantum state tunnels through energy barriers
abovetheenergyofthestateitself[1]. Tunnelinghasbeenhypothesizedasanadvantageousphysical
resourceforoptimization[2–6]. Herewepresentthefirstexperimentalevidenceofacomputational
role of multiqubit quantum tunneling in the evolution of a programmable quantum annealer. We
developed a theoretical model based on a NIBA Quantum Master Equation to describe the multi-
qubitdissipativecotunnelingeffectsunderthecomplexnoisecharacteristicsofsuchquantumdevices.
We start by considering a computational primitive, the simplest non-convex optimization problem
consisting of just one global and one local minimum. The quantum evolutions enable tunneling to
theglobalminimumwhilethecorrespondingclassicalpathsaretrappedinafalseminimum. Inour
study the non-convex potentials are realized by frustrated networks of qubit clusters with strong
intra-clustercoupling. Weshowthatthecollectiveeffectofthequantumenvironmentissuppressed
in the “critical” phase during the evolution where quantum tunneling “decides” the right path to
solution. In a later stage dissipation facilitates the multiqubit cotunneling leading to the solution
state. The predictions of the model accurately describe the experimental data from the D-Wave
II quantum annealer at NASA Ames. In our computational primitive the temperature dependence
of the probability of success in the quantum model is opposite to that of the classical paths with
thermalhopping. Specifically,weprovideananalysisofanoptimizationproblemwithsixteenqubits,
demonstrating eight qubit cotunneling that increases success probabilities. Furthermore, we report
results for larger problems with up to 200 qubits that contain the primitive as subproblems.
I. INTRODUCTION tunneling in Ref. [8] (see below).
Quantumtunnelingwasdiscoveredinthelate1920s
to explain field electron emission in vacuum tubes.
Today this phenomenon is at the core of many es-
sential technological innovations such as the tunnel
field-effect transistor, field emission displays and the
scanning tunneling microscope. Tunneling is also at
the heart of energy and charge transport in biological
and chemical processes. Recently cotunneling effects
involvingmultiplequantummechanicalparticleshave
been used to develop single electron transistors and
hypersensitive measurement instruments.
Quantum tunneling, in particular for thin but high
energy barriers, has been hypothesized as an advan-
tageous mechanism for quantum optimization [2–6].
In classical simulated annealing or cooling optimiza-
tion algorithms, the corresponding temperature pa-
rameter must be raised to overcome energy barriers.
But if there are many potential local minima with
smaller energy differences than the height of the bar-
rier, the temperature must also be lowered to distin-
guish between them so the algorithm can converge to FIG.1. QuantumannealingfunctionsA(s)andB(s). The
the global minimum. Quantum tunneling is present functionA(t)isdefinedastwicethemedianenergydiffer-
even at zero-temperature. Therefore, for some energy encebetweenthetwolowesteigenstatesoftheexperimen-
tallysuperconductingfluxqubit. ThefunctionB(s)isde-
landscapes, one might expect that quantum dynam-
finedas1.41picohenriestimesthesquareofthepersistent
ical evolutions can converge to the global minimum
current I2(s), as explained in App. IIB.
faster than classical optimization algorithms. Quan- p
tumannealing[2,3]isatechniqueinspiredbyclassical
annealingbutdesignedtotakeadvantageofquantum
tunneling. Single qubit quantum tunneling for a pro-
grammable annealer has been demonstrated experi-
mentally in Ref. [7], and two qubit cotunneling has Thestateevolutionintransversefieldquantuman-
been detected indirectly using microscopic resonant nealing is governed by a time dependent Hamiltonian
2
of the form [9] the lowest energy state should decrease with increas-
ing temperature for a quantum system but should in-
H0(s)=A(s)HD+B(s)HP (1) crease for a classical system. This is indeed observed,
HD =−(cid:88)σx (2) demonstrating eight qubit cotunneling that increases
µ success probabilities. We compare with physically
µ
plausible models of the hardware that only employ
(cid:88) (cid:88)
HP =− hµσµz − Jµνσµzσνz . (3) product states which do not allow for multiqubit tun-
µ µν neling transitions. Experimentally we find that for
this situation the D-Wave II processor returns the so-
HereHDisthedriverHamiltonian,HPistheproblem
lutionthatminimizestheenergywithhigherprobabil-
Hamiltonian whose ground state is the solution of an
ity than these models. Beyond the original 16 qubit
optimization problem of interest, {σx,σz} are Pauli
µ µ probe problem we also explore larger problems that
matrices acting on spin µ, s = t/t is the annealing
qa contain multiple weak-strong cluster pairs.
parameter, and t is the duration of the quantum
qa
annealing process. The functions A(s) and B(s) used
in the rest of the paper are shown in Fig. 1. The
II. A PRIMITIVE “PROBE” PROBLEM
Hamiltonian path H (s) describes an evolution of ef-
0 CHARACTERIZED BY A DOUBLE WELL
fective 2-level spin systems (qubits) from an initial
POTENTIAL
phase with a unique ground state to a final Hamilto-
nian with eigenstates aligned with the z quantization
A. The quantum Hamiltonian
axis. In the initial unique ground state all the qubits
arealignedwiththeeffectivetransversemagneticfield
in the x direction. See Appendix A for a more com- Thearchetypalprimitivetostudyquantumtunnel-
plete derivation of the single qubit Hamiltonian and ing is a double-well potential: two local minima sepa-
the parameteres the experimental system considered rated by an energy barrier. Our aim is to distinguish
in this paper. quantumtunnelingfromthermalactivationinamodel
Closed system quantum adiabatic evolutions gov- using classical paths. Classical paths are limited to
erned by the time-dependent Hamiltonian given in local spin vector dynamics over product states to get
Eq. (1) will arrive at the final ground state of the overtheenergybarrier. Incontrast,thesignatureofa
problem Hamiltonian if the total evolution time is quantumsystemisthatentangledstatesareavailable
large compared to the inverse minimum energy gap as well. We utilize qubit networks of the D-Wave 2
alongtheHamiltonianpath[9]. Inthispaperweshall quantumannealerchipatNASAAmestodesigntime
analyze the performance of a quantum annealing de- dependent asymmetric double well potentials where a
vice with superconducting flux qubits [7, 10, 11]. The classical path continuously connects the initial global
qubitsarecoupledinductivelyinaconnectivitygraph minimum to the final false minimum. In this way one
thatisformedbyagridofcellswithhighinternalcon- can study how the system escapes the local minimum
nectivity. The qubits are subject to interaction with and traverses the energy barrier to reach the global
theenvironmentwiththedominantnoisesourcebeing optimum. Wewillseethatquantumtunnelingresults
spin diffusion at the superconductor insulator inter- in a different final probability of success than the cor-
face [12–14]. This is known to produce control errors, respondingclassicaldynamicsoverclassicalpaths. We
energy level broadening as well as thermal excitation compare the experimental data from the device with
and relaxation [15, 16]. The noise characteristics of numerical simulations of classical paths and with the
individual qubits have been studied in macroscopic predictions of a comprehensive analytical model for
resonant tunneling experiments [17]. We show nev- dissipative multiqubit cotunneling. Based on the re-
ertheless that even under such conditions the device sults of this comparison, we establish the functional
performance can benefit from multiqubit cotunneling role of tunneling in the evolution on a programmable
of strongly interacting qubit clusters. This is of rele- quantum annealer.
vance for current programmable quantum annealers, We now detail how the double well potential and
such as the D-Wave II chip at NASA Ames. time evolution can be constructed in the case of net-
In this work we design an Ising model implementa- workgraphswithfiniteconnectivity. Wewillfocuson
tion with 16 qubits of a computational primitive, the theparticularcaseoftheso-calledChimeragraphthat
simplest non-convex optimization problem consisting connectsthequbitsinthecurrentD-WaveIIarchitec-
of just one global and one local minimum. The final ture, althoughsimilarconstructionscanbeappliedto
globalminimum canonlybereachedbytraversingan more general network architectures. We choose our
energy barrier. We develop a NIBA Quantum Master doublewellprimitiveprobeproblemtobetheonede-
Equation which takes high and low frequency noise picted in Fig. 2. We use two Chimera cells, each with
into account. Our comprehensive open quantum sys- n = 8 qubits. We find it useful during our analysis
temmodeling,withcloseagreementwithD-Wavema- to keep n explicit. We will choose equal local fields
chineoutput,demonstrateshowcotunnelingcanexist for the spins within each cell. We also choose all the
and play a functional role in presence of both Ohmic couplings to be equal and ferromagnetic. There are
and strong 1/f noise components affecting coherence n2/4intra-cellcouplingsandn/2inter-cellcouplings.
of the flux qubits. Independent of specific choices in The spins within each cell tend to move together as
the quantum models or the classical models, the pre- an homogenous cluster because flipping only one spin
diction is that the probabilities to find the system in rises the energy by an amount ∝ nJ which is much
3
J = 1 quantumannealingtheeffectofthelocalz-fieldsdom-
inatesthatoftheinter-cellIsingcouplings. According
to (7), because h and h have the opposite signs so
1 2
j=1 j=5 8 12 will the z-projections of the spins (cid:104)σz (cid:105) in the two
j,k
clusters early in the evolution.
A key observation is that in the absence of quan-
j=2 j=6 9 13
tumtunnelingandthermalhopingthespinprojection
h h = -1
1 2
of the two clusters stay opposite during the evolution
j=3 j=7 10 14 arriving to the false minimum with residual energy
relative to the global minimum equal to n(J −2h ).
1
j=4 j=8 11 15 The system will get trapped into the false minimum.
To escape it all spins in the left cluster must flip the
k=1 k=2
sign which requires going over the barrier top. At the
tippingpointofthebarriertheleftclusterhaszeroto-
talz-polarizationandthereforethebarriergrowswith
FIG. 2. The probe problem under study, consisting of
the ferromagnetic energy of the cluster (n/2)2J. For
16qubitsintwounitcellsoftheso-calledChimeragraph.
Allqubitsareferromagneticallycoupledandevolveastwo sufficiently large n the barrier height O(n2) is much
distinct qubit clusters. At the end of the annealing evo- greater than the residual energy O(n). It will be
lution the right cluster is strongly pinned downward due shown below that in certain region of the annealing
to strong local fields acting on all qubits in that cell. The parameter s all qubits in the left cluster will tunnel
local magnetic field h1 in the left cluster is weaker how- in a concerted motion under the energy barrier sepa-
ever,andservesasabifurcationparameter. Forh1 <J/2 rating the two potential wells that correspond to the
the left cluster will reverse its orientation during the an-
opposite z-polarizations of the cluster.
nealing sweep and eventually align itself with the right
Next, we discuss an approximation which reduces
cluster. Note the permutation symmetry in each column
the size of the Hamiltonian matrix for the 2 unit cell
which allows us to adopt the large spin description.
problem from 22n to (n/4 + 1)4. We introduce to-
tal spin operators for each column of a unit cell (cf.
Fig. 2)
greater than the energy of the inter-cell bond ∝ 2J.
The problem Hamiltonian is of Ising form
n/2 n
1(cid:88) 1 (cid:88)
Sα = σα , Sα = σα , (8)
HP =HP+HP+HP (4) k,1 2 j,k k,2 2 j,k
1 2 1,2 j=1 j=n/2+1
n
(cid:88) (cid:88)
HkP =−J σjz,kσjz(cid:48),k− hkσjz,k (5) where α ∈ x,y,z, and k ∈ {1,2} denotes the left and
(cid:104)j,j(cid:48)(cid:105)∈intra j=1 right Chimera cells. Because the intra-cell Hamilto-
(cid:88) nians (5) and the driver Hamiltonian are symmetric
HP =−J σz σz . (6)
1,2 j,1 j,2 with respect to qubit permutations they can be writ-
j∈inter
ten in terms of the total spin operators
The index k ∈ {1,2} denotes the Chimera cell, the HP =−4JSz Sz +2h Sz (9)
first sum in (5) goes over the intra-cell couplings de- k k,1 k,2 k k
(cid:88)
picted in Fig. 2, and the second sum goes over the HD =−2 Sx .
k,m
inter-cell couplings corresponding to j ∈(n/2+1,n) k,m=1,2
in Fig. 2; h denotes the local fields within each cell.
k
Wenotethattheinter-cellHamiltonianHP inEq.(6)
Toward the end of quantum annealing the prob- 1,2
lem Hamiltonian HP is dominating the evolution and does not possess the qubit permutation symmetry.
(cid:104)σz (cid:105) (cid:39) ±1 where (cid:104)...(cid:105) denotes a quantum mechan- However, as explained above, the qubits in each cell
j,k
tend to evolve as homogenous clusters. Therefore one
ical average. We choose the local fields to point in
can approximate the inter-cell Hamiltonian in terms
opposite directions. For simplicity we set the largest
magnitude field h = −J (the largest magnitude of the total spin operators for the columns
2
field). We further choose J > h1 > 0. If h1 < J/2 8
then the global minimum corresponds to the qubits H1P,2 (cid:39)−nJS1z,2S2z,2 . (10)
inbothclusterspolarizedinthesamedirectionasthe
largest field, h . There also exists a “false” minimum WeobservethatthesystemHamiltoniancommutes
corresponding2to the clusters oriented in the opposite withthetotalspinoperatorsSk2,m =(cid:80)α∈x,y,z(Skα,m)2.
directions (each along its own local field). Given that all qubits in the initial state are polarized
We now explain the onset of the frustration in this alongthex-axisthisrestrictstheevolutiontothesub-
system. Weobservethatatthebeginningofquantum space of maximum total spin values n/4 for each col-
annealing umn. This subspace is spanned by the basis vectors
|n/4,m (cid:105) corresponding to the definite projections
k,m
(cid:104)σjz,k(cid:105)(cid:39)hkB(s)/A(s), B(s)/A(s)(cid:28)1. (7) of column spins on the z axis mk,m =−n/4,...,n/4.
As a measure of the error incurred by this approxi-
The Ising coupling terms in the problem Hamiltonian mation, it can be shown that the two lowest energy
(5) are quadratic while the local field terms are lin- levels are within 0.1% of the exact values for the case
ear in z-polarizations. Therefore at the beginning of of interest, n=8.
1
4
order in 1/n the Hamiltonian becomes
(cid:88) (cid:113)
HWKB(q ,q ,p ,p ,s)=−nA(s) 1−q2cosp
1 2 1 2 k k
k=1,2
(cid:88) (cid:16) (cid:88) (cid:17) n
+nB(s)J h q −nq2/4 − B(s)Jq q .
k k k 2 1 2
k=1,2 j
(13)
WKB theory based on this Hamiltonian describes
eigenstatesandeigenvalueswithlogarithmicaccuracy
in the asymptotic limit n(cid:29)1. It also gives a reason-
able estimates already for n=8 (see App. B).
We will now consider the potential corresponding
to a low energy description with very low momentum
U(q ,q ,s)=HWKB(q ,q ,0,0,s).
1 2 1 2
FIG. 3. Gap of the quantum Hamiltonian for h = 0.44
1 The same potential is obtained in Ref. [5] projecting
asafunctionoftheannealingparameter. Thecontinuous
the Hamiltonian of large spin operators Eqs. (9), (10)
lineistheenergydifferencebetweenthegroundstateand
over spin coherent states, which are product states.
the first excited state. The avoided crossing at s = 0.26
corresponds to a minimum gap of 180 MHz. The dashed The different panels in Figure 4 depict the potential
lineistheenergydifferencebetweenthegroundstateand U(q1,q2) for different values of the annealing param-
the second excited state. eter s with local field h = 0.44. Initially (s = 0)
1
there is only a global minimum at q = q = 0 corre-
1 2
sponding to all spins aligned with the x-direction. As
s grows the minima begins to move to the left corner
(−1,1) corresponding to the opposite orientations of
B. Effective energy potential the classical paths
the clusters. This effect was already mentioned above
of product states
in the general context. The terms in the effective po-
tential corresponding to the local fields h are linear
k
inq ,anddominatetheIsingcouplingenergybetween
We now derive the effective potential over product k
the large spins that is is quadratic in q for |q |(cid:28)1.
states used to study the difference between thermal k k
For larger values of s the Ising terms begin to com-
hoping among classical paths, and quantum tunnel-
pete with the local fields and the plateau is formed in
ing in the quantum models. This will also serve to
the vicinity of q = 0 following the local bifurcation
clarifythetunnelingpicturedescribedabove. Wefirst 1
ofU andgivingrisetoanewminimumcorresponding
transformtoarepresentationthatcontainsanexplicit
to the ferromagnetic alignment of the two clusters.
momentum operator. We think of each cluster k as
At some value s the two minima coexists with the
a “particle” with the coordinate proportional to its c
total spin z projection (cid:80) Sz . Because the x- same energy. For s > sc the minimum correspond-
j=1,2 k,j ingtotheferromagneticclusteralignmenthasalower
component of the total spin of the cluster does not
energy, smoothly connecting to the solution state at
commute with the z-component it is naturally associ-
the end of quantum annealing corresponding to the
ated with the momentum that causes the particle to
global minimum of the potential U(q ,q ,s = 1) at
move. Thisallowsustothinkofthez-componentofa 1 2
q =q =1.
largespinasaparticlemovinginaslowlytime-varying 1 1
In a closed system quantum adiabatic evolution al-
potential, formalizing the cartoon pictures sometimes
gorithm the system tunnels from the old global mini-
drawn to illustrate quantum annealing that show a
mum to the new one in the vicinity of s [5]. It tun-
particleescapingalocalminimuminacontinuouspo- c
nels under the barrier whose top approximately cor-
tential via tunneling.
responds to zero z-magnetization q =0. The tunnel-
1
ing corresponds to the avoided crossing between the
The canonically conjugated coordinate and mo-
two lowest energy levels of the quantum Hamiltonian
menta operators can be naturally introduced within
H (s) of Eq. (1) shown in Fig. 3. During the tun-
the WKB framework (see App. B) 0
neling the total spin of the left cluster switches its
direction. TheswitchingmanifestsitselfinFig.5asa
steepchangeins-dependenceofthequantummechan-
n
S1z,k+S2z,k = 2qk (11) ical average of the the left cluster polarization (cid:104)q1(cid:105) in
the instantaneous ground state. In contrast, the right
n(cid:113)
Sx +Sx ≈ 1−q2cosp . (12) cluster which does not tunnel displays a smooth be-
1,k 2,k 2 k k havior of its average polarization (cid:104)q (cid:105).
2
Inclassicaldynamicalevolutions,whentunnelingis
not possible, the system will continue to reside in the
where [q,p] = i(2/n) and 2/n (cid:28) 1 plays a role of initially global minimum emanating from the point
Planck constant in traditional WKB. To the leading (0,0). Fig. 6 shows the energies of the two classical
5
FIG. 4. The five snapshots in this figure show how the energy landscape evolves and a double-well potential is formed
duringannealingschedule. The3Dplotalsodepictssuchevolutionasafunctionofaneffectiveorientationangleforthe
largespins. Theminimumthatformsfirstwouldtrapaclassicalparticlemovinginthispotential. Laterintheannealing
evolution a second minimum forms and eventually becomes the global minimum. To reach this global minimum the
system state has to traverse the energy barrier between them. The origin of this bifurcation is explained in the text.
FIG. 6. Bifurcation for h = 0.44. The energies of the
1
paths of the double-well minima and the height of the ef-
fectiveenergybarrierinGHzareplottedinGHz. Wesub-
tract the instantaneous energy of the final false minima.
Thepathcorrespondingtotheglobalminimumappearsas
FIG. 5. Plot of (cid:104)q (cid:105) for the quantum ground state and
1
a bifurcation with higher energy. When this path crosses
firstexcitedstateasafunctionoftheannealingparameter
withthepathofthefinalfalseminimum,theheightofthe
s for h = 0.44. We also plot the value of the parameter
1
energy barrier is substantial.
q along the paths corresponding to the false and global
1
minima for the effective energy potential. The correspon-
dence is good, except at the avoided crossing where the
quantum states are entangled.
the value of the parameter q corresponding to the
1
local minima of U(q ,q ,s) as a function of the an-
1 2
nealing parameter s. The parameter q for the clas-
1
paths corresponding to the two minima of the energy sical path connecting to the final false minimum first
potential U(q ,q ,s) for h = 0.44. Classical dynam- alignswiththequantumgroundstateandlateraligns
1 2 1
ical evolutions will get trapped in the false minimum withthequantumfirstexcitedstate. Thispathstarts
path due to the bifurcation seen in this figure. Clas- at the point (0,0). The bifurcation path correspond-
sical trajectories can only reach the global minimum ing to the global minimum first appears with a defi-
throughthermalexcitations. InFigure5wealsoshow nite but small q , and later aligns with the quantum
1
6
ground state after the avoided crossing. wheren(ω)=[exp(βω)−1]−1 isthePlancknumberat
temperatureT =(cid:126)/k β andJ(ω)isthetemperature-
B
independentspectralfunction[20]thatcanbetreated
as an antisymmetric function of frequency, J(ω) =
III. MODELING THE ANNEALING
DYNAMICS [S(ω)−S(−ω)]/(cid:126)2.
The effect of the Ohmic noise on multiqubit quan-
tum annealing was studied numerically in [21]. This
A. Open quantum system models
workassumedweaksystem-environmentcouplingand
utilizedtheRedfieldformalismtoderivethequantum
Under realistic conditions, the performance of a Markovianmasterequationinthebasisofthe(instan-
quantum annealer as an optimizer can be strongly in- taneous) adiabatic eigenstates of the qubit Hamilto-
fluencedbythecouplingtotheenvironment. Inorder nian. It was built on earlier studies of open-system
to capture this effect we present a phenomenological quantum annealing in [22–26] where similar assump-
openquantumsystemmodelbyincorporatingtheex- tions were made.
perimentalcharacterizationofthenoisethatwasper- In addition to Ohmic noise, an important role
formed to date on D-Wave devices (the experimental is played by a low-frequency noise of the 1/f type
platform for the present investigation). [27,28]producedbythespinsintheamorphousparts
We shall assume that each flux qubit is coupled of the qubit device [12, 14, 29]. In current D-Wave
to its own environment with an independent noise chips this noise is coupled to the flux qubit relatively
source; this assumption is consistent with experimen- stronglyaswasshowninrecentexperiments[17]. Ad-
tal data [8]. We separate the bath excitations into ditionally, our analysis shows that noise effects are
two parts. Excitations with frequencies slower than significantly enhanced by collective effects associated
the annealing rate in our experiments (5 KHz) will with multiqubit cotunneling. While future genera-
be treated as a “static noise” whose effect can be in- tions of quantum annealer chips will hopefully have
cludedbyanappropriateaveragingofthesuccessrate reduced levels of flux qubit noise, it will still produce
over the local field errors (∼ 5% for the D-Wave II ahighlynonlineareffectforsufficientlyalargenumber
chip). The excitations with higher frequencies will be of cotunneling qubits.
modeled as a bath of harmonic oscillators. This ap- In recent years, the noise spectrum of flux qubits
proach is quite general, independent of the true phys- wasstudiedusingavarietyofapproachesthatincludes
ical source of the noise, being valid under the condi- dynamical decoupling schemes [30], free-induction
tion that the bath is in thermal equilibrium, can be Ramsey interference [31], coherent spectroscopy with
treatedwithintheconditionsoflinearresponsetheory, strong microwave driving [32], and macroscopic reso-
and has Gaussian fluctuations [19]. The correspond- nant tunneling (MRT) techniques [17, 33].
ing system-bath Hamiltonian is In MRT experiments the qubit state is probed in a
way that is most similar (compared with other meth-
1 (cid:88)2n ods) to the quantum annealing process itself, with-
H(s)=H (s)+ σzQ (s)+H , s=t/t .
0 2 µ µ B qa out using high-frequency driving and involving slow
µ=1 tunneling dynamics of each qubit within its group-
(14)
state manifold. While the exact microscopic models
Hereafter we will use for brevity a single index µ for
oflow-frequencynoisearenotwellunderstood[12,14]
single-qubit Pauli matrices instead of the double in-
its effect on system evolution in MRT experiments
dexation employed in Sec. II. In the equation above,
appeared to be well described by phenomenological
t is the duration of the quantum annealing process,
qa models [17, 33]. There, the quantity of interest is
H (t) is the Hamiltonian for an N-qubit system (as
0 the incoherent-tunneling rate between the “up” and
in Eq. 1), H is the standard Hamiltonian of the
B “down” eigenstates of the single flux qubit Hamilto-
bosonic bath, and Qµ(s) is a bosonic noise operator nian−1((cid:15)σz+∆σx)asafunctionofthebias(cid:15). In[34],
that couples the µth qubit to its environment. The 2
theGaussianformoftheMRTlineisdescribedwitha
coupling parameters of bosonic bath operator Q (s)
µ noise model whose spectral density is sharply peaked
dependontheannealingparametersthroughtheper-
at low frequency. In [17], this model is extended by
sistent current (see App. A). In what follows we use
attributingthelinearformofthetailsintheMRTline
the notation Q and make the dependence on s im-
µ shapetothehigh-frequency(Ohmic)partofthenoise
plicit for simplicity.
spectral density. The MRT data collected for small
Thepropertiesofthesystem’snoisearedetermined tunneling amplitudes (∆ < 1 MHz) and in a broad
by the noise spectral density S(ω): rangeoffrequencies(0.4MHz−4GHz)andtempera-
tures(21mK−38mK)aresurprisinglywelldescribed
(cid:90) ∞
eiωt(cid:104)eiHBtQ e−iHBtQ (cid:105)=S(ω)δ , (15) by a phenomenological “hybrid” noise model
µ ν µν
0
S(ω)=S (ω)+S (ω), (17)
LF HF
where the inclusion of the Kronecker delta function
wherethehigh-frequencypartofthenoisehasastan-
δ isaconsequenceoftheassumptionofindependent
µν
dard Ohmic form, [35]
baths. The temperature dependence of the spectral
density is given by J (ω)=ηωexp(ω/ω ), (18)
HF c
S(ω)=(cid:126)2J(ω)(n¯(ω)+1), (16) and the low-frequency part is described only by its
7
first two cumulants related to the width W and shift
(cid:15) of the MRT line:
p
(cid:90) ∞ dω (cid:126)ω
W2 = J (ω)coth , (19)
2π LF k T
0 B
(cid:90) ∞ dωJ (ω)
(cid:15) = LF (20)
p 2π ω
0
often called the reorganization energy, i.e. the energy
change of the bath degrees of freedom during an in-
coherent tunneling process. The width W and shift
(cid:15) measured in [17] satisfy the thermodynamic rela-
p
tion W2 = 2k(cid:126)BT(cid:15)p which requires that the character-
isticcutoffofthelow-frequencynoise, ω (cid:28)k T/(cid:126).
FL B
Thehybridmodelusesnoinformationaboutthenoise
spectrum at frequencies below W except the assump-
tion that (cid:126)ω (cid:28) W. It is justified because in the
LF
experiments W ∼K T [17].
B
The noise linewidth W = W(s), as probed by a
qubit, depends on the point at which it is taken the
annealling process via a persistent current. In our
study we will use the data from [17] collected on a FIG. 7. The minimum gap between the ground- and first
D-Wave I device. Similar measurements done in the excited-state levels ∆Emin = mins∈(0,1)[E1(s) − E0(s)]
(21) during quantum annealing as a function of the
D-Wave II chip used here yield similar parameters.
rescaled bias h /J is shown by the solid green line. The
Duringthequantumannealingthesystemtypically 1
horizontal boundary of thered-filled area at324MHz cor-
follows several stages identified in [21] depending on
responds to 15.5mK, which is the lowest temperature in
the magnitude of the instantaneous energy gap be-
our experiments. The bias value h /J = 0.5 corresponds
1
tween the ground-state energy E0(s) and the excited- to zero energy gap, which is achieved at the end of quan-
state energy E1(s) of the control Hamiltonian H0(s), tumannealingwhentheeigenstatesofH0(tqa)=HP with
whose energy spectrum evolves during quantum an- parallel and anti-parallel cluster orientations are degen-
nealing according to erate. The upper inset shows the avoided crossing be-
tween the energy levels E (s), E (s) in the weak-strong
1 0
H (s)|ψ (s)(cid:105)=E (s)|ψ (s)(cid:105). (21) cluster problem at h /J = 0.48. Dashed lines show the
0 γ γ γ 1
energy levels corresponding to the diabatic basis of states
Intheproblemstudiedherewithtwocoupledclusters, |φ∗(s)(cid:105),|φ∗(s)(cid:105) (29) formed by the rotation of the adia-
1 i
the separation between the ground state to the sec- batic eigenstates (21) that maximizes the average Ham-
ond excited state (> 3GHz) is much greater than the ming distance (25), (28) between the spin configurations.
temperature during the entire course of the quantum Thelowerinsetshowsthespinconfigurationsinbothclus-
annealing. Therefore the dynamics of the 2n qubit ters that dominate the characteristics of the eigenstates
beforeandaftertheavoidedcrossing. Duringthepassage
system can be described using the two lowest-energy
of the avoided crossing, spins in the left cluster (shown in
instantaneous eigenstates {|ψ (s)(cid:105),|ψ (s)(cid:105)}.
0 1 black) reverse their orientations.
At the beginning of quantum annealing, the qubit
Hamiltonian H (as in Eq. 1) is dominated by
0
the driver Hamiltonian term and the energy gap be-
tween the ground state and the first excited state
rate v . In our experiments, the quantum anneal-
((cid:39) 2A(s) ∼ 5GHz) is very large compared to the qa
ing rate is 50 kHz and the error due to Landau-Zener
temperature. At that stage the system resides in a
transitions is very small; it is 1% for the smallest
ground state |ψ (s)(cid:105) with overwhelming probability.
0 gap in our study of approximately 5MHz achieved
For local fields h < J/2, the system evolution goes
1 at h /J = 0.48, and it becomes completely neg-
through the so-called “avoided-crossing” region at in- 1
ligible for smaller values of h /J (e.g. 10−10 for
termediate times where E (s) and E (s) approach 1
1 0 h /J =0.46). Theminimumgapislessthanthetem-
closely to, and then repel from, each other. (See in- p1erature, ∆E (cid:46) k T, for most of the h /J values
set in Fig. 3.) This level repulsion occurs due to the min B 1
under study (cf. Fig. 7). Therefore the population of
cotunneling of qubits in the left cluster between the
the excited state at the end of quantum annealing is
opposite z-polarizations. At the point where the gap
entirely due to thermal excitation [15, 16, 36].
E (s)−E (s) reaches its minimum ∆E the cor-
1 0 min
responding adiabatic eigenstates are formed by the When the instantaneous energy gap E (s)−E (s)
1 0
symmetric and anti-symmetric superpositions of the is sufficiently large compare to W(s), the coupling to
cluster orientations (cf. Fig. 7). environment can be treated as a perturbation. This
The success of closed-system quantum annealing is regimecorrespondstotheinitialstageofevolutionup
limited by the Landau-Zener transitions away from to approximately the avoided crossing. In the per-
the adiabatic group state whose probability decreases turbative regime, the transition rate between the first
exponentially with the ratio ∆E2 /v of the square excitedstateandagroundstateofthecontrolHamil-
min qa
of minimum energy gap to the quantum annealing tonian H (s) is given by Fermi’s golden rule for a
0
8
FIG. 8. The solid black line shows the dependence of the
FIG. 9. Dependence of the overlap factor b on s for
effectivelinewidthofthelow-frequencynoiseh1/2(s)W(s)
h = 0.44J. The two plots correspond to the overlap
on annealing parameter s. The dashed blue line is a plot 1
coefficient in the adiabatic basis {|ψ (cid:105),|ψ (cid:105)} and rotated
of the transition rate (57) ΓFGR vs. s calculated in sec- 0 1
10 (diabatic) basis (29). Vertical dashed line indicates the
ond order in spin-boson interaction using Fermi’s golden
pointwheretheminimumenergygapisreached(avoided-
rule (22). The solid blue is a plot of the transition rate
crossing). Red,blue,andgraycolorsindicatethedifferent
Γ calculated using NIBA (57). The green line shows
10 stages of the quantum annealing described in the text.
thes-dependenceoftheenergygapE (s)−E (s)between
1 0
the two lowest-energy levels during the annealing at 15.5
mK. The lowest temperature in our experiments is shown
crossing correspond to opposite total spin z- projec-
by the horizontal red dotted line. All plots correspond to
tionsoftheleftcluster. Whensincreases,theaverage
biash =0.44J. Attheearlystageofquantumannealing
1
h1/2(s)W(s)(cid:28)E (s)−E (s)andΓFGR givesanaccurate “Hamming distance” between the eigenstates,
1 0 10
description of the dynamics. The inset show the transi-
2n
tion rates ΓF10GR, Γ10 vs. s using the same units as in the h(s)= 1 (cid:88)|Z11(s)−Z00(s)|2, (25)
main plot but showing a greater range of values. Vertical 4 µ µ
white lines in the inset and main plot mark the value of µ=1
s =0.269 where the deviation of the system state from
eq also increases very steeply as shown in Fig. 10. (We
instantaneous Gibbs distribution is 1%. This deviation
note that the maximum value of h(s) is proportional
grows very rapidly for greater values of s due to the sup-
ton.) Inthatregionthetransitionbetweentheeigen-
pressionofthetransitionrates(cf. Fig.9). Onlythestage
ofthequantumannealingwiths(cid:38)s determinesthefinal states requires qubit cotunneling of a progressively
eq
ground state population. One can see that at that stage higher order, leading to an exponential decay of the
h1/2(s)W(s) (cid:29) Γ , justifying the NIBA approximation overlap coefficient b(s) with s and a steep decelera-
10
discussed in the main text. tion of the environment-induced transitions between
thetwostates,ascanbeseenfromtheplotofΓFGR(s)
10
in Fig. 8.
single-boson process:
ΓFGR(s)=b(s)J(ω )[n¯(ω )+1], (22)
10 10 10
where ω =(E (s)−E (s))/(cid:126) is the Bohr frequency
10 1 0
for the transition and
2n
1 (cid:88)
b(s)= |Z10(s)|2. (23)
4 µ
µ=1
Here and below we use the following notation for the
matrix elements:
Zµγγ(cid:48)(s)=(cid:104)ψγ(s)|σµz|ψγ(cid:48)(s)(cid:105), γ,γ(cid:48) =0,1. (24)
Thisisanoverlapfactorthatdetermineshowstrongly
the transition 1 ↔ 0 is coupled to the environment.
Fig. 9 shows the dependence of b(s) on the annealing
FIG. 10. Dependence of the average Hamming distance
parameter s. We observe a steep exponential fall-off
between the two lowest-energy eigenstates on s for adi-
off this coefficient after the avoided crossing. This abatic basis (|ψ (cid:105),|ψ (cid:105)) and rotated diabatic basis (29)
0 1
happens because, starting from the avoided cross- for h = 0.44J. The vertical dashed line indicates the
1
ing region, intra-cell ferromagnetic interaction plays pointwheretheminimumenergygapisreached(avoided-
a substantial role by causing the spins in each unit crossing). Red,blue,andgraycolorsindicatethedifferent
celltomoveinunison,formingtwoclusterswithtotal stages of the quantum annealing described in the text.
spin value n/2 each. As can be seen from Fig. 5, the
first two eigenstates |ψ (cid:105) and |ψ (cid:105) after the avoided The fact that qubits within each unit cell tend to
0 1
9
move together, forming large spins, amplifies the ef-
fect of the environment on their quantum dynam-
ics. In particular, we will show below that the ef-
fective linewidth of the low-frequency noise as seen
by the two-state system {|ψ (s)(cid:105),|ψ (s)(cid:105)} system is
1 0
h1/2(s)W(s)andthattheeffectiveOhmiccoefficientis
h(s)η(s). Thisamplificationbecomesimportantwhen
clustersincreasetheirz-polarizationsandh∼n>>1
(see Fig. 8) at later stages of the quantum anneal-
ing. For sufficiently large h1/2(s)W(s) (cid:38) ∆E ,
min
the description of the system dynamics becomes sub-
stantially non-perturbative in the spin-boson inter-
action. Equilibria of the environmental degrees of
freedom shift depending on the collective qubit-state
FIG. 11. Optimal rotation angles ϑ∗ vs. s for diabatic
evolution, which in turn affects the state itself caus-
eigenstates(29)atdifferentvaluesofh /J =0.35+0.01i,
ing the polaronic effect. In this case, the adiabatic 1
with i = 1,2,....,13 corresponding to the curves in the
basis of instantaneous eigenstates {|ψ (s)(cid:105),|ψ (s)(cid:105)}
1 0 figure numbered from the bottom to the top. Red corre-
formed by the superposition of up and down clus- sponds to the thermalization phase, blue corresponds to
ter orientations loses its physical significance. In- the loss of thermalization, and black corresponds to the
stead, the dynamics occurs between the two diabatic frozenphasedescribedinthemaintext. Optimalrotation
states {|φ (s)(cid:105),|φ (s)(cid:105)} with the predominately op- anglesonlyhavephysicalmeaningstartingfromaboutthe
1 0
posite cluster orientations corresponding (roughly) to end of the thermalization phase where the spins in the
the bottoms of the wells of the classical potential in cells start to behave in a concerted manner. All angles
approach π/2 in a frozen phase corresponding to the adi-
Fig. 4 separated by the barrier.
abatic eigenstates.
Weintroduceaunitaryrotationonangleϑdefining
a new basis of states |φ (ϑ,s)(cid:105):
0,1
(cid:18) (cid:19) Wewillshowbelowthatthenon-perturbativetreat-
(cid:88) ϑ π
|φ (cid:105)= (−1)ij+j+1 cos −(−1)i+j |ψ (cid:105), ment of the effects of noise and dissipation does not
i j
2 4
change the Markovian nature of the system dynamics
j=0,1
(26) but modifies the instantaneous transition rate Γ10(s)
where i,j =0,1, and corresponding matrix elements compared to its value ΓFGR(s) (22) from a single-
10
boson process. The non-perturbative analysis in the
Zµγγ(cid:48)(ϑ,s)=(cid:104)φγ(ϑ,s)|σµz|φγ(cid:48)(ϑ,s)(cid:105), γ,γ(cid:48) =0,1. rotatedbasisisjustifiedwithinthecontextofthethe-
(27) ory of spin-boson interaction developed in [35]. The
We find the angle ϑ (s) that maximizes the average individual transitions due to the coupling to bosons
∗
Hamming distance between the states areassociatedwithso-called“blip-cojourn”pairswith
blips forming a dilute gas if the duration of the blip
1 (cid:88)2n τb =1/h1/2(s)W(s) is much shorter than the charac-
h∗(s)=mϑax4 |Zµ11(ϑ,s)−Zµ00(ϑ,s)|2. (28) teristic inter-blip distance ∼ 1/Γ10(s), with
µ=1
Γ (s)(cid:28)h1/2(s)W(s). (30)
10
Then our new instantaneous basis will be
Intheadiabaticbasis,h(s)canreachverysmallvalues
|φ∗(s)(cid:105)=|φ (ϑ (s),s)(cid:105), γ =0,1. (29) neartheavoidedcrossing(seeFig.10)duetothespuri-
γ i ∗
ousquantuminterferenceeffects. Intheoptimallyro-
The system dynamics is mostly a hopping between tated(diabatic)basis,τ−1 =(h∗(s))1/2W(s)ismono-
b
these states associated with the incoherent cotunnel- tonically increasing during the annealing to its max-
ing of the spins in the left cluster connecting states imum value (n=8 in the problem of interest). When
with predominantly “up” (|φ∗(cid:105)) and “down” (|ketφ∗) the condition (30) is satisfied at the last stages of the
1 0
qubit configurations that are in average a Hamming annealing, the gas of blips is dilute and we can apply
distance h (s) apart. (h (s) is shown in Fig.10.) In the Noninteracting-blip Approximation (NIBA) [35].
∗ ∗
essence, this approach is related to the pointer basis While this method was applied to analyze a Landau-
idea, that the system tends to be localized in states Zener problem in a driven spin-1/2 system coupled
induced by environmental coupling [37]. to a finite temperature bath in a number of papers
The rotation angle ϑ∗(s) during the annealing is [15, 16, 38], its application in the present context of
shown in Fig.11 for different vales of h /J. In the strongly correlated qubit dynamics is novel and leads
1
later stages of the annealing the angle ϑ always ap- to qualitatively different features.
proaches the value π/2 that corresponds to the adia- Before we proceed further with NIBA analysis, we
batic (energy) basis with state |φ∗(cid:105) being an excited emphasize that our theory will not provide an accu-
1
state 26. In other words, quantum annealing along ratetreatmentoftheregionofswhereh1/2(s)W(s)∼
the pointer states arrives at the encoded solution of Γ , representing a crossover between the pertur-
10
the computational problem because toward the end bative treatment based on single-boson processes
of the annealing the pointer basis converges back to and Ohmic spectral density (ΓFGR(s)) and non-
10
eigenstates of the problem Hamiltonian. perturbativeNIBAtheorythatincludeslow-frequency
10
noise (Γ (s)). In this intermediate region the low- exp(−(E −E )/k T). Whenthetemperaturegrows,
10 1 0 B
frequency noise cannot be described by the two low- thetransitionratesdominatedbylow-frequencynoise
est cumulants of its special density and we are simply changeverylittle,whiletheGibbsfactorincreasesap-
luckingtheexperimentaldataforthetheoreticalanal- preciably as can be seen from Fig. (12). This in turn
ysis. increasesthepopulationoftheexcitedstate,reducing
Our main observation is that this region does not the success of the quantum annealing as will be seen
affectthepopulationofthegroundstateattheendof later.
thequantumannealing. Sincethequantumannealing
rate 1/t =50 kHz is constant, the system stays very
qa
close to an instantaneous thermal equilibrium while B. Non Interacting Blip Approximation
ΓF10GR(s), Γ10 (cid:29)1/tqa (31) To implement NIBA in our problem we need to ex-
plicitlyrepresenttheenvironmentaldegreesofthebo-
(see Fig. 8). This is a “thermalization phase” of the
son operators that diagonalize the free-boson Hamil-
annealing process. On the other hand, due to the tonian H =(cid:80)2n (cid:80) (cid:126)ω (b† b +1/2):
strongeffectoflow-frequencynoiseonD-Wavequbits B µ=1 ν µν µν µν
the condition Q =(cid:88)λ (b† +b ), µ=1,...,2n. (33)
µ µν µν µν
1/t (cid:28)h1/2(s)W(s) (32) ν
qa
Here λ and ω are the microscopic parameters
µν µν
is held almost everywhere except for the very early that will not enter into any observable directly but
stages (cf. Fig. 8). Therefore the non-perturbative will do so only via the spectral function (16) J(ω) =
regime (30) is established well within the thermaliza- 2π(cid:80) (λ /(cid:126)ω )2δ(ω−ω ),identicalforeachqubit.
µ µν µν µν
tion phase where the (Gibbs) distribution function is Weproceedbymakingasmallpolarontransformation
not sensitive to the noise model. of the original Hamiltonian with unitary operator
(cid:34) (cid:35)
(cid:88)
U(s)=exp −Λ (s)|φ∗(s)(cid:105)(cid:104)φ∗(s)|
γ γ γ
γ
. Here we use the notation
Λ (s)=(cid:88)ξγγ(b −b† ), ξγγ(cid:48) = λµν Zγγ(cid:48), (34)
γ µν µν µν µν (cid:126)ω µ
µν
µν
wheretheindexesµ,ν andγ runovertheirrespective
ranges defined above and matrix elements Zγγ(cid:48) are
µ
defined in the diabatic basis as
Zγγ(cid:48)(s)=(cid:104)φ∗(s)|σz|φ∗ (s)(cid:105). (35)
µ γ µ γ(cid:48)
The system Hamiltonian after the transformation
H = UH U−1 can be written (up to identity opera-
0 0
FIG. 12. In main plot, the solid line with red, blue, and tor) as
blackcoloringshowsthedependenceoft Γ (57)onthe
qa 10 H =(cid:126)Ω(|φ∗(cid:105)(cid:104)φ∗|−|φ∗(cid:105)(cid:104)φ∗|), (36)
annealing parameter s. Red corresponds to the thermal- 0 1 1 0 0
ization phase, blue corresponds to the loss of thermaliza-
where
tion, and black corresponds to the frozen phase described
in the main text. The solid green line is a plot of the Ω(s)=ε(s)sinϑ (s)−(cid:126)(cid:15) d(s) (37)
∗ p
Gibbsfactorfortheinstantaneousenergydifference(cid:126)Ω(s)
and ε = (E − E )/(cid:126). The first term above cor-
between the diabatic states given in (37). Both lines are 1 0
plottedat15.5mK.Dashedlinescorrespondtothetemper- responds to the system energy gap in a diabatic
ature 35mK. The inset shows the evolution of the instan- (pointer) basis and the second term gives the pola-
taneouspopulationofthediabaticeigenstate(cid:104)φ∗(s)|ρ|φ∗(cid:105). ronic shift due to reorganization energy of the envi-
0 0
The colors correspond to the same stages of quantum an- ronment. Similar to the linewidth, it is renormalized
nealing as in the main plot. All data are for the value of with respect to its single qubit MRT value by a coef-
the local field h1/J =0.44. ficient d(s) reflective of a collective qubit behavior:
2n
At a later stage the thermal distribution can no d= 1 (cid:88)(cid:2)(Z11)2−(Z00)2(cid:3). (38)
longer be supported due to the steep slowdown of the 4 µ µ
µ=1
transition rates as shown in Fig. (12). Then the sys-
The Hamiltonian of the coupling to environment
tem enters the final “frozen” stage where the transi-
afterthepolarontransformH =H10 |φ∗(cid:105)(cid:104)φ∗|+h.c.
tions are suppressed over the period of the annealing int int 1 0
is strictly non-diagonal in the basis states and
and part of the system population remains trapped
ionf Fthige.e1x2c,ittehdesstuactcee.ssApsrocbaanbiblietyseisendeitnertmheineindsebryt Hγinγt(cid:48) =(cid:32)(cid:126)2∆ +(cid:88)Zµγγ(cid:48)λµν(b†µν +bµν −2ξµγνγ(cid:48))(cid:33) eΛγ(cid:48)γ,
the occupation of the ground state at the beginning µ,ν
of the frozen phase determined by the Gibbs factor (39)