Close
About
FAQ
Home
Collections
Login
USC Login
Register
0
Selected
Invert selection
Deselect all
Deselect all
Click here to refresh results
Click here to refresh results
USC
/
Digital Library
/
University of Southern California Dissertations and Theses
/
Open-system modeling of quantum annealing: theory and applications
(USC Thesis Other)
Open-system modeling of quantum annealing: theory and applications
PDF
Download
Share
Open document
Flip pages
Contact Us
Contact Us
Copy asset link
Request this asset
Transcript (if available)
Content
Open-system modeling of quantum annealing: theory and applications by Ka Wa Yip A Dissertation Presented to the FACULTY OF THE GRADUATE SCHOOL UNIVERSITY OF SOUTHERN CALIFORNIA In Partial Fulllment of the Requirements for the Degree DOCTOR OF PHILOSOPHY (Physics) May 2021 Copyright 2021 Ka Wa Yip Acknowledgements I would like to thank my advisor Prof. Daniel Lidar for his supports, teachings and guidance throughout my Ph.D. This dissertation would be impossible without him, and his wisdom and leadership. Prof. Lidar is a great mentor and fatherly gure. He taught me a lot and I am honored to be his student. Next I would like to thank Dr. Tameem Albash for his guidance and discussions especially regarding all the details in simulations and physics, at the early stage of my doctoral studies. Finally, without naming anyone, I would like to thank all the friends, sta, and professors at USC. The sta at the department provided excellent administrative assistance for my doctoral studies. ii Table of Contents Acknowledgements ii List of Tables viii List of Figures ix Abstract xviii Chapter 1: Introduction 1 Chapter 2: Quantum trajectories for time-dependent adiabatic master equa- tions 6 2.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6 2.2 Adiabatic master equation in Lindblad form . . . . . . . . . . . . . . . . . . 9 2.3 Stochastic Schr odinger Equation . . . . . . . . . . . . . . . . . . . . . . . . . 10 2.3.1 Unravelling the master equation . . . . . . . . . . . . . . . . . . . . . 11 2.3.2 Deterministic evolution and jump process . . . . . . . . . . . . . . . . 13 2.4 Simulation procedure for adiabatic quantum trajectories . . . . . . . . . . . 16 2.5 V: Case studies . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18 2.5.1 8-qubit chain . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21 2.5.2 8-qubit non-adiabatic example . . . . . . . . . . . . . . . . . . . . . . 23 2.5.3 16-qubit \tunneling-probe" Hamiltonian . . . . . . . . . . . . . . . . 26 2.6 Simulation cost comparison . . . . . . . . . . . . . . . . . . . . . . . . . . . 27 2.6.1 Number of trajectories . . . . . . . . . . . . . . . . . . . . . . . . . . 28 2.6.2 Cost comparison . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 28 2.6.3 Parallel implementation . . . . . . . . . . . . . . . . . . . . . . . . . 31 2.7 Conclusion of this chapter . . . . . . . . . . . . . . . . . . . . . . . . . . . . 32 2.8 Acknowledgments of this chapter . . . . . . . . . . . . . . . . . . . . . . . . 33 Chapter 3: 1=f noise 34 3.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 34 3.2 1=f noise: Bloch vector-Magnetic eld approach . . . . . . . . . . . . . . . . 35 3.2.1 One qubit and one uctuator . . . . . . . . . . . . . . . . . . . . . . 35 3.2.1.1 Qubit . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 36 3.2.1.2 Fluctuator . . . . . . . . . . . . . . . . . . . . . . . . . . . . 36 iii 3.2.1.3 Comparison of analytical solutions and stochastic simulations 37 3.2.2 One qubit and an ensemble of uctuators . . . . . . . . . . . . . . . . 39 3.2.2.1 An ensemble of uctuators to generate an overall 1=! power spectrum: S(!) =A=!. . . . . . . . . . . . . . . . . . . . . 39 3.2.2.2 Distributions of the uctuators . . . . . . . . . . . . . . . . 41 3.2.2.3 Simulation results: The importance of the number of decades and min . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 41 3.2.2.4 Simulation results: Dependence on b=hbi . . . . . . . . . . 42 3.2.2.5 Simulation results: Non-exponential decay . . . . . . . . . . 43 3.3 1=f noise: Stochastic Hamiltonian formulation . . . . . . . . . . . . . . . . . 44 3.3.1 Single uctuator . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 45 3.3.2 Many uctuators . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 45 3.3.3 Case studies . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 46 3.3.3.1 Single uctuator: A(t) = 0 and B(t) = 0 for all t . . . . . . 46 3.3.3.2 Multiple uctuators with a time-dependent system Hamiltonian 47 3.4 Eect of temperature, T 1 , T 2 decay time, and strength of uctuators . . . . . 49 3.4.1 Two temperature limits . . . . . . . . . . . . . . . . . . . . . . . . . 49 3.4.1.1 Innite temperature T!1 . . . . . . . . . . . . . . . . . . 50 3.4.1.2 Zero temperature T! 0 . . . . . . . . . . . . . . . . . . . . 50 3.4.2 T 1 and T 2 time . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 51 3.4.2.1 T 1 time . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 51 3.4.2.2 T 2 time . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 53 3.4.3 Dierent uctuator strengths . . . . . . . . . . . . . . . . . . . . . . 53 3.5 Modeling 1=f noise for superconducting qubit: Capacitively shunted ux qubit (CSFQ) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 55 3.5.1 Capacitively shunted ux qubit (CSFQ) circuit Hamiltonian . . . . . 56 3.5.2 CSFQ: Persistent current operator as uctuator . . . . . . . . . . . . 56 3.5.3 Case studies . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 58 3.5.4 Closed-system annealing evolution . . . . . . . . . . . . . . . . . . . . 59 3.5.5 Adiabatic master equation simulation . . . . . . . . . . . . . . . . . . 59 3.5.6 1=f noise simulation . . . . . . . . . . . . . . . . . . . . . . . . . . . 59 3.6 Modeling 1=f noise for superconducting qubit: Transmon . . . . . . . . . . . 61 3.6.1 Transmon Circuit Hamiltonian . . . . . . . . . . . . . . . . . . . . . 61 3.6.2 Transmon: Form of uctuators . . . . . . . . . . . . . . . . . . . . . 61 3.6.3 Case studies . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 62 3.6.4 Adiabatic master equation simulation . . . . . . . . . . . . . . . . . . 63 3.6.5 1=f noise simulation . . . . . . . . . . . . . . . . . . . . . . . . . . . 64 3.7 Fitting experimental data and reverse engineering of the spectrum . . . . . . 65 3.7.1 Fitting results with IBM machine . . . . . . . . . . . . . . . . . . . . 66 3.7.1.1 Free induction decay: Fluctuation A =I p . . . . . . . . . . 67 3.7.1.2 Dynamical decouplings: Fluctuation A = [I p ] 22 . . . . . . . 67 3.8 Conclusion of this chapter . . . . . . . . . . . . . . . . . . . . . . . . . . . . 70 3.9 Acknowledgments of this chapter . . . . . . . . . . . . . . . . . . . . . . . . 70 iv Chapter 4: Quantum feedback control 71 4.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 71 4.1.1 Feedback and control . . . . . . . . . . . . . . . . . . . . . . . . . . . 72 4.1.2 Weak measurement . . . . . . . . . . . . . . . . . . . . . . . . . . . . 72 4.2 Thermal excitation as error: detection . . . . . . . . . . . . . . . . . . . . . 73 4.3 Delay between detection and feedback . . . . . . . . . . . . . . . . . . . . . 74 4.4 Feedback superoperatorL f : Lindbladian . . . . . . . . . . . . . . . . . . . . 75 4.5 Case studies . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 76 4.6 Feedback superoperatorL f : Hamiltonian and basis of feedback . . . . . . . . 77 Chapter 5: Reverse Annealing: Iterative reverse annealing 81 5.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 82 5.2 Iterative reverse annealing: the protocol . . . . . . . . . . . . . . . . . . . . 82 5.3 Ferromagnetic p-spin model . . . . . . . . . . . . . . . . . . . . . . . . . . . 85 5.4 Unitary dynamics . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 88 5.5 Open system dynamics subject to dephasing-induced relaxation . . . . . . . 91 5.6 Open system dynamics with a pause . . . . . . . . . . . . . . . . . . . . . . 97 5.7 Conclusion of this chapter . . . . . . . . . . . . . . . . . . . . . . . . . . . . 98 5.8 Acknowledgments of this chapter . . . . . . . . . . . . . . . . . . . . . . . . 99 Chapter 6: Experimental iterative reverse annealing 104 6.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 104 6.2 Problem denition and protocols of experiment . . . . . . . . . . . . . . . . 106 6.2.1 p-spin model with p = 2 . . . . . . . . . . . . . . . . . . . . . . . . . 107 6.2.2 Reverse annealing protocol . . . . . . . . . . . . . . . . . . . . . . . . 108 6.3 Experimental results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 110 6.3.1 Dependence on initial conditions . . . . . . . . . . . . . . . . . . . . 110 6.3.1.1 m 0 = 0 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 110 6.3.1.2 m 0 = 0:8; 0:9: up/down symmetry breaking for 0:3.s inv . 0:5112 6.3.1.3 s inv & 0:5>s : freezing . . . . . . . . . . . . . . . . . . . . 112 6.3.1.4 s inv . 0:3<s : spin-bath polarization . . . . . . . . . . . . 113 6.3.1.5 Total success probability . . . . . . . . . . . . . . . . . . . . 113 6.3.2 Dependence on annealing time . . . . . . . . . . . . . . . . . . . . . . 114 6.3.3 Eects of pausing . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 115 6.3.4 Size dependence . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 116 6.3.5 Eects of iteration . . . . . . . . . . . . . . . . . . . . . . . . . . . . 117 6.4 Numerical simulations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 119 6.4.1 Closed system model . . . . . . . . . . . . . . . . . . . . . . . . . . . 119 6.4.1.1 Dependence on system size and annealing time . . . . . . . 120 6.4.1.2 Dependence on initial state and number of iterations . . . . 121 6.4.2 Open quantum system model . . . . . . . . . . . . . . . . . . . . . . 123 6.4.2.1 Dependence on annealing time . . . . . . . . . . . . . . . . 126 6.4.2.2 Dependence on initial condition and the number of cycles . 128 6.4.2.3 Size dependence and partial success probability . . . . . . . 129 6.5 Semi-classical simulations . . . . . . . . . . . . . . . . . . . . . . . . . . . . 130 v 6.5.1 Total success probability . . . . . . . . . . . . . . . . . . . . . . . . . 131 6.5.2 Partial success probability . . . . . . . . . . . . . . . . . . . . . . . . 132 6.6 Conclusion of this chapter . . . . . . . . . . . . . . . . . . . . . . . . . . . . 134 6.7 Acknowledgments of this chapter . . . . . . . . . . . . . . . . . . . . . . . . 136 Chapter 7: Adiabatic reverse annealing 137 7.1 Introduction and terminology . . . . . . . . . . . . . . . . . . . . . . . . . . 137 7.2 Preliminary I: Adiabatic reverse annealing (ARA) Hamiltonian spectrum . . 138 7.2.1 Controlling the parameters: . . . . . . . . . . . . . . . . . . . . . . 139 7.2.2 Controlling the parameters: Initial states i or (H init ) . . . . . . . . . 140 7.2.3 Controlling the parameters: N . . . . . . . . . . . . . . . . . . . . . . 141 7.3 Preliminary II: Time evolution . . . . . . . . . . . . . . . . . . . . . . . . . . 141 7.3.1 Closed system . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 142 7.3.2 Open system: Lindbladian dynamics . . . . . . . . . . . . . . . . . . 142 7.4 Preliminary III: Success metrics and classes of problems . . . . . . . . . . . . 143 7.4.1 Success probability . . . . . . . . . . . . . . . . . . . . . . . . . . . . 143 7.4.2 Time to solution (TTS) . . . . . . . . . . . . . . . . . . . . . . . . . 143 7.4.3 Classes of problems . . . . . . . . . . . . . . . . . . . . . . . . . . . . 143 7.5 Illustrative examples . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 144 7.5.1 Results of 8-qubit simulation: = 40ns. = 1. . . . . . . . . . . . . 144 7.5.2 Results of 4-qubit simulation: = 4ns. = 1. . . . . . . . . . . . . . 145 7.5.3 Results of 4-qubit simulation: = 400ns. = 0:3. . . . . . . . . . . . 146 7.6 ARA: Random initial state c = 0:5. . . . . . . . . . . . . . . . . . . . . . . . 146 7.6.1 QA . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 149 7.7 ARA performance in terms of and H init . . . . . . . . . . . . . . . . . . . . 149 7.7.1 Min. Gap vs H init vs . N = 4. . . . . . . . . . . . . . . . . . . . . . 150 7.7.2 Dierent Anneal times + Dierent =f1; 2; 5g. N = 4. ARA Open . . 150 7.7.3 Dierent Anneal times + Dierent < 1. N = 4. ARA Closed and ARA Open . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 152 7.7.4 Dierent Anneal times. = 0:1, N # = 2. N =f8; 10g. ARA Open . . . . 153 7.8 Independent vs collective SB coupling . . . . . . . . . . . . . . . . . . . . . . 154 7.8.1 Long anneal times . . . . . . . . . . . . . . . . . . . . . . . . . . . 154 7.8.2 Intermediate anneal times . . . . . . . . . . . . . . . . . . . . . . . 155 7.9 Conclusion of this chapter . . . . . . . . . . . . . . . . . . . . . . . . . . . . 156 Chapter 8: Conclusion 158 References 161 Appendices 166 A Quantum trajectories for time-dependent adiabatic master equations . . . . 167 A.1 Error estimates . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 167 A.1.1 Error associated with the norm squared of a no-jump trajectory167 A.1.2 Error associated with probability elements . . . . . . . . . . 171 vi A.2 Proof of equivalence between the master equation and trajectories for- mulations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 173 A.2.1 To jump or not to jump . . . . . . . . . . . . . . . . . . . . 173 A.2.2 Averaging over trajectories . . . . . . . . . . . . . . . . . . 174 A.2.3 Upper bound on t . . . . . . . . . . . . . . . . . . . . . . 176 A.3 On the validity of waiting times (quantum time-dependent operators) 176 A.4 Eect of the annealing schedule on the 8-qubit chain example . . . . 178 A.5 Proof that eigenstates of H S are modied only under jumps . . . . . 179 A.5.1 Whenj (t)i is an eigenstate of H S (t) . . . . . . . . . . . . . 180 A.5.2 No degeneracy in " b (t) . . . . . . . . . . . . . . . . . . . . . 181 A.6 Derivation of Eq. (2.31) . . . . . . . . . . . . . . . . . . . . . . . . . 182 B Experimental iterative reverse annealing . . . . . . . . . . . . . . . . . . . . 183 B.1 Spectrum of thep-spin problem withp = 2 and scaling of the minimum gap with N . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 183 B.2 Bound on success probabilities, for closed system and open-system dynamics with collective system-bath interaction . . . . . . . . . . . . 183 B.3 SVMC and SVMC-TF for reverse annealing . . . . . . . . . . . . . . 187 vii List of Tables 3.1 Form of uctuation (Qubit vs CSFQ) . . . . . . . . . . . . . . . . . . . . . . 57 3.2 Form of uctuation. z =j0ih0jj1ih1j. ; I p = I p ( x ) is the persistent current operator. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 62 viii List of Figures 1.1 D-Wave 2000Q NASA annealing schedules. . . . . . . . . . . . . . . . . . . . 2 2.1 A depiction of the stochastic evolution of the statej i by an innitesimal time state at time t. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14 2.2 Graphs of (a) the 8-qubit chain, (b) the 8-qubit Hamiltonian exhibiting a small gap, and (c) the 16-qubit \tunneling-probe" Hamiltonian of Ref. [47]. (a) Only the rst qubit is subjected to an applied eld and each qubit is fer- romagnetically coupled withJ =1. (b) Solid lines corresponds to ferromag- netic coupling and dashed lines corresponds to antiferromagnetic coupling. The thickness denotes the strength of the coupling. Local elds are shown inside the circles. Full parameters are given in Eq. (2.32). (c) The left 8-qubit cell and right 8-qubit cell are each subjected to applied elds with opposite direction. Each qubit is ferromagnetically coupled to others as shown by the lines, with J =1. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 19 2.3 The evolution of the population in the instantaneous ground state for the 8-qubit problem in Eqs. (2.27) and (2.30) for a total time of t f = 10s and temperature 2:6GHz (in ~ 1 units), as a function of the normalized time s = t=t f . The quantum trajectories results with 10 4 trajectories (QT, blue circles) are in excellent agreement with the adiabatic master equation (ME, red solid line). Inset: the convergence of the ground state population (averaged over quantum trajectories) towards the master equation result as a function of the number of trajectories, ats = 1. The error bars represent 2 condence intervals, where is the standard deviation of the mean generated by taking 10 3 bootstraps over the number of trajectories. . . . . . . . . . . . . . . . . 20 ix 2.4 The overlap squared of the (normalized) state with the instantaneous ground state ofH S (t) for a typical single trajectory of the 8 qubit chain in Sec. 2.5.1, with t f = 10s and temperature 2:62GHz, as a function of the normalized times =t=t f . The sudden changes in the overlap are due to the action of the jump operatorsfA i (t)g, taking the state from one eigenstate to another. This is to be contrasted with the smooth behavior of Fig. 2.3 when we average over dierent trajectories. Inset: a histogram of the net number of jumps to the instantaneous ground state (GS). A negative number indicates a jump out of the ground state, and a positive number indicates a jump towards the ground state. The change from negative to positive net jumps occurs at the minimum gap point. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 22 2.5 The evolution of the population in the instantaneous ground state for the 8-qubit problem in Eqs. (2.27) and (2.32) for a total time of t f = 10s and temperature 1:57GHz, as a function of the normalized time s = t=t f . The quantum trajectories results with 5 10 3 trajectories (blue circles) are in excellent agreement with the master equation (red solid line). We also show the closed-system evolution (yellow dashed line) to highlight that the evolution is not adiabatic. Inset: the instantaneous energy gap between the rst excited state and the ground state during the anneal. The minimum occurs at s = 0:46, coinciding with the sharp discontinuity observed in the instantaneous ground state population. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 23 2.6 The overlap squared of the (normalized) state with the instantaneous ground state of H S (t) (blue solid curve) and the sum of the rst three instantaneous excited states (red dashed curve) for a typical single trajectory of the 8-qubit problem in Sec. 2.5.2, witht f = 10s and temperature 1:57GHz, as a function of the normalized time s = t=t f . A diabatic transition occurs at s = 0:46. The continuous decay immediately afterward and between the `incomplete' jumps is to be contrasted with the step-function like single trajectories of the adiabatic case seen in Fig. 2.4. . . . . . . . . . . . . . . . . . . . . . . . . . 25 2.7 Quantum trajectory results for a temperature of 12mK and 20mK using the DW2X annealing schedule. Results are averaged over 5000 trajectories. A revival of instantaneous ground state population occurs after the minimum gap (shown by the dashed vertical line). . . . . . . . . . . . . . . . . . . . . 26 3.1 Bloch vector representing the state of a qubit in the rotating (with angular fre- quencyB 0 ) frame of reference. In the laboratory frame of reference it precesses around the z-axis with (time-dependent) angular velocity B. Fluctuations of this velocity is b(t) Photo credit: [50]. . . . . . . . . . . . . . . . . . . . . . . 35 3.2 Blue lines: results reproduced from averaging 8000 realizations. The initial state of the uctuator is either state 1 or1 with equal probability. Green lines: analytical solutions. Note that b = 5 for both curves. . . . . . . . . . . 38 x 3.3 Saturation of decay asb further increases. Note that = 1 switch/(time unit) for all curves. The decay rate is saturated (cannot decay faster) afterg = 2b= reaches 1. The curves are produced from averaging 10000 realizations with a random initial uctuator state for each realization. . . . . . . . . . . . . . . 39 3.4 (a) A cartoon illustrating a single qubit coupled to a series of uctuators with dierent individual ipping frequencies (b) Distribution of frequencies of 12000 uctuators. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 41 3.5 Reproduced gures of 2 to 5 decades. Couplings distributed withhbi=hbi = 0:2 and aroundhbi = 4:6 10 7 Hz. max = 10 THz with dierent min . n d = 1000. dp 0j =1 are distributed according tohdp 0j i =dp eq = 0:08. . . . 42 3.6 Top: min = 10 6 GHz, max = 100 GHz, n d = 100,hbi = 10 4 =4; Bottom: min = 10 4 GHz, max = 100 GHz, n d = 100,hbi = 0:01=4. Results are averaged over 8k realizations. Each uctuator is initially in a pure state randomly sampled according tohp 0 i =p eq = 0:08. b=hbi 0:2. . . . . . . 43 3.7 12 decades of noise: max = 10 12 Hz with min = 10 0 = 1 Hz. Couplings distributed withhbi=hbi = 0:2 and aroundhbi = (4:6=)10 7 Hz. n d = 1000. dp 0j =1 are distributed according tohdp 0j i = dp eq . As time doubles, the reduced decay factor drops more than 2 times. . . . . . . . . . . . . . . . . . 44 3.8 g = 2b= . Averaged over 10000 trajectories. . . . . . . . . . . . . . . . . . . 47 3.9 Pure dephasing noise: min = 0:1 GHz, max = 10 GHz. Number of uctuators is 10.hbi = 0:2 GHz, b=hbi = 0:2. Results are averaged over 30k trajectories. The initial value of each uctuator is randomly sampled according to the thermal equilibrium distribution p eq = 0:08 (i.e. each uctuator is 0:08 more likely to be initialized in state 1 than state1). (a) Instantaneous Ground state population, (b) (Absolute value) of o-diagonal elements in the computational basis. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 48 3.10 High temperature limit, pure dephasing noise: min = 0:01 GHz, max = 1 GHz. Number of uctuators is 20.hbi = 0:2 GHz, b=hbi = 0:2. Results are averaged over 30k trajectories. The initial value of each uctuator is randomly sampled according to thermal equilibrium distribution p eq = 0 (i.e. each uctuator is equally likely to be initialized in state 1 than state1), since T!1. (a) Instantaneous Ground state population, (b) (Absolute value) of o-diagonal elements in the computational basis. . . . . . . . . . . . . . . . . 50 3.11 Zero temperature limit, pure dephasing noise. Same uctuator parameters as in Fig. 3.10. Results are averaged over 30k trajectories. p eq = 1, since T! 0. (a) Instantaneous Ground state population, (b) (Absolute value) of o-diagonal elements in the computational basis. . . . . . . . . . . . . . . . 51 xi 3.12 Transverse noise: min = 0:01 GHz, max = 1 GHz. Number of uctuators is 200.hbi = 0:046 GHz, b=hbi = 0:2. Results are averaged over 1k trajectories. The initial value of each uctuator is randomly sampled according to thermal equilibrium distributionp eq = 0:08 (i.e. each uctuator is 0:08 more likely to be initialized in state 1 than state1). Error bars are 2 over 1k trajectories. 52 3.13 Longitudinal noise: min = 0:01 GHz, max = 1 GHz. Number of uctuators is 100. hbi = 0:046 GHz, b=hbi = 0:2. Results are averaged over 1k trajecto- ries. The initial value of each uctuator is randomly sampled according to the thermal equilibrium distributionp eq = 0:08 (i.e. each uctuator is 0:08 more likely to be initialized in state 1 than state1). Here we denote the decay prole ase t=T 2 since there is no longitudinal noise (i.e. T 1 =1). Error bars are 2 over 1k trajectories. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 54 3.14 Longitudinal noise: min = 0:01 GHz, max = 1 GHz. Number of uctuators is 20. Results are averaged over 3k trajectories. The initial value of each uctuator is randomly sampled according to the thermal equilibrium distri- butionp eq = 0:08 (i.e., each uctuator is 0:08 more likely to be initialized in state 1 than state1). Blue line: hbi = 0:2 GHz; Green line: hbi = 2 GHz. b=hbi = 0:2. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 55 3.15 Schematic of the experimental set up. Capacitively shunted ux qubit (CSFQ) in the center, with curved arrows showing x and z uxes threading their corresponding loops. Source: Ref. [59]. . . . . . . . . . . . . . . . . . . . . . 57 3.16 Spectrum with z (0)= =1:21. . . . . . . . . . . . . . . . . . . . . . . . . 58 3.17 The eigenstate populations for z (0)= =1:21. (a) Closed-system evolution, (b) Evolution under the adiabatic master equation. (c) 1=f noise: min = 0:01 GHz, max = 1 GHz. Number of uctuators is 100. hbi = 0:46 GHz, b=hbi = 0:2. Results are averaged over 1k trajectories. The initial value of each uctuator is randomly sampled according to the thermal equilibrium distribution p eq = 0:08. Some of the population is lost to higher excited states. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 60 3.18 Spectrum with x = 0:5. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 62 3.19 The population inj+i for coupling (a) g = 0:008(2). (b) g = 0:08(2). Evolution under (adiabatic) master equation. . . . . . . . . . . . . . . . . . . 63 3.20 A = z . The population ofj+i for x = =0:5. min = 1 GHz, max = 10 GHz. Number of uctuators is 20.hbi = 0:4600 GHz, b=hbi = 0:2. Results are averaged over 300 trajectories. p eq = 0:08. . . . . . . . . . . . . . . . . . 64 3.21 A =I p . The population ofj+i for x = =0:5. min = 0:1 GHz, max = 10 GHz. Number of uctuators is 20.hbi = 0:0092 GHz, b=hbi = 0:2. Results are averaged over 300 trajectories. p eq = 0:08. . . . . . . . . . . . . . . . . 65 xii 3.22 The processes of determining the spectrum [56] and performing simulations based on the spectrum. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 66 3.23 t f = 5405ns,b = 0:12(9:2)10 3 ; min = 0:05 GHz; max = 5 GHz, uctuations in I p , 20 uctuators, 200 trajectories. Free induction decay. Rotated frame. . 68 3.24 Fluctuations in [I p ] 22 , 20 uctuators, 300 trajectories. Rotated frame. min = 0:01 GHz; max = 0:1 GHz. (a), (b): b = 9:2929 10 4 GHz, (c) Modied averaged uctuator strength: b = 1:05 10 3 GHz. . . . . . . . . . . . . . . 69 4.1 By successively interacting the system with environmental (probe) qubits which are subsequently measured, we can obtain the weak measurement tra- jectories. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 73 4.2 Photodetection current I(t). . . . . . . . . . . . . . . . . . . . . . . . . . . . 74 4.3 Evolution of instantaneous ground state population for a total time of t f = 100s and temperature 2:6 GHz, with feedback applied with some dierent delays (0; 5; 10; 30s). Results are obtained from averaging 5k quantum trajectories of the quantum feedback algorithm. Error bars are 2 over 5k trajectories. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 77 4.4 Linear Schedule. Final ground state population. Single qubit. Total annealing time t f = 100s. Dashed purple line is the nal ground state population without feedback. Average over 5k trajectories with Lindbladian feedback. There is a non-small (thus experimentally allowed) 60us where the Final GS population is maximized. . . . . . . . . . . . . . . . . . . . . . . . . . . . 78 4.5 Evolution of ground state population under feedback error correction of dier- ent feedback Hamiltonians ( = 0). The results of quantum jump trajectories (averaged over 5k trajectories) of quantum feedback algorithm agree with the feedback master equation solutions. Error bars are 2 over 5k trajectories. . 80 5.1 (a) Annealing schedules (in units such that ~ = 1) as a function of the an- nealing fraction s(t), chosen to be similar to the schedules of the D-Wave processors. (b) Annealing fraction s(t). The blue solid curve represents stan- dard forward quantum annealing of total annealing time = 500ns. . . . . . 83 5.2 (a) Success probability in unitary reverse annealing as a function of the in- version point s inv , for several values of the magnetization of the initial state. The dashed vertical line indicates s inv = s 0:309. The annealing time is = 100 ns. We sampled the interval s inv 2 (0; 1) using a step size of s = 0:002, and repeated the dynamics for each choice ofs inv . (b) Maximum success probability achievable with unitary reverse annealing, as a function of the magnetization of the initial state, for two annealing times: = 100 ns and = 1000 ns. The inset zooms in on the region m2 [1:0; 0:0]. Note the logarithmic scale on the vertical axis. . . . . . . . . . . . . . . . . . . . . . . 87 xiii 5.3 Success probability in reverse annealing as a function of the inversion point s inv , for several values of the magnetization of the initial state. Thep-spin sys- tem is coupled to a collective dephasing bosonic environment as in Eq. (5.10), and the coupling strength is = 1 10 3 . The dashed vertical line denotes the time s of the avoided crossing between the ground state and the rst excited state. In (a) the annealing time is = 100 ns, in (b) = 500 ns. We sampled the intervals inv 2 (0; 1) using a step s = 0:02. All other parameters are given in the main text. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 90 5.4 Success probability in reverse annealing as a function of the inversion point s inv , for N 2f3;:::; 8g. The initial state is the rst excited state of the maximum spin subspace (m 0 = 1 2=N). The dashed vertical line denotes the time s of the avoided crossing between the ground state and the rst excited state. The annealing time is = 100 ns. We sampled the interval s inv 2 (0; 1) using a step size of s = 0:005. . . . . . . . . . . . . . . . . . . 101 5.5 Maximum success probability achievable with reverse annealing as a function of the number of qubitsN, using the collective and the independent dephasing models of Eqs. (5.10) and (5.12), respectively. The annealing time is = 100 ns.102 5.6 Success probability in paused reverse annealing for the collective dephasing model, as a function of the inversion point s inv , for several values of the mag- netization of the initial state. The coupling strength is = 1 10 3 . The annealing time is = 100 ns. In (a) a pause of durationt p = 100 ns is inserted at the inversion point, while in (b) t p = 400 ns. All other parameters are given in the main text. Dotted lines with empty symbols refer to open sys- tem reverse annealing of time = 100 ns and no pauses. The dashed vertical line denotes s inv = s . The interval s inv 2 (0; 1) is sampled using a step size s = 0:02. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 102 5.7 Comparison of the success probability in reverse annealing for the collective and independent dephasing models, as a function of the inversion point s inv , forN2f3;:::; 8g. The initial state is the rst excited state of the maximum spin subspace (m 0 = 1 2=N). The dashed vertical line denotes the time s of the avoided crossing between the ground state and the rst excited state. The annealing time is = 100 ns, and a pause of duration t p = 100 ns is inserted at the inversion point. . . . . . . . . . . . . . . . . . . . . . . . . . . 103 6.1 (a) Annealing schedule of D-Wave 2000Q with the \DW 2000Q 6\ solver. (b) Forward (FA) and reverse annealing (RA) protocols are dened in Eq. (6.4)) and Eq. (6.5)), respectively. The latter incorporates a pause of duration t p . The parameter values for this plot aret a = 1:2 [s] for forward annealing, and = 1 [s] and s inv = 0:2 for reverse annealing. These symbols are dened below Eq. (6.5). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 107 xiv 6.2 Success probability for dierent initial conditions m 0 (magnetization) in re- verse annealing on the D-Wave 2000Q device as a function of the inversion points inv . The dashed line is the minimum-gap point ats 0:36 forN = 20. Whens inv <s the reverse direction of the anneal goes through and past the minimum gap point, and then crosses it again during the forward anneal. When s inv > s there is no crossing of the minimum gap point. (a) Partial success probability for m 0 = 0; 0:8, and 0:9. (b) Total success probability for the same set of m 0 as in (a). Here and in all subsequent gures the labels \up" and \down" mean the populations of the all-up state and the all-down state, respectively. Error bars denote standard deviations. . . . . . . . . . . 111 6.3 Success probability for dierent annealing times in reverse annealing on the D-Wave 2000Q device as a function ofs inv . (a) Partial success probability for = 5; 10, and 20[s]. (b) Total success probability for the same set of as in (a). The dashed line is the minimum-gap point at N = 20, s 0:36. . . . . 114 6.4 Success probability for dierent pausing times t p in reverse annealing as a function ofs inv . (a) Partial success probability for pauset p = 0; 2, and 10 [s]. (b) Total success probability for the same set of t p as in (a). The dashed line is the minimum-gap point at N = 20, s 0:36. . . . . . . . . . . . . . . . . 115 6.5 Success probability for dierent system sizesN in reverse annealing as a func- tion ofs inv . (a) Partial success probability for system sizeN = 4; 8; 12; 16; 20, and 24. (b) Total success probability for the same set of N as in (a). (c) The position of the minimum gap for 4N 22, along with the value of s inv at which the success probability =:5 for eachN. Inset: the minimum gap for each N. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 116 6.6 Success probability for dierent number of iterations r in iterated reverse annealing as a function of s inv . (a) and (b) are partial and total success probabilities for number of iterationsr = 1; 10; 25, and 50 with the initial state m 0 = 0:9 (the rst excited state), respectively. (c) and (d) are partial and total success probabilities with the initial state m 0 = 0 (the highest excited state). The dashed line is the minimum-gap point at N = 20, s 0:36. . . 118 6.7 Total success probability for closed system with r = 1 (single cycle) and a single-spin down as the initial state. (a) N = 4, (b) N = 8. Here and in the subsequent gures the nanoscecond timescale is set by the energy scale of the D-Wave annealing schedule shown in Fig. 6.1. . . . . . . . . . . . . . . . . . 121 6.8 Total success probability with two spins down as the initial state. (a) N = 4 and r = 1, (b) N = 8 and r = 1, (c) N = 4 and r = 2, (d) N = 8 and r = 2. 122 6.9 Top row: Success probability as a function of s inv with (a) independent and (b) collective dephasing (SB denotes system-bath coupling). Initial state: j0001i (m 0 = 0:5) andr = 1, a single cycle. Bottom row: Success probability with dierent . Initial state:j0011i (m 0 = 0). (a) r = 1, (b) r = 2. . . . . . 127 xv 6.10 Success probability as a function ofs inv for independent dephasing. The initial state has one spin down and is 1[s]. The main plot shows the results for N = 4 and 8. The inset shows the partial success probabilities for N = 8. . . 129 6.11 Total success probability for SVMC and SVMC-TF. (a) = 10 3 sweeps, (b) = 10 4 sweeps. Error bars are 2 over 10 4 samples. . . . . . . . . . . . . . . 132 6.12 Partial success probability by (a) SVMC, N = 4and (b) SVMC-TF, N = 4. = 10 3 sweeps. (c) SVMC-TF, N = 8. = 10 3 sweeps. Error bars are 2 over 10 4 samples. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 133 6.13 Partial success probability for SVMC-TF. (a) N = 16, = 10 3 sweeps, (b) N = 16, = 10 4 sweeps, (c) N = 32, = 10 3 sweeps, (b) N = 32, = 10 4 sweeps. Error bars are 2 over 10 4 samples. . . . . . . . . . . . . . . . . . . 135 7.1 The 16 lowest lying energies of 10 qubits, with H init having N # = 2. . . . . . 139 7.2 (s) for =f1; ; 5g in units where h = 1. . . . . . . . . . . . . . . . . . 140 7.3 (s) for dierent i . = 1 in units where h = 1. . . . . . . . . . . . . . . . . 140 7.4 for N =f3; ; 20g in units where h = 1. Inset: the value of s min for each N. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 141 7.5 Open-system and close-system adiabatic reverse annealing (ARA) simulation results. = 40ns. = 1. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 145 7.6 Open-system and close-system adiabatic reverse annealing (ARA) simulation results. = 4ns. = 1. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 146 7.7 Open-system and close-system adiabatic reverse annealing (ARA) simulation results. = 400ns. = 0:3. . . . . . . . . . . . . . . . . . . . . . . . . . . . 147 7.8 ARA performance for (a) closed system, (b) open system with independent system-bath coupling, (c) open system with collective system-bath coupling, (d) comparison of the 3 situations in terms of TTS (N = 10). . . . . . . . . . 148 7.9 Standard QA in an open system with two dierent system-bath coupling model.149 7.10 The value of minimum gap as a function of H init and . N = 4. . . . . . . . 150 7.11 ARA open results (a) = 1, (b) = 2, (c) = 5. Independent system-bath coupling. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 151 7.12 ARA open results, with initial state beingj##""i. (a) success probability, (b) TTS. Independent system-bath coupling. . . . . . . . . . . . . . . . . . . . . 152 7.13 ARA Closed (a) = 0:1, (b) = 0:5. N = 4. . . . . . . . . . . . . . . . . . . . 153 7.14 ARA Open . (a) = 0:1, (b) = 0:5. N = 4. Independent system-bath coupling. 153 xvi 7.15 (a) N = 8, (b) N = 10. = 0:1 and N # = 2. Independent system-bath coupling. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 154 7.16 Dependence of success probability as a function of N due to the two dierent SB coupling model. c = 0:5, = 1 and = 10ns. . . . . . . . . . . . . . . . 156 8.1 (a): Dotted lines: Linear schedule; Solid lines: D-Wave schedule used in Fig. 2.3; (b): Same as in Fig. 2.3 (the evolution of the population in the instantaneous ground state for the 8-qubit problem), but with a linear sched- ule. Inset: the convergence of the ground state population towards the AME results. (c): A histogram of the net number of jumps to the instantaneous ground state (GS). It shares the same pattern as in the non-linear schedule case (the inset of Fig. 2.4), but the net number of jumps out of the ground state is smaller due to the smaller energy scale in the linear schedule. . . . . 179 8.2 Full spectrum of 4 qubits. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 184 xvii Abstract In this dissertation, we explore how quantum annealing (QA) and its applications behave in an open system setting. We give derivations and numerical recipes for eective parallel simulation methods for time-dependent open dynamics of quantum annealing devices. We consider the weak-coupling limit to an environment and also the case of 1=f noise. The stochastic time-dependent quantum trajectories technique can be utilized in studying weak measurements and feedback error correction in quantum annealing. Then we focus on open- system descriptions of reverse annealing (RA), which is a promising variant and application of quantum annealing. We show that, with various simulation tools, reverse annealing can benet from the interaction between the annealer and its environment. xviii Chapter 1 Introduction Quantum annealing (QA) [1, 2] and the broader class of adiabatic quantum computing (AQC) [3] are quantum approaches to solving combinatorial optimization problems and sampling problems. Quantum annealing obtains the global solution(s) by using quantum uctuations to escape from local minima [1, 3{7]. Quantum annealing is generally considered to be more noise resilient than gate-based quantum computing, and there are many ongoing developments (e.g. [8]) on how to further improve the coherence of quantum annealers. Accurate open quantum system theories of QA along with ecient simulation procedures serve at least two purposes: 1) to aid in the design of experimental quantum annealers and derivation of annealing protocols/applications such as advanced annealing schedules and reverse annealing, 2) to help the development of error correction techniques tailored to noisy quantum annealers. The dynamics in quantum annealers is driven by a time-dependent system Hamiltonian, which in this thesis We will refer to specically as the annealing Hamiltonian. In standard forward quantum annealing, it takes the following form: H(s) = A(s) 2 V TF + B(s) 2 H p : (1.1) 1 0 0.2 0.4 0.6 0.8 1 0 2 4 6 8 10 12 14 A(s) B(s) Figure 1.1: D-Wave 2000Q NASA annealing schedules. V TF = P N i i x is the transverse eld, N is the number of qubits, and H p is the problem Hamiltonian. Here s = t=t f is the normalized time, with t f the total physical annealing time. The functions A(s) and B(s) are the annealing schedules (see for example Fig. 1.1.) One metric for measuring the success of quantum annealing is the ground state population at the end of the anneal. This is equivalent to the success probability of reaching the correct solution, often simply called the success probability. One can model and estimate this population by a variety of open-system analysis tools. There are many previous studies about the open-system modeling of quantum annealing [9{12]. Notably, a time-dependent adiabatic master equation was developed to describe decoherence of quantum annealers in the weak-coupling limit [12, 13]. Under this master equation the decoherence occurs in the instantaneous eigenbasis of H(s). However, it suers from two problems: 1. Solving the master equation for a density matrix requires storing and updating N 2 1 real numbers, which becomes infeasible for largeN; 2. The assumed weak-coupling limit does not apply to slow noise such as 1=f noise, which is present, e.g., in superconducting-qubit based quantum annealers. While there exists literature [10, 14] addressing the modeling of slow noise in current annealing device, knowledge about how in particular the 1=f power spectrum aects QA and ecient modeling of 1=f-noise for QA is lacking. The rst two chapters of this 2 dissertation propose new and original solutions to these two problems. The rst chapter proposes a time-dependent quantum trajectories method adapted to QA, which gives a factor ofN reduction in the simulation cost. The cost of running many trajectories can further be minimized with parallel computing. The quantum trajectories method for QA also provides insight into individual quantum jump trajectories and their jump statistics, thus shedding light on open system quantum adiabatic evolution beyond the master equation. The second chapter proposes stochastic trajectories methods for the slow noise of 1=f power spectrum, with the time-dependent system Hamiltonian of QA. We also analyze how dierent elements such as temperature and 1=f noise range play a role in the decoherence in QA. Moreover, past research in weak measurement and feedback-based error correction fo- cused on time-independent Hamiltonians and gated-based quantum computing. Based on the time-dependent trajectories techniques we develop, we further investigate how to perform weak measurement and feedback error correction on quantum annealers. Standard forward annealing with the annealing Hamiltonian (Eq. 6.27) has the major drawback of suering from thermal excitations, which take a certain portion of population out of the solution states at the end of the anneal. There care at least two remedies to this drawback. One approach is to start in a random excited state of the problem Hamiltonian H p and simply make use of relaxation to get to the ground state of H p , which has been studied in [15]. However, one can do more. Unlike forward annealing, where s runs from 0 to 1, one can start from an excited state of the problem Hamiltonian H p at s = 1, 1 reverse anneal to an inversion point s = s , and then forward anneal back to s = 1. This approach combines both relaxation from excited states and annealing, and allows most of the relaxation to occur at the smaller gap at s =s . As for the second remedy, we observe that the annealing Hamiltonian (Eq. 6.27) is not the unique form of a Hamiltonian drive for quantum annealers. One can proceed with standard forward annealing, but with an extra term in Eq. (6.27), sometimes called initialization Hamiltonian or catalyst Hamiltonian [16{ 1 In Eq. (6.27), we do not assume that s goes from 0 to 1. 3 18]. This additional term provides an easier annealing path to the solution and thus leads to a higher population in the solution state at the end of the anneal. The success of adding the initialization/catalyst Hamiltonian often depends on how close it is to the correct solution. The two remedies can be summarized within the general framework of reverse annealing. The rst is \reversed" in the sense that it starts at s = 1; the second is \reversed" in the sense that the catalyst term depends on our knowledge about the solution states at s = 1. Reverse annealing [16, 19{25] in general can therefore be separated into two formalisms: iterated reverse annealing (IRA) and adiabatic reverse annealing (ARA). As mentioned be- fore, one can start at s = 1 with an initial state that is usually the excited state of H p , anneal reversely to an inversion point s = s , and then anneal forwardly to s = 1. The procedure can be performed more than once and may be iterated with the last output state as the new input state. It is called iterative reverse annealing (IRA) and its function has been included in the current generation of quantum annealers built by D-Wave [26]. Adia- batic reverse annealing (ARA) is similar to standard forward annealing, but with an extra initialization/catalyst Hamiltonian. The two formalisms have been studied previously in a closed system setting. In this dissertation we address the problem of how to model and predict their performance in an open-system settings, which has so far been lacking in the literature. We show that both reverse annealing variants can benet from the interaction between the annealer and the bath, in that the nal ground state population can be higher than in a closed system setting. Therefore, a quantum annealer can take advantage of noise to improve its performance in the sense of increased success probabilities. This dissertation is organized as follows. In chapter 2, we derive the time-dependent quantum trajectories solution to adiabatic master equation, which yields a quadratic sav- ing in time cost compared to the master equation approach. In chapter 3, we extend the stochastic trajectory formulation to the simulation of slow noise with a 1=f power spectrum. In chapter 4, we describe and study the weak measurement and feedback error correction protocol tailored to QA. Chapters 5-7 focus on reverse annealing applications. In Chapter 5, 4 we focus on the open-system simulation of iterative reverse annealing (IRA) and show how it benets from environmental-induced relaxation. In chapter 6, we study the iterative reverse annealing protocol in D-Wave and perform open-system simulations that agree with the ex- perimental data. In chapter 7, we focus on the open-system description of another variant of reverse annealing called adiabatic reverse annealing (ARA). In chapter 8 we provide a brief summary along with an outlook and perspective on the signicance of the results obtained in this thesis. 5 Chapter 2 Quantum trajectories for time-dependent adiabatic master equations 2.1 Introduction In this chapter we describe a quantum trajectories technique for the unraveling of the quan- tum adiabatic master equation in Lindblad form. The goal is that by evolving a complex state vector of dimension N instead of a complex density matrix of dimension N 2 , simulations of larger system sizes become feasible. The cost of running many trajectories, which is required to recover the master equation evolution, can be minimized by running the trajectories in parallel, making this method suitable for high performance computing clusters. In general, the trajectories method can provide up to a factor N advantage over directly solving the master equation. In special cases where only the expectation values of certain observables are desired, an advantage of up to a factor N 2 is possible. We test the method by demon- strating agreement with direct solution of the quantum adiabatic master equation for 8-qubit quantum annealing examples. We also apply the quantum trajectories method to a 16-qubit example originally introduced to demonstrate the role of tunneling in quantum annealing, which is signicantly more time consuming to solve directly using the master equation. The quantum trajectories method provides insight into individual quantum jump trajectories and 6 their statistics, thus shedding light on open system quantum adiabatic evolution beyond the master equation. This chapter is based on [27]. With the growing ability to control and measure ever-larger quantum systems, under- standing how to model the interactions between open quantum systems and their environ- ment has become exceedingly important [28]. The open system dynamics is often described in terms of a master equation in Lindblad form, describing the eective dynamics of the quantum system after the environmental degrees of freedom have been traced out [29]. An equivalent approach is that of quantum trajectories (also known as the Monte Carlo wave- function method) [30{32], which can be understood as an unraveling of the master equation in Lindblad form, and which generates a stochastic process whose average is fully equivalent to the master equation (for a review, see Ref. [33]). Each trajectory in this approach can also be viewed as the result of continuous indirect measurements of the environment in a certain basis [34]. A quantum trajectories approach exists also for non-Markovian master equations [35, 36]. While a vast literature exists on the topic of quantum trajectories for time-independent master equations, much less is known for the case of time-dependent master equations (see, e.g., Ref. [37]), which is our focus here. Specically, we focus on the case of open systems evolving adiabatically according to a time-dependent Hamiltonian, weakly coupled to the environment [13, 38]. This is particularly relevant in the context of quantum annealing and more generally adiabatic quantum computing, whereby the computation proceeds via a time-dependent Hamiltonian and the result of the computation is encoded in the ground state of the nal Hamiltonian (for reviews see Refs. [5, 17]). A large body of literature exists on the use of quantum Monte Carlo methods in the context of quantum annealing (see, e.g., Refs. [4, 39, 40]), but these methods focus on equilibrium properties, while here we are interested in dynamics. The interplay between the key quantities that determine adiabaticity and non-unitary dynamics has not been previously explored in the setting of Monte Carlo wavefunction methods, and here we resolve this question by nding an upper bound on the 7 size of the Monte Carlo time-step. Nor has the question of reducing the computational cost of simulations of the adiabatic master equation via quantum trajectories been discussed so far, and we address this here. Thus, here we develop the rst treatment of a quantum trajectories unravelling of a time-dependent adiabatic master equation (AME). We make a formal comparison between the quantum trajectory unraveling of the Lindblad master equation with time-independent and time-dependent operators, and discuss the validity of applying it to the unraveling of the AME. While our analysis closely follows the standard time-independent approach, the time-dependent case results in new additional validity conditions that must be satised. The individual trajectories of the AME shed new light on how the average case captured by the AME emerges. When the quantum state is a pure energy eigenstate and the unitary evolution is adiabatic, the drift term in the quantum trajectories approach vanishes, and the quantum state follows the instantaneous eigenstates until a jump occurs. We can associate these jumps with an excitation or relaxation process, depending on the direction of the jump. This provides an intuitive (yet rigorous) picture for how the averaged dynamics of the AME arises. An important advantage of the quantum trajectories approach is that for anN-dimensional system, one quantum trajectory requires storing and updating 2N 1 real numbers, while solving the master equation for a density matrix requires storing and updating N 2 1 real numbers. This quadratic saving allows simulations of systems with sizes that are infeasible by directly solving the master equation. The tradeo is that many trajectories must be run in order to accurately approximate the solution of the master equation, but this tradeo can be reduced by using many parallel processes to represent each trajectory. Our presentation is organized as follows. In Sec. 2.2 we brie y review the AME. We unravel the AME in Sec. 2.3 into quantum trajectories taking the form of quantum jumps, allowing for an arbitrary time-dependence of the Hamiltonian and Lindblad operators. In Sec. 2.4 we provide an algorithmic implementation for our adiabatic quantum trajectories 8 and in Sec. 2.5 we present three case studies. We perform a cost comparison between the direct simulation of the AME and the quantum trajectories method in Sec. 2.6. Additional technical details and proofs are provided in the Appendices. 2.2 Adiabatic master equation in Lindblad form We focus on the AME in Lindblad form, which can be derived with suitable approximations (in the weak coupling limit after performing the Born-Markov, rotating wave, and adiabatic approximation) from rst principles starting from the system Hamiltonian H S , the environ- ment Hamiltonian H B , and the interaction Hamiltonian H I = g P A B , with system operators A , environment operators B , and system-bath coupling strength g [13]. The adiabatic (Lindblad) master equation describes the evolution of the system density matrix (t) and has d dt (t) =i [H S (t) +H LS (t);(t)] +L WCL [(t)] ; (2.1) where H LS (t), which commutes with H S (t), is a Lamb shift Hamiltonian arising from the interaction with the environment. The dissipative termL WCL takes the form: L WCL [(t)] X ; X ! (!) L ;! (t)(t)L y ;! (t) 1 2 L y ;! (t)L ;! (t);(t) ; (2.2) where the sum over ! is over the Bohr frequencies (eigenenergy dierences) of H S , (!) is an element of the positive matrix , and satises the Kubo-Martin-Schwinger (KMS) condition if the bath is in a thermal state with inverse temperature = 1=T : (!) =e ! (!) : (2.3) 9 The time-dependent Lindblad operators are given by: L ;! (t) = X a;b !;" b (t)"a(t) h" a (t)jA j" b (t)ij" a (t)ih" b (t)j ; (2.4) wherej" a (t)i is the a-th instantaneous energy eigenstate of H S (t) with eigenvalue " a (t). With this form for the Lindblad operators, decoherence can be understood as occurring in the instantaneous energy eigenbasis [12]. For the purpose of unravelling the above master equation into quantum trajectories, it is convenient to diagonalize the matrix by an appropriate unitary transformation u(!): X ; u i; (!) (!)u j; (!) y = 0 B B B B @ 0 1 (!) 0 ::: 0 0 2 (!) ::: . . . . . . . . . 1 C C C C A i;j ; (2.5) and to dene new operators A i;! (t) given by L ;! (t) = X i u i; (!)A i;! (t) : (2.6) In this basis, we can write the dissipative part in diagonal form as: L WCL [(t)] = X i X ! 0 i (!) A i;! (t)(t)A y i;! (t) 1 2 n A y i;! (t)A i;! (t);(t) o : (2.7) 2.3 Stochastic Schr odinger Equation With Eq. (2.7), the master equation Eq. (2.1) is in diagonal form and can be unravelled into quantum trajectories. The trajectory is described by a stochastic dierential equation (SDE) in the form of jumps or diusion. Let us consider the case where the coecients 10 0 i (!) in Eq. (2.7) also depend on time. If all 0 i (!;t) 0, then the dynamics is completely positive (CP)-divisible [41], and the master equation can be unravelled by using the known unravelling of the the time-independent SDE case [29, 34, 42], simply by replacing the time- independent operators and coecients by the time-dependent ones. Such an unravelling is also possible, but with modications, when the dynamics is positive (P)-divisible, i.e., where 0 i (!;t) need not be all positive. 1 This can be in the form of: • Jump trajectories: the master equation is unravelled via the non-Markovian quantum jump method (NMQJ) [43{46], where terms with negative coecients 0 i (!;t) describe the negative channel. • Diusive trajectories: recent work on diusive trajectories [37] replaces 0 i (!;t) and the operators by the eigenvalues and eigenvectors of a positive transition rate operator W (Eq. (11) in [37]). P-divisible dynamics can be unravelled into a SDE in terms of such eigenvalues and eigenvectors. In the following, we focus on the case of CP maps [with all 0 i (!) 0] and unravel the master equation in the quantum jumps picture. 2.3.1 Unravelling the master equation First we absorb the 0 coecients into the denition of A i : p 0 i (!)A i;! (t)!A i (t) : (2.8) 1 The condition on 0 i (!;t) such that the map is P-divisible can be found in the proof given in [43] or Eq. (25) in [37]. 11 In this redenition, the index i now includes the Bohr frequencies. We write Eq. (2.1) in terms of an eective non-Hermitian Hamiltonian H e : d dt (t) =i H e (t) S (t) S (t)H y e (t) + X i A i (t) S (t)A y i (t); (2.9) where H e (t) =H S (t) +H LS (t) i 2 X i A y i (t)A i (t): (2.10) Equation (2.9) can be unravelled into quantum trajectories in the quantum jumps picture, where each trajectory describes the stochastic evolution of a pure state (if the initial state is mixed, = P i p i j i ih i j, then the evolution can be performed on each initial pure state). The stochastic evolution of the pure state can be written in terms of a stochastic Schr odinger equation (in It^ o form), the ensemble average of which is equivalent to the master equation: dj (t)i = iH e (t) + 1 2 X i hA y i (t)A i (t)i ! dtj (t)i + X i dN i (t) 0 @ A i (t) q hA y i (t)A i (t)i 1 1 A j (t)i (2.11) where hA y i (t)A i (t)ih (t)jA y i (t)A i (t)j (t)i =kA i (t)j (t)ik 2 : (2.12) We give a derivation of Eq. (2.11) below. The rst term on the r.h.s. of Eq. (2.11) gives a deterministic evolution composed of a Hermitian component [i(H S (t) +H LS (t))] and a \drift" component D(t) 1 2 X i A y i (t)A i (t)hA y i (t)A i (t)i ; (2.13) and the second term describes the stochastic jump process. The stochastic variabledN i (t) N i (t +dt)N i (t) is the number of jumps of typei in the intervaldt, where we have denoted 12 byN i (t) the number of jumps of typei up to timet. The expectation value of the stochastic variable is given by [29]: E[dN i (t)] =hA y i (t)A i (t)idt: (2.14) Since the probability of a jump occurring scales linearly with dt, the probability of having more than one jump vanishes faster than dt, so as dt! 0 only one jump out of all possible types during dt is permitted. Therefore we can write [42]: dN i (t) = 8 > < > : 1 with prob. hA y i (t)A i (t)idt 0 with prob. 1hA y i (t)A i (t)idt (2.15) with the It^ o table: dN j (t)dN k (t) = jk dN j (t) (2.16a) dN j (t)dt = 0: (2.16b) From Eq. (2.15), the probability of any jump occurring, P i hA y i (t)A i (t)idt, is small compared to the probability of no jump occurring, so P i dN i (t) = 0 most of the time. During the innitesimal time-step dt, if P i dN i (t) = 0, then only the deterministic evolution takes places; if however P i dN i (t) = 1, then a jump occurs. When a jump occurs, it dominates over the deterministic evolution, which is proportional to dt, and the deterministic part can be ignored. 2.3.2 Deterministic evolution and jump process We now derive Eq. (2.11) by explaining how each probability element appears. Let us denote byj (t)i andj ~ (t)i the normalized and unnormalized state vectors respectively, and assume they are equal at time t, i.e.,j ~ (t)i =j (t)i. For the innitesimal time-step from t to t +dt, the state vector evolution fromj (t)i toj (t +dt)i involves two possibilities: either 13 |ψ(t)i |ψ(t+dt)i |ψ(t+dt)i no jump (1−dp(t)) jump (dp(t)) Figure 2.1: A depiction of the stochastic evolution of the statej i by an innitesimal time state at time t. no jump occurring (with probability 1dp) or a jump occurring (with probabilitydp). This is depicted in Fig. 2.1. When no jump occurs, the evolution is described by the Schr odinger equation associated withH e , and since the eective Hamiltonian is non-Hermitian the norm of the state vector is not preserved during the evolution: dj ~ (t)i dt =iH e (t)j ~ (t)i : (2.17) The resulting state after one innitesimal time-step dt is: j ~ (t +dt)i = exp [iH e (t)dt]j (t)i (2.18a) = IidtH e (t) +O(dt 2 ) j (t)i : (2.18b) The norm squaredk ~ (t +dt)k 2 is the probability of the conditional evolution under H e , so that (as we show explicitly in Appendix A.1) the jump probability is given by 1k ~ (t +dt)k 2 =dt X i hA y i (t)A i (t)i +O(dt 2 ) (2.19) 14 [recall Eq. (2.12) and note that H S (t) +H LS (t) cancels out to rst order]. Therefore we can identify the innitesimal jump probability dp(t) with the r.h.s. of Eq. (2.19), i.e., to rst order in dt: dp(t) = X i dp i (t) =dt(t) (2.20a) dp i (t) =dthA y i (t)A i (t)i ; (2.20b) where (t) = _ p(t) is the jump rate, and dp i (t) is the probability of the jump of type i. Note that since our denition of theA i operators includes the rates 0 [recall Eq. (2.8)], the jump rate depends on the instantaneous Bohr frequencies and the KMS condition. When the jump of type i occurs the state is updated as: j ~ (t +dt)i =A i (t)j (t)i : (2.21) We can unify the two possibilities in Eq. (2.18) and Eq. (2.21) as a stochastic Schr odinger equation for the unnormalized state vector where only terms of order dt are kept: dj ~ (t)i =j ~ (t +dt)ij ~ (t)i (2.22a) =idtH e (t)j ~ (t)i + X i dN i (t) (A i (t)1)j ~ (t)i ; (2.22b) where we usedj ~ (t)i = j (t)i. The stochastic element dN i (t) has the properties given in Eqs. (2.15) and (2.16); 1 is subtracted since when the jump occurs P i dN i (t)j ~ (t)i = j ~ (t)i and the term involving H e (t) is absent, so in this manner we ensure thatj ~ (t)i is appropriately subtracted from the r.h.s. 15 We can write a similar expression for the normalized state vectorj (t)i by normalizing Eqs. (2.18) and (2.21). If a deterministic evolution occurs, we have j (t +dt)i = exp [iH e (t)dt]j (t)i kexp [iH e (t)dt]j (t)ik (2.23a) = (1idtH e (t) +O(dt 2 ))j (t)i q 1dt P i hA y i (t)A i (t)i +O(dt 2 ) (2.23b) = 1idtH e (t) + 1 2 X i hA y i (t)A i (t)idt ! j (t)i +O(dt 2 ) : (2.23c) If a jump of type i occurs, we have j (t +dt)i = A i (t)j (t)i kA i (t)j (t)ik = A i (t)j (t)i q hA y i (t)A i (t)i : (2.24) Therefore, in analogy to Eq. (2.22) we can write the stochastic Schr odinger equation for the normalized state as in Eq. (2.11). 2.4 Simulation procedure for adiabatic quantum trajectories In this section we formulate an algorithm for implementing adiabatic quantum trajectories. We start by noticing that the update in Eq. (2.23a) corresponds to the evolution by the rst part of the stochastic Schr odinger equation. Therefore, the deterministic evolution in the rst term of Eq. (2.11) is equivalent to propagating the state vector via the Schr odinger equation with H e (t) and then renormalizing it. When a jump occurs, one of the operators A i (t) is applied. The relative weight of each A i (t) is dp i (t), given in Eq. (2.20b). In this case, the state is evolved as in Eq. (2.21), and the normalized state is given in Eq. (2.24). The update in Eq. (2.24) corresponds to the evolution by the second term of Eq. (2.11). 16 This provides a direct way to algorithmically implement the quantum trajectories method. Starting from a known normalized initial state, the state is evolved via a sequence of deter- ministic evolutions and jumps, as in Eqs. (2.18) and (2.24), by drawing a random number at each nite but small time-step t and determining which of the two choices to take. Compared to the standard time-independent case, the size of the time-step must satisfy additional conditions in order for the approximations to hold: t min t 2kH e (t)k k _ H e (t)k ; 1 kH e (t)k ; (t) 2 (t) _ (t) ; (2.25) wherekk is the operator norm (largest singular value). We give a proof of this new bound in Appendix A.2. While the second and third terms reduce to the known conditions for the time-independent case, the rst term in Eq. (2.25) is unique to the time-dependent case and re ects the error associated with time evolution under a time-dependent eective Hamiltonian. This term highlights the fact that the faster the eective Hamiltonian and its eigenstates vary in time, the smaller is the time-step required to properly follow the trajectory. Eq. (2.25) can also be viewed as the physical timestep upper bound in weak measurement. However, drawing a random number at each time-step is computationally expensive, so it is more ecient to use the waiting time distribution [29] to determine the rst jump event. As we mentioned before [Eq. (2.19)] the square norm of the unnormalized wavefunction at t +dt gives the probability of no jump during the innitesimal interval [t;t +dt]. We show in Appendix A.3 that starting from the normalized statej (t)i, the probability of no jump occurring in the nite (not necessarily small) time interval [t;t +] is given by k ~ (t +)k 2 = exp Z t+ t (s)ds ; (2.26) where the jump rate (t) is given in Eq. (2.20a). With this, the simulation procedure for one single trajectory is as follows, starting from t: 17 • Draw a random number r. • Propagate the unnormalized wavefunction by solving the Schr odinger equation with H e [Eq. (2.17)] until the jump condition is reached at t +, i.e., for such that h ~ (t +)j ~ (t +)ir. (Recall that the norm of the unnormalized wavefunction will keep decreasing in this process.) • Determine which jump occurs by drawing another random number and update the wavefunction by applying jump operators, and renormalize. • Repeat the above steps with the new normalized state. • Repeat until the nal simulation time is reached. We prove that averaging over quantum trajectories recovers the master equation in Ap- pendix A.2. Specically, we show there that if we denote the state of the k-th trajec- tory at time t byj k (t)i, then we can approximate the master equation solution for the density matrix (t) as 1 n P n k=1 j k (t)ih k (t)j for large n. Choosing a basisfjz i ig for the system Hilbert space, we can thus approximate the density matrix elementhz i j(t)jz j i as 1 n P n k=1 hz i j k (t)ih k (t)jz j i for large n. 2.5 V: Case studies We consider a system of N qubits with a transverse-eld Ising Hamiltonian given by H S (t) = A(t)H X S +B(t)H Z S ; (2.27a) H X S N X i=1 x i ; (2.27b) H Z S N X i=1 h i z i + N X i>j=1 J ij z i z j : (2.27c) 18 J=-1 J=-1 J=-1 J=-1 J=-1 J=-1 J=-1 h 1 = 1 4 (a) -2/3 -2/3 2/3 -1 1/3 1 -1 1 (b) h L h R =−1 j=1 j=5 j=2 j=6 j=3 j=7 j=4 j=8 J=−1 (c) Figure 2.2: Graphs of (a) the 8-qubit chain, (b) the 8-qubit Hamiltonian exhibiting a small gap, and (c) the 16-qubit \tunneling-probe" Hamiltonian of Ref. [47]. (a) Only the rst qubit is subjected to an applied eld and each qubit is ferromagnetically coupled with J = 1. (b) Solid lines corresponds to ferromagnetic coupling and dashed lines corresponds to antiferromagnetic coupling. The thickness denotes the strength of the coupling. Local elds are shown inside the circles. Full parameters are given in Eq. (2.32). (c) The left 8-qubit cell and right 8-qubit cell are each subjected to applied elds with opposite direction. Each qubit is ferromagnetically coupled to others as shown by the lines, with J =1. We assume that the qubit-system is coupled to independent, identical bosonic baths, with the bath and interaction Hamiltonian being H B = N X i=1 1 X k=1 ! k b y k;i b k;i ; (2.28a) H I =g N X i=1 z i X k b y k;i +b k;i ; (2.28b) 19 0 0.2 0.4 0.6 0.8 1 0.5 0.6 0.7 0.8 0.9 1 Instantaneous GS Population QT ME 10 2 10 3 10 4 0.3 0.4 0.5 0.6 Figure 2.3: The evolution of the population in the instantaneous ground state for the 8-qubit problem in Eqs. (2.27) and (2.30) for a total time of t f = 10s and temperature 2:6GHz (in ~ 1 units), as a function of the normalized time s = t=t f . The quantum trajectories results with 10 4 trajectories (QT, blue circles) are in excellent agreement with the adiabatic master equation (ME, red solid line). Inset: the convergence of the ground state population (averaged over quantum trajectories) towards the master equation result as a function of the number of trajectories, at s = 1. The error bars represent 2 condence intervals, where is the standard deviation of the mean generated by taking 10 3 bootstraps over the number of trajectories. where b y k;i and b k;i are, respectively, raising and lowering operators for the k-th oscillator mode with natural frequency ! k . The bath correlation functions appearing in Eq. (2.2) are given by ij (!) = 2g 2 !e j!j=!c 1e ! ij ; (2.29) arising from an Ohmic spectral density [13], and satisfy the KMS condition (2.3). 20 2.5.1 8-qubit chain As a rst illustrative example and as a consistency check, we reproduce the master equa- tion evolution of the 8-qubit ferromagnetic Ising spin chain in a transverse eld studied in Ref. [13]. For this problem, the Ising parameters are given by [also shown in Fig. 2.2(a)] h 1 = 1 4 ; h i>1 = 0 ; J i;i+1 =1 ; i = 1;:::; 8 : (2.30) The functional form of the functions A(t) and B(t) is given in Appendix A.4 (they are the annealing schedule of the D-Wave One \Rainier" processor, described in detail, e.g., in Ref. [48], and were also used in the original AME study of the 8-qubit chain [13]). The eect of changing the schedule is small; for comparison we provide in Appendix A.4 results for a linear annealing schedule with the same bath parameters. As shown in Fig. 2.3, we recover the master equation solution within the error bars. The initial state is the ground state of H S (0), which is the uniform superposition state. It is illustrative to see how a single trajectory diers from the averaged case, and we show this in Fig. 2.4. Instead of the smooth change in the population as observed in the averaged case, the single trajectory behaves like a step-function. This is explained by the fact that the drift term vanishes ifj (t)i is a nondegenerate eigenstate, as shown in Appendix A.5. Therefore changes in the state's overlap with the instantaneous ground state occur only due to the jump operators. In this picture, the ground state population revival observed after the minimum gap is crossed is associated with jumps from the rst excited state (or higher states for largeT ) to the ground state. After the minimum gap, there are more transitions back to the ground state than out of the ground state (see the inset of Fig. 2.4). Such a dierence (divided by the number of trajectories) leads to the rise of the ground state population. 21 Figure 2.4: The overlap squared of the (normalized) state with the instantaneous ground state ofH S (t) for a typical single trajectory of the 8 qubit chain in Sec. 2.5.1, witht f = 10s and temperature 2:62GHz, as a function of the normalized time s = t=t f . The sudden changes in the overlap are due to the action of the jump operatorsfA i (t)g, taking the state from one eigenstate to another. This is to be contrasted with the smooth behavior of Fig. 2.3 when we average over dierent trajectories. Inset: a histogram of the net number of jumps to the instantaneous ground state (GS). A negative number indicates a jump out of the ground state, and a positive number indicates a jump towards the ground state. The change from negative to positive net jumps occurs at the minimum gap point. Using Eq. (2.20), we can give an explicit expression for the jump rate from the rst excited state back to ground state. As shown in Appendix A.6, this is given by: 1!0 (t) = 8 X =1 (! 10 )jh 0 (t)j z j 1 (t)ij 2 : (2.31) 22 0 0.2 0.4 0.6 0.8 1 0 0.2 0.4 0.6 0.8 1 Instantaneous GS Population QT ME Closed 0 0.5 1 10 -2 10 0 10 2 Instantaneous gap Figure 2.5: The evolution of the population in the instantaneous ground state for the 8-qubit problem in Eqs. (2.27) and (2.32) for a total time of t f = 10s and temperature 1:57GHz, as a function of the normalized times =t=t f . The quantum trajectories results with 5 10 3 trajectories (blue circles) are in excellent agreement with the master equation (red solid line). We also show the closed-system evolution (yellow dashed line) to highlight that the evolution is not adiabatic. Inset: the instantaneous energy gap between the rst excited state and the ground state during the anneal. The minimum occurs at s = 0:46, coinciding with the sharp discontinuity observed in the instantaneous ground state population. 2.5.2 8-qubit non-adiabatic example We now consider an 8-qubit problem with a suciently small minimum gap such that the closed-system evolution is not adiabatic even witht f = 10s and using the DW2X annealing schedule (described in detail, e.g., in Ref. [40]). While this strictly violates the assumptions under which the AME is derived, 2 we can ask about the dynamics associated with the master 2 Equation (27) in Ref. [13] is a necessary condition for the validity of the AME. It states that h 2 t f 1, where is the ground state gap and h = max s2[0;1];a;b jh a (s)jj@ s H S (s)j b (s)ij, with s = t=t f andj a (s)i being the instantaneous a-th eigenstate of the system Hamiltonian H S (s). We nd that for t f = 10s, the l.h.s. 5. 23 equation irrespective of its origins. We are interested in this example since it illustrates some aspects of the quantum trajectories picture which are not visible in the adiabatic limit, as explained below. The Ising Hamiltonian H Z S is dened with parameters: 3 ~ h = (2;2; 2;3; 1; 3;3; 3) ; (2.32a) 3J 0;4 = 1; 3J 0;5 =3; 3J 0;6 =1; 3J 0;7 =1; (2.32b) 3J 1;4 =1; 3J 1;5 =2; 3J 1;6 = 2; 3J 1;7 =3; (2.32c) 3J 2;4 =2; 3J 2;5 =1; 3J 2;6 =3; 3J 2;7 =2; (2.32d) 3J 3;4 =3; 3J 3;5 =3; 3J 3;6 =1; 3J 3;7 =3; (2.32e) Figure 2.5 shows our simulation results, obtained by solving the AME directly and by using the trajectories approach. Reassuringly, the agreement between the two is excellent. Also plotted are the closed system results for this problem, which exhibit a sharp diabatic transi- tion out of the ground state at the minimum gap point (the small gap is shown in the inset). The AME and trajectories results show that the ground state population loss starts before the diabatic transition, due to thermal excitations, but that the ground state population loss is partially mitigated by the presence of the thermal bath, with the open system ending up with a higher ground state population than the closed system. The diabatic transition results in dierent trajectories than those observed for the adia- batic case in Sec. 2.5.1. We show such a case in Fig. 2.6. Instead of the pulse-like structure seen in Fig. 2.4, we observe a combination of both drifts and jumps. Because the diabatic transition generates a non-eigenstate that is a coherent superposition of the ground state and rst excited state, drifts caused by the environment show up in the subsequent evolu- tion. Furthermore, this superposition also means that the Lindblad operator associated with ! = 0, if having dierent component weights in Eq. (2.4), can also induce jumps (e.g., the jumps around s = 0:6 in Fig. 2.6), an eect that is completely absent in the adiabatic case. 24 These jumps need not project the state completely onto an instantaneous energy eigenbasis state, but they can change the relative weights on the dierent occupied eigenstates, which manifest themselves as `incomplete' jumps in the trajectories. We reemphasize that due to the violation of the conditions under which the AME is derived, the observations we have reported for this example are strictly valid only when the AME is taken at face value, and do not necessarily re ect actual physical dynamics. 0 0.2 0.4 0.6 0.8 1 0 0.2 0.4 0.6 0.8 1 Population Trajectory, GS Trajectory, ESs Figure 2.6: The overlap squared of the (normalized) state with the instantaneous ground state of H S (t) (blue solid curve) and the sum of the rst three instantaneous excited states (red dashed curve) for a typical single trajectory of the 8-qubit problem in Sec. 2.5.2, with t f = 10s and temperature 1:57GHz, as a function of the normalized time s = t=t f . A diabatic transition occurs at s = 0:46. The continuous decay immediately afterward and between the `incomplete' jumps is to be contrasted with the step-function like single trajec- tories of the adiabatic case seen in Fig. 2.4. 25 0 0.2 0.4 0.6 0.8 1 0 0.2 0.4 0.6 0.8 1 Population Figure 2.7: Quantum trajectory results for a temperature of 12mK and 20mK using the DW2X annealing schedule. Results are averaged over 5000 trajectories. A revival of in- stantaneous ground state population occurs after the minimum gap (shown by the dashed vertical line). 2.5.3 16-qubit \tunneling-probe" Hamiltonian In order to demonstrate the computational utility of the trajectories approach over the master equation approach, we now give results for a 16-qubit system rst studied in Ref. [47] for the purpose of probing tunneling in quantum annealing. For this problem, the parameters of Eq. (2.27) are [also shown in Fig. 2.2(b)]: h L = 0:44 ; h R =1 ; J i;i+1 =1 ; i = 1;:::; 16 : (2.33) where the setsL andR range overi = 1;:::; 8 andi = 9;:::; 16, respectively. Ref. [47] chose the value ofh L to ensure that the minimum ground state gap is lower than the temperature 26 T = 15:5mK, in order to study a non-trivial interplay between tunneling and thermal ac- tivation. 3 This parameter choice means that incoherent eects play a relatively strong role in this problem, which are not well captured by the AME. Thus, similarly to the previous example, the AME is being used here outside of its strict validity domain. We are inter- ested in testing whether it can nevertheless qualitatively capture the correct physical eects. Moreover, direct master equation simulations for such a large system take longer than 24 hours (which is a standard time-window on high-performance clusters), while each quantum trajectory takes less than 24 hours. We can then exploit many CPU cores to perform many trajectories in parallel. To this end we used 320 CPU cores and repeated the simulation 16 times for a total of over 5000 trajectories. Our simulations (see Fig. 2.7) show how population is lost from the instantaneous ground state to the rst excited state near the minimum gap point s 0:308. It also shows a small population revival after the minimum gap is crossed. As in Sec. 2.5.1, this revival is associated with jumps from the rst excited state (or higher states for largeT ) back to the ground state. Encouragingly, despite the perturbative nature of the AME, this revival is qualitatively in agreement with the results of Ref. [47] (see their Fig. 4). The latter work found a stronger revival on the basis of the non-perturbative, non-interacting blip approximation (NIBA), which more accurately captures additional transitions that occur when the energy level broadening is larger than the energy gap between energy levels. 2.6 Simulation cost comparison We now provide a cost comparison between the simulations cost of directly solving the AME and the quantum trajectories method. The rst two subsections in this section follow Ref. [29] closely (while adding some details), and we borrow the notation used in that reference. In 3 The problem features one global minimum and one local (false) minimum, which are separated by a tall energy barrier at around the minimum gap point (see Fig. 2 in Ref. [47]); to reach the global minimum from the false minimum, the system state has to transverse the barrier. Such transitions can be modeled as quantum jumps in the quantum jump picture. 27 subsection 2.6.3 we provide a new analysis that reveals that the quantum trajectories method can exhibit a scaling advantage ranging between O(N) and O(N 2 ) over the direct solution of the AME. 2.6.1 Number of trajectories The number of trajectories needed can be found from the standard error of the sample mean. As an example, let us consider the standard error ^ t associated with the instantaneous ground state populationh 0 (t)j(t)j 0 (t)i: ^ 2 t = 1 R(R 1) R X r=1 jh 0 (t)j r (t)ij 2 ^ M t 2 ; (2.34) wherej r (t)i denotes the state associated with trajectoryr at timet, ^ M t = 1 R P R r=1 jh 0 (t)j r (t)ij 2 , and R is the total number of trajectories. By xing the value of the standard error ^ t , the number of necessary trajectories R can then be determined. 2.6.2 Cost comparison Since we expect ^ t 1 p R , let us write ^ 2 t = B (N) R 1 R ; (2.35) where B (N) = 1 R1 P R r=1 (h r (t)jBj r (t)i ^ M t ) 2 for an observable B and mean value ^ M t = 1 R P R r=1 h r (t)jBj r (t)i. The factor B (N) is a non-increasing function of the sys- tem dimension N [29]: B (N)N x ; (2.36) where the scaling x depends on the observable: 0 (not self-averaging)x 1 (strongly self-averaging) : (2.37) 28 Thus, to obtain the same standard error for increasing dimension, the number of trajectories need not be increased in general. This is another advantage of the trajectories method for growing system dimension. Such a phenomenon has also been observed in time-dependent stochastic density functional theory [49]. From Ref. [29], the total serial CPU time required for the simulation of the master equation, denotedT AME , versus the stochastic method with R trajectories, denoted T StS , is: T AME =k 1 s 1 (N)N ; (2.38a) T StS =k 2 R(N)s 2 (N)N ; (2.38b) where k 1 and k 2 are constants depending on the specic implementation of each method, s 1 (N) is the total number of evaluations ofL WCL [(t)] [Eq. (2.2)] using the master equation method, ands 2 (N) is the total number of evaluations ofH e (t)j (t)i [Eq. (2.10)] in a single trajectory. R(N) in Eq. (2.38b) is the minimum number of trajectories needed to obtain a standard error lower than a particular chosen value. To account for the constraint that R(N) 1, we rewrite Eq. (2.36) as B (N) = B N x , and Eq. (2.35) as: R(N) = 8 > > < > > : l B N x ^ 2 t m for N <N 1 for NN (2.39) where N =d( B =^ 2 t ) 1 x e. For x > 0, the required number of trajectories decreases with N untilN , after which one trajectory gives the expectation value within the desired accuracy. 29 In general, the number of operations needed to evaluateL WCL [(t)] relative to the number needed to evaluate H e (t)j (t)i diers by a factor of N, so that + 1, and Eq. (2.38) becomes T AME =k 1 s 1 (N)N +1 ; (2.40a) T StS = 8 > < > : k 0 2 s 2 (N)N x for N <N k 2 s 2 (N)N for NN (2.40b) where k 0 2 =k 2 B ^ 2 t ; (2.41) andk 0 2 hence depends on the required accuracy as well. In many situation s 1 (N) ands 2 (N) grow withN, but they are roughly equal or grow in same manner withN. By dividing these two expressions, we can obtain the ratio of T AME =T StS , T AME T StS = 8 > > < > > : k 1 k 0 2 N 1+x for N <N k 1 k 2 N for NN : (2.42) Since 0x 1, we can write k 1 k 0 2 N T AME T StS k 1 k 0 2 N 2 for N <N ; T AME T StS = k 1 k 2 N for NN : (2.43) We say that the trajectories method has an advantage over the direct master equation solution if T AME T StS > 1. The constant factor k 1 k 0 2 is typically a small number because it is proportional to the required standard error squared [Eq. (2.41)]. Therefore there is an advantage for the trajectories method when eitherN is suciently large or when a sucient number of CPU cores C is available (see the next subsection). Equation (2.43) shows that an advantage beyond linear in N is attainable for N <N on a single CPU. The reason is 30 that the number of trajectories needed to achieve a xed accuracy decreases with increasing system dimension. ForN >N , only one trajectory is required, and the advantage scales as O(N). We note that the larger-than-linear advantage only holds if we are interested in estimating operators with the same self-averaging property. This is in contrast to evolving the entire density matrix, as in the AME, which allows the expectation value of any observable to be calculated. If we demand this same capability from the trajectories approach, then only the linear advantage holds. 2.6.3 Parallel implementation The stochastic wave function method is very well-suited for parallel computing implemen- tations. The communication needed between each core is minimal since each trajectory is simulated independently. AssumingC CPU cores are used, whereCR(N), we can adjust the time cost [Eq. (2.38b)] for the stochastic method to T StS =k 2 R(N) C s 2 (N)N : (2.44) Note that the number of coresC is held constant, i.e., is independent of the system dimension N. Therefore, we can update Eq. (2.43) to: k 1 k 0 2 CN T AME T StS k 1 k 0 2 CN 2 for N <N ? ; T AME T StS = k 1 k 2 N for NN ? : (2.45) where N ? =d( B = (C^ 2 t )) 1 x e. Here N ? is the system dimension where R(N ? ) =C, and one execution of the C parallel CPU cores is enough to obtain the desired standard error. 31 Again, the larger-than-linear advantage inN only holds if we are interested in estimating operators with the same self-averaging property, and otherwise we can only expect a linear advantage in N. 2.7 Conclusion of this chapter In this chapter, we have shown how quantum trajectories (in the form of quantum jumps) can be unravelled from the adiabatic master equation. We have described and demonstrated a simulation procedure in terms of the waiting time distribution that reproduces the re- sults of the master equation for examples involving 8 and 16-qubit systems. Direct master equation simulations for the 16-qubit example would take a long time, but the simulation of the quantum trajectories remains computationally feasible for larger system dimensions by allowing us to simulate many trajectories in parallel. A scaling cost comparison of the two methods shows that, generically, the quantum trajectories method yields an improve- ment by a factor linear in the system dimensionN over directly solving the adiabatic master equation. However, the trajectories method can be expected to be up to a factor cN 2 faster than a direct simulation of the master equation if only the expectation value of specic self- averaging observables is desired. Here c is a constant proportional to the number of parallel processes and the target standard error. We therefore believe this approach will be particularly useful in enabling the study of larger systems than has been possible using a direct simulation of the AME. In addition, the quantum trajectories method oers fresh physical insight into the nature of individual trajectories and their statistics, which may become a helpful tool in interpreting computational bottlenecks in quantum annealing and adiabatic quantum computing. Finally, while we did not address this in the present chapter, the quantum trajectories approach is well known to be a convenient path towards continuous measurement and the inclusion of quantum feedback control [28]. This approach might in the future provide a 32 path towards error correction of adiabatic quantum computing, e.g., by formulating control targets that push the system back to the ground state after diabatic or thermal transitions. We address this in Chap. 4. 2.8 Acknowledgments of this chapter The research is based upon work (partially) supported by the Oce of the Director of Na- tional Intelligence (ODNI), Intelligence Advanced Research Projects Activity (IARPA), via the U.S. Army Research Oce contract W911NF-17-C-0050. The views and conclusions contained herein are those of the authors and should not be interpreted as necessarily rep- resenting the ocial policies or endorsements, either expressed or implied, of the ODNI, IARPA, or the U.S. Government. The U.S. Government is authorized to reproduce and distribute reprints for Governmental purposes notwithstanding any copyright annotation thereon. Computation for the work described in this chapter was supported by the Uni- versity of Southern California's Center for High-Performance Computing. The codes in this chapter are available at: https://github.com/USCqserver/Adiabatic-master-equation-and-quantum-trajectories 33 Chapter 3 1=f noise 3.1 Introduction The trajectories method in the previous Chapter can be applied to non-Markovian noise, in particular 1=f (telegraph) noise. 1=f noise means that the spectrum of the environment S(f) / 1=f. To model 1=f (telegraph) noise, one needs to add stochastic elements to trajectories simulations and perform ensemble averaging. We start with the known Bloch vector-magnetic eld approach [50{53], which in many ways is similar to spin vector dynam- ics [54, 55]. Then we formulate a stochastic Hamiltonian approach for 1=f noise. The latter approach allows an easier adaption to the time-dependent annealing Hamiltonian, and also to high energy-level modeling of superconducting qubits. There are also experimental studies of 1=f noise on superconducting qubits [56, 57]. Note that in [56], a stochastic Hamiltonian has been considered, but only one stochastic Hamiltonian is included in the model. In this thesis we formulate a stochastic Hamiltonian approach for 1=f noise, that includes a series of uctuator terms with dierent ipping frequencies and couplings, which can be further added to the annealing system Hamiltonian. The new and original contributions of this chapter include the followings: 34 • We derive from the Bloch vector-magnetic eld approach a stochastic Hamiltonian approach that includes a series of uctuator terms with dierent ipping frequencies and Gaussian couplings that altogether constitute the 1=f power spectrum; • We, for the rst time, show how to incorporate the many uctuators into the inher- ently time-dependent annealing Hamiltonians, and construct an overall time-dependent stochastic Hamiltonian. We thoroughly analyze how noise with 1=f power spectrum aects quantum annealing; • We introduce exibility into the description of the noise uctuators. For example, we allow the noise uctuators to take any form of multi-level operators and any noise-axis direction, depending on the superconducting qubits and circuit decoherence process. 3.2 1=f noise: Bloch vector-Magnetic eld approach 3.2.1 One qubit and one uctuator Figure 3.1: Bloch vector representing the state of a qubit in the rotating (with angular frequency B 0 ) frame of reference. In the laboratory frame of reference it precesses around the z-axis with (time-dependent) angular velocity B. Fluctuations of this velocity is b(t) Photo credit: [50]. 35 3.2.1.1 Qubit One can express the qubit under Hamiltonian drive in the form of a Bloch vector M = (M x ;M y ;M z ) under the in uence of a sum of magnetic elds B. The dynamics of the qubit is described by the precession equation (Schrodinger equation) _ M = B M (3.1) 3.2.1.2 Fluctuator Assume for simplicity that the qubit is subjected to a time-independent magnetic eld point- ing in the z-direction B = B = B z . However, a stochastic magnetic eld b(t) exists and is pointing parallel toB such that the total magnetic eld isB +b(t). In this section we focus on one stochastic magnetic eld b(t). The stochastic magnetic eld b(t) originates from a so-called uctuator, which can nd itself in one of two (meta)stable states, 1 and1, and once in a while makes a switch between them. We consider the symmetric random telegraph (RT) process: 1!1 = 1!1 = =2. When the uctuator is in state 1, the stochastic magnetic eld is b(t) = +b; similarly when it is in state1, the stochastic magnetic eld is b(t) =b. Dene m + (t) = M x (t) +iM y (t). Its ensemble averagehm + (t)i is called free induction decay (FID) [50]. The decay constant of the ensemble averagehm + (t)i gives T 2 . The dierential equation satised byhm + (t)i and its analytical solution is given in [50] as: h m + i + h _ m + i =b 2 hm + i (3.2) 36 with initial conditionshm + (0)i = 1,h _ m + (0)i = 0 (The initial state of the qubit is (j0i + j1i)= p 2, such that the initial Bloch vector pointing in the x-direction: M(0) = (1; 0; 0)). The solution of Eq. (3.2) with these initial conditions is hm + i = (2) 1 e t=2 ( + 1)e t=2 + ( 1)e t=2 ; p 1 (2b= ) 2 : (3.3) 3.2.1.3 Comparison of analytical solutions and stochastic simulations We develop new and ecient parallelizable simulation techniques based on the algorithm developed in the previous chapter on quantum trajectories (Chap. 2). Using this algorithm, one does not need to draw a random number at each timestep, but instead calculate the waiting times for each ipping of the stochastic b-eld. We compare the analytical solution of Eq. (3.3) with the stochastic random eld process simulation results we produced, which are averaged over many trajectories. The simulation procedure is detailed next. Dene the ratio of the strength of the uctuator over the rate of the uctuator as: g =b=( =2) = 2b : (3.4) Using Gillispie's algorithm [58] in the previous chapter of quantum trajectories (Chap. 2) to determine the ipping (waiting) times of uctuator (thus the sign of the stochastic magnetic eld), and averaging over many realizations, we plot the time evolution of the ensemble averagehm + (t)i withg = 0:8 andg = 5. Although in each realization the initial state of the qubit is xed, the initial state of the uctuator is not. If the uctuator starts in either state with equal probability, the initial conditionh _ m + (0)i = 0 above is fullled. We numerically simulatehm + (t)is with dierent g and plot the results as blue lines in Fig. 3.2. In the same gure, we also plot the analytical solutions of Eq. (3.3) as green lines. Note that the unit of is (number of switches)/(time unit) [51]. 37 Figure 3.2: Blue lines: results reproduced from averaging 8000 realizations. The initial state of the uctuator is either state 1 or1 with equal probability. Green lines: analytical solutions. Note that b = 5 for both curves. We perform simulations for a wider range of g and plot the results in Fig. 3.3. We have several observations. The resulting ensemble averagehm + i of the qubit decays faster (or T 2 is shorter) when g is larger (the stronger the uctuator is or the slower the uctuator switches). However, its decay rate saturates as the ratio g reaches 1, where the behavior of hm + (t)i transitions from strictly decay to damped oscillation. 38 Figure 3.3: Saturation of decay as b further increases. Note that = 1 switch/(time unit) for all curves. The decay rate is saturated (cannot decay faster) after g = 2b= reaches 1. The curves are produced from averaging 10000 realizations with a random initial uctuator state for each realization. 3.2.2 One qubit and an ensemble of uctuators 3.2.2.1 An ensemble of uctuators to generate an overall 1=! power spectrum: S(!) =A=!. The following terminology is summarized from [50{53]. A single two-state uctuator is described by a stochastic variable (t) =1 with rates 1!2 = 2!1 = =2. The power spectrum of the noise generated by the i-th uctuator is: S i (!) =b 2 i L i (!); (3.5) 39 whereL i (!) is a Lorentzian function of frequency, L i (!) 1 i ! 2 + 2 i : (3.6) The couplings b i are distributed with a small dispersion around an average valuehbi. If a distribution P ( )/ 1= is assumed for the switching rates i 2 [ min ; max ], the total uctuation after summing up all the uctuating elds (t) = P i i (t) exhibits a 1/! total power spectrum of the form S(!) =A=!,A > 0. In this case, let the constant of proportionality in the distribution function P ( ) be P ( ) = P 0 2 : (3.7) Under these conditions the total power spectrum S(!) reads S(!) =hb 2 i Z max min d P 0 2 L (!) A ! : (3.8) The amplitudeA can be expressed in terms of the number of uctuators per noise decade, n d =N T ln (10)= ln ( max = min ), as followsA =hb 2 iP 0 =4 =hb 2 in d =(2 ln(10)). Other notes: • n d means the number of uctuators per decade. • The strength b i of the uctuators has Gaussian distribution with average ashbi, and standard deviation b. • The initial value of each uctuator is randomly sampled according to thermal equilib- rium distribution p eq , withhp 0 i> 0 means more uctuators start at state 1 than in 1. 40 3.2.2.2 Distributions of the uctuators Fig. 3.4(a) illustrates the couplings between the qubit and uctuators with strength b i and ipping frequency i . The (1= ) distribution function (Eq. (3.7)) of the frequencies of the uctuators corresponds to a uniform distribution of ln . Fig. 3.4(b) shows the frequency distribution of 12000 uctuators, where there are 12 decades of noise ([10 0 ; 10 12 ] Hz) and the number of uctuators per decade is n d = 100. 1 2 n b n b 2 b 1 Qubit Fluctuators (a) 0 200 400 600 800 1000 0 2000 4000 6000 8000 10000 (b) Figure 3.4: (a) A cartoon illustrating a single qubit coupled to a series of uctuators with dierent individual ipping frequencies (b) Distribution of frequencies of 12000 uctuators. 3.2.2.3 Simulation results: The importance of the number of decades and min Unlike the case of a single uctuator, there are no analytical solutions for an ensemble of uctuators. We have to perform numerical simulations. Similar to the case of a single uctuator, we use Gillispie's algorithm to determine the ipping times of uctuators. At each iteration, we randomly select a uctuator based on the uniform distribution of ln 41 (Sec. 3.2.2.2). Again one can measure how fast the coherence is lost by how fast the FID hm + (t)i decays (or oscillates). Due to the 1= distribution function, most of the uctuators have their ipping frequency close to min . By xing max , the more decades of noise, the smaller the min . We perform simulations to investigate the sensitivity of qubit coherence to the number of decades of uctuator frequencies (and thus min ) and plot the results in Fig. 3.5. For example, if max = 10 THz, 2 decades of noise means min = 100 GHz and 4 decades of noise means min = 1 GHz. As we can see in Fig. 3.5, the decay ofhm + (t)i is faster with more decades of noise (and thus a smaller min ). 0 1 2 3 4 5 10 -9 0 0.2 0.4 0.6 0.8 1 Figure 3.5: Reproduced gures of 2 to 5 decades. Couplings distributed withhbi=hbi = 0:2 and aroundhbi = 4:6 10 7 Hz. max = 10 THz with dierent min . n d = 1000. dp 0j =1 are distributed according tohdp 0j i =dp eq = 0:08. 3.2.2.4 Simulation results: Dependence on b=hbi We plot in Fig. 3.6 another case where the decay time is sensitive to the average strength of uctuators min . In general, the stronger the uctuators are the fasterhm + (t)i decays. We observe two interesting points in Fig. 3.6. Comparing the top subgure with the bottom subgure, the bottom subgure has b=hbi increased 100 times. 42 • The decay in the bottom subgure is 100 times faster. • to obtain a similar decay envelope as in the top subgure, min also has to increase 100 times. 0 5000 10000 0 0.5 1 h m + ( t ) i 0 50 100 t i m e ( s ) 0 0.5 1 h m + ( t ) i Figure 3.6: Top: min = 10 6 GHz, max = 100 GHz, n d = 100,hbi = 10 4 =4; Bottom: min = 10 4 GHz, max = 100 GHz, n d = 100,hbi = 0:01=4. Results are averaged over 8k realizations. Each uctuator is initially in a pure state randomly sampled according to hp 0 i =p eq = 0:08. b=hbi 0:2. 3.2.2.5 Simulation results: Non-exponential decay We have seen in Fig. 3.5 that the more decades of noise are included, the more rapid the coherence is lost. For exponential decay, lnhm + (t)i would decrease linearly. We can dene the reduced Decay factor of the qubit as (t) = lnhm + (t)i [50]. We plot in Fig. 3.7 the eect of the 12000 uctuators on the (t) evolution with frequencies distributed according to Fig. 3.4(b). As seen in Fig. 3.7, the reduced decay factor does not decrease in a linear manner. As time doubles from 1ns to 2ns, the reduced decay factor drops more than 2 times. 43 0 0.5 1 1.5 2 t ( s e c ) #10 -9 -3 -2.5 -2 -1.5 -1 -0.5 0 ! ( t ) R e d u c e d D e c a y F ac t o r ! ( t ) Figure 3.7: 12 decades of noise: max = 10 12 Hz with min = 10 0 = 1 Hz. Couplings distributed withhbi=hbi = 0:2 and aroundhbi = (4:6=) 10 7 Hz. n d = 1000. dp 0j =1 are distributed according tohdp 0j i =dp eq . As time doubles, the reduced decay factor drops more than 2 times. 3.3 1=f noise: Stochastic Hamiltonian formulation We consider the stochastic Hamiltonian formulation of 1=f noise, which can also be applied to the case where the Hamiltonian drive of the qubit is time-dependent. A Hamiltonian is only dened up to a unitary transformation (change of basis). Therefore we can dene the system Hamiltonian with x and z components only, and with real A(t) and B(t) (for example, they can be the DWave NASA schedule in Fig. 1.1): ^ H sys = 1 2 (A(t) x +B(t) z ) : (3.9) 44 It can be cast into the spin-(time dependent) magnetic eld description in Sec. 3.2 as ^ H sys = (1=2)B; B z B(t); B x A(t): (3.10) The Schr odinger equation with this system Hamiltonian drive turns out to be equivalent to the precession equation for the Bloch vector: _ M = B M: (3.11) 3.3.1 Single uctuator A single uctuator induces an extra stochastic magnetic eld b 1 (t) = b 1 1 (t) along the B z axis. Its eect can be included in the ^ H sys , as: ^ H sys = 1 2 (A(t) x +B(t) z ) + 1 2 b 1 1 (t) z ; (3.12) where 1 (t) switches with rate 1 =2. The 1 2 convention follows from Sec. 3.2.2.1. 3.3.2 Many uctuators For many uctuators, ^ H sys = 1 2 (A(t) x +B(t) z ) + 1 2 X i b i i (t) z ; (3.13) where i (t) switches with rate i =2. The ipping frequencies i (t) of the uctuator-i is randomly drawn from the distribution of 1= i . The initial thermal distribution of i (t), and Gaussian distribution ofb i follows from Sec. 3.2.2.1. 45 3.3.3 Case studies 3.3.3.1 Single uctuator: A(t) = 0 and B(t) = 0 for all t Consider the following system Hamiltonian, where the system drive is inherently zero (A(t) = 0 and B(t) = 0). Therefore the system Hamiltonian only has a single uctuator pointing along the z-axis. ^ H sys =b 1 1 (t) z ; (3.14) is the case where there is no system driving. (where 1 (t) switches with 1 =2. ) The qubit's initial state isj+i = j0i+j1i p 2 . Denotej (t)i =c(t)j0i +d(t)j1i the qubit state at timet. With this initial condition, the Bloch vector can be expressed in terms of thej+i's coecients as: M x (t) =cd +c d; (3.15) M y (t) =icd ic d: (3.16) Denotej (t)i the state vector obtained in one trajectory of the stochastic Schr odinger equation: i~ d dt j (t)i = ^ H sys (t)j (t)i: (3.17) After averaging many trajectories, we recover the density matrix evolution. Dene (t) as (t) =hj (t)ih (t)ji = 1 2 0 B @ 1 +hM z (t)i hM x (t)iihM y (t)i hM x (t)i +ihM y (t)i 1hM z (t)i 1 C A = 1 2 0 B @ 1 +hM z (t)i hM x (t)iihM y (t)i hm + (t)i 1hM z (t)i 1 C A : (3.18) The o-diagonal element of the density matrix coincides with the denition ofhm + (t)i in Sec. 3.2.1.2. We plot in Fig. 3.8 the evolution the o-diagonal element 01 (t)=2 of the density 46 matrix, by performing a trajectories simulation of the stochastic Schr odinger equation, for two values of g = 2b 1 = 1 . Note that in this time-independent drive (Eq. (3.14)) the 01 (t)=2 is a real number (hM y (t)i = 0). We see in Fig. 3.8 that the simulation data agrees with Fig. 3.3 obtained by the Bloch vector-Magnetic eld approach. 0 2 4 6 8 -1 -0.5 0 0.5 1 Figure 3.8: g = 2b= . Averaged over 10000 trajectories. 3.3.3.2 Multiple uctuators with a time-dependent system Hamiltonian Consider the following example. ^ H sys = A 2 1 t t f x + t t f z + 1 2 10 X i=1 b i i (t) z : (3.19) where: A = 1 GHz; The sum over uctuator-i with the uctuator i (t) distributed in 1= i with i 2 [ min ; max ]; min = 0:1 GHz, max = 10 GHz. i (t) switches with i =2 and takes value off1; 1g; 47 The strengths of the uctuators has meanhbi = 0:2 GHz and standard deviation b=hbi = 0:2. For the ground state population (where the ground state is dened w.r.t the system Hamiltonian without the stochastic noise), its evolution with an annealing time t f = 1s is plotted in Fig. 3.9(a). 0 2 4 6 8 10 time (ns) 0.4 0.5 0.6 0.7 0.8 0.9 1 Instantaneous GS Population (a) 0 2 4 6 8 10 time (ns) 0 0.2 0.4 0.6 0.8 1 (b) Figure 3.9: Pure dephasing noise: min = 0:1 GHz, max = 10 GHz. Number of uctuators is 10. hbi = 0:2 GHz, b=hbi = 0:2. Results are averaged over 30k trajectories. The initial value of each uctuator is randomly sampled according to the thermal equilibrium distribution p eq = 0:08 (i.e. each uctuator is 0:08 more likely to be initialized in state 1 than state1). (a) Instantaneous Ground state population, (b) (Absolute value) of o- diagonal elements in the computational basis. For the coherence, because of the time-dependent system Hamiltonian, the o-diagonal part of the density matrix (in the computational basis) now has a nonzero imaginary com- ponent. Therefore, the absolute valuej 01 (t)j = 1 2 jhM + (t)ij is plotted. Its evolution with an annealing time t f = 1s is plotted in Fig. 3.9(b). 48 3.4 Eect of temperature,T 1 ,T 2 decay time, and strength of uctuators We investigate the temperature dependence of 1=f telegraph noise simulation on an annealing Hamiltonian; and look into the T 1 , T 2 decay due to the 1=f telegraph noise. 3.4.1 Two temperature limits The initial distribution of uctuators depends on temperature. The initial value of each uctuator is randomly sampled according to the thermal equilibrium distributionp eq , which is related to the temperature T by [50]: p eq = 1 (e =2k B T )'=2k B T: (3.20) We now perform simulations in the two temperature limits (T!1 and T! 0). Here we study the system evolution with the following annealing Hamiltonianm whose initial state is also its ground state: ^ H sys = A 2 1 t t f x + t t f z + 1 2 X i b i i (t) z : (3.21) where: A = 1 GHz; The sum over uctuator-i with the uctuator i (t) distributed in 1= i with i 2 [ min ; max ]; min = 0:01 GHz, max = 1 GHz. i (t) switches with i =2 and takes value off1; 1g; The strengths of the uctuators has meanhbi = 0:2 GHz and standard deviation b=hbi = 0:2. 49 3.4.1.1 Innite temperature T!1 Take T ! 1, the thermal equilibrium distribution of p eq ! 0. In this case, the ini- tial distribution of uctuators has equal probability of being in spin up or spin down, i.e. P ( i (0) = 1) = P ( i (0) =1) = 1 2 . Fig. 3.10 shows that the steady state is a maxi- mally mixed state, with the (instantaneous) ground state approaching 1=2, and o-diagonal elements in computational basis decaying to 0. 0 20 40 60 80 100 time (ns) 0.4 0.5 0.6 0.7 0.8 0.9 1 Instantaneous GS Population (a) 0 20 40 60 80 100 time (ns) 0 0.2 0.4 0.6 0.8 1 (b) Figure 3.10: High temperature limit, pure dephasing noise: min = 0:01 GHz, max = 1 GHz. Number of uctuators is 20. hbi = 0:2 GHz, b=hbi = 0:2. Results are averaged over 30k trajectories. The initial value of each uctuator is randomly sampled according to thermal equilibrium distribution p eq = 0 (i.e. each uctuator is equally likely to be initialized in state 1 than state1), since T !1. (a) Instantaneous Ground state population, (b) (Absolute value) of o-diagonal elements in the computational basis. 3.4.1.2 Zero temperature T! 0 For T ! 0, the thermal equilibrium distribution of p eq ! 1. In this case, the initial distribution of uctuators are all collectively spin up, i.e. P ( i (0) = 1) = 1. Fig. 3.11 shows that the steady state is also a maximally mixed state, with (instantaneous) ground state approaching 1=2, and o-diagonal elements in computational basis decaying decaying to 0. Assuming the strengths and frequencies of the uctuators do not depend on temperature, 50 with temperature T! 0, we can see that the uctuators destroy the coherence faster and lead to oscillations in the (instantaneous) GS population early in the anneal. 0 20 40 60 80 100 time (ns) 0 0.2 0.4 0.6 0.8 1 Instantaneous GS Population (a) 0 20 40 60 80 100 time (ns) 0 0.2 0.4 0.6 0.8 1 (b) Figure 3.11: Zero temperature limit, pure dephasing noise. Same uctuator parameters as in Fig. 3.10. Results are averaged over 30k trajectories. p eq = 1, since T ! 0. (a) Instantaneous Ground state population, (b) (Absolute value) of o-diagonal elements in the computational basis. 3.4.2 T 1 and T 2 time 3.4.2.1 T 1 time To measure the T 1 time, we initialize the state in thej1i state and evolve it with the following Stochastic Hamiltonian with uctuators in the transverse direction x (Note that with uctuator as z , the T 1 time is essentially innity since it commutes with the system Hamiltonian): ^ H sys = E 2 z + 1 2 X i b i i (t) x ; (3.22) where i (t) switches with i =2. The ipping frequencies i (t) of the uctuator-i is randomly drawn from the distribution of 1= i At large timetT 1 , after averaging many trajectories, 51 the resulting density matrix under this initial condition and Hamiltonian can be expressed as: (t) = 0 B @ 1 2 e t T 1 + 1 2 0 0 1 2 e t T 1 + 1 2 1 C A (3.23) The T 1 time can then be found by nding the decay prole of the diagonal elements of Eq. (3.23). For the example of a qubit with energy of E=2 = 1 GHz and the uctuator strength beinghbi = 0:046 GHz, the T 1 time and decay prole is plotted in Fig. 3.12. 0 0.2 0.4 0.6 0.8 1 0.4 0.5 0.6 0.7 0.8 0.9 1 1.1 Figure 3.12: Transverse noise: min = 0:01 GHz, max = 1 GHz. Number of uctuators is 200. hbi = 0:046 GHz, b=hbi = 0:2. Results are averaged over 1k trajectories. The initial value of each uctuator is randomly sampled according to thermal equilibrium distribution p eq = 0:08 (i.e. each uctuator is 0:08 more likely to be initialized in state 1 than state 1). Error bars are 2 over 1k trajectories. 52 3.4.2.2 T 2 time To measure the T 2 time (without presence of Longitudinal noise), we initialize the state in thej+i state and evolve it with the following Stochastic Hamiltonian with uctuators in the longitudinal direction z : ^ H sys = E 2 z + 1 2 X i b i i (t) z ; (3.24) where i (t) switches with i =2. After averaging many trajectories, the resulting density matrix under this initial condition and Hamiltonian can be expressed as: (t) = 0 B @ 1 2 1 2 e i!t e t T 2 1 2 e i!t e t T 2 1 2 1 C A with (0) = 0 B @ 1 2 1 2 1 2 1 2 1 C A : (3.25) The T 2 time can then be found by nding the decay prole of expectation value of the j+i state. Its time evolution according to Eq. (3.25) is: h+j(t)j+i = 1 2 + 1 2 e t T 2 cos!t: (3.26) As t T 2 , the qubit goes to the maximally mixed state. For the example of a qubit with energy of E=2 = 1 GHz and uctuator strength ofhbi = 0:046 GHz, the expectation value and decay prole is plotted in Fig. 3.13. 3.4.3 Dierent uctuator strengths We proceed to study the dependence on qubit evolution with the average uctuator strength hbi. Consider the following system Hamiltonian and uctuator parameters. ^ H sys = E 2 z + 1 2 X i b i i (t) x ; (3.27) 53 0 0.02 0.04 0.06 0.08 0.1 0 0.2 0.4 0.6 0.8 1 1.2 Figure 3.13: Longitudinal noise: min = 0:01 GHz, max = 1 GHz. Number of uctuators is 100. hbi = 0:046 GHz, b=hbi = 0:2. Results are averaged over 1k trajectories. The initial value of each uctuator is randomly sampled according to the thermal equilibrium distribution p eq = 0:08 (i.e. each uctuator is 0:08 more likely to be initialized in state 1 than state1). Here we denote the decay prole as e t=T 2 since there is no longitudinal noise (i.e. T 1 =1). Error bars are 2 over 1k trajectories. where: E = 1 GHz; min = 0:01 GHz, max = 1 GHz. i (t) switches with i =2 and takes value off1; 1g; We compare the two cases wherehbi = 2 GHz andhbi = 0:2 GHz. Fig. 3.14, we see that a higher coupling strength of the uctuators speeds up the evolution to a maximally mixed state. 54 0 1 2 3 4 5 time (ns) 10 -7 0.4 0.5 0.6 0.7 0.8 0.9 1 Instantaneous GS Population Figure 3.14: Longitudinal noise: min = 0:01 GHz, max = 1 GHz. Number of uctuators is 20. Results are averaged over 3k trajectories. The initial value of each uctuator is randomly sampled according to the thermal equilibrium distribution p eq = 0:08 (i.e., each uctuator is 0:08 more likely to be initialized in state 1 than state1). Blue line: hbi = 0:2 GHz; Green line:hbi = 2 GHz. b=hbi = 0:2. 3.5 Modeling 1=f noise for superconducting qubit: Capacitively shunted ux qubit (CSFQ) In this and the next section, we study how to model 1=f noise for a superconducting circuit Hamiltonian, and how the open-system dynamics caused by 1=f noise is dierent from that of the weak-coupling limit. We study two kinds of superconducting circuit Hamiltonians, namely Capacitively shunted ux qubit (CSFQ) and Transmon. 55 3.5.1 Capacitively shunted ux qubit (CSFQ) circuit Hamiltonian The capacitively shunted ux qubit (CSFQ) circuit HamiltonianH S (s) takes the form of [59]: H S (s) = E C 2 ^ n 2 2E J cos ( ^ ') + 2E J s 1 +d 2 tan 2 x (s) 2 cos x (s) 2 cos (2 ^ ' +' z (s)' 0 (s)) ; which can be rewritten as displacement operators in the charge basis: H S (s) = E J E C 2E J ^ n 2 ( ^ d 1 + ^ d y 1 )+ cos x (s) 2 r 1 +d 2 tan 2 x (s) 2 e i(z (s) 0 (s)) ^ d 2 +e i(z (s) 0 (s)) ^ d y 2 ! ; where (as in Fig. 3.15): E J =I z 0 2 is the Josephson junction energy; E C is the mean Capacitive energy = e 2 C sh ; ^ n is the charge operator (with the same dimension as H S ); ^ d 1 is the charge displacement operator by 1 cooper pair (with the same matrix dimension as H S ); < 1 is ratio of the current in one of the small x-loop junctions to the current in the large z-loop junction. d = E J1 E J2 E J1 +E J2 is the x-loop junction asymmetry; x (s) is the time-dependent x (barrier) bias phase; z (s) is the time-dependent z (tilt) bias phase; 0 = tan 1 d tan x(s) 2 ; ^ d 2 is the charge displacement operator by 2 cooper pairs (with the same dimension as H S ). 3.5.2 CSFQ: Persistent current operator as uctuator One important question is how to include the uctuators and form the overall stochas- tic Hamiltonian for 1=f noise simulations. The total Stochastic Hamiltonian H(s) can be 56 Figure 3.15: Schematic of the experimental set up. Capacitively shunted ux qubit (CSFQ) in the center, with curved arrows showing x and z uxes threading their corresponding loops. Source: Ref. [59]. formulated as uctuations in the persistent current operator, i.e., the uctuators are all proportional to the persistent current operator. H(s) =H S ( x (s); z (s)) + 1 2 X i b i i (s)I p (s); (3.28) where H S (s) is the time-dependent circuit Hamiltonian; I p (s) = I p (( x (s); z (s))) is the persistent current operator. Qubit CSFQ uctuation z I p (s) Table 3.1: Form of uctuation (Qubit vs CSFQ) Table. 3.1 compares the form of uctuators for qubit and superconducting circuits. 57 3.5.3 Case studies We study the time evolution under the circuit Hamiltonian, with the parameters (4 decimal gures) and ux schedules given in Eq. (3.29) (which also appears in the s-curve studies in [59]). E J = 592:9434=2 GHz; E c = 5:4092GHz; = 0:46; d = 0:1; x (s) = 0:96s + 1:04; z (s) = (0:326s + z (0)): (3.29) 0 0.2 0.4 0.6 0.8 1 -180 -160 -140 -120 -100 GHz 0th 1st 2nd 3rd 4th Figure 3.16: Spectrum with z (0)= =1:21. The initial tilt bias z (0) is an external parameter whose value we can set. Here we focus on z (0)= =1:21 (the spectrum of the rst 5 levels is plotted in Fig. 3.16). We consider the cases of closed-system evolution, open-system evolution from the adiabatic master equation, and open-system evolution from the uctuator 1=f simulation. For closed system evolution and the adiabatic master equation, we truncate the energy levels of the system to 10, and plot the populations of the rst 5 levels. For uctuator 1=f simulation, again beneting from the fact that in the trajectories approach we simulate the state vector evolution rather than the density matrix's , we can eciently simulate the populations of all 121 levels of the CSFQ circuit Hamiltonian. 58 3.5.4 Closed-system annealing evolution We simulate the closed-system annealing evolution for t f = 6 ns with value of z (0)= = 1:21, and plot the results in Fig. 3.17(a). In the later stage of the anneal, there is a sequence of diabatic transitions of population to excited states due to the small gaps at that point. The population rst transfers from the ground state to the rst excited state, then to the 2-nd excited state, then to the 3-rd excited state, and ultimately ends up in the 4-th excited state. 3.5.5 Adiabatic master equation simulation We perform the adiabatic master equation simulations for t f = 6 ns with z (0)= =1:21. We set the temperature to be 20mK, g 2 to be 1e 4 and the cuto frequency ! c to be 4GHz. Since the Lindblad operator is dened based on energy dierences, we also additionally set the absolute tolerance of energy dierence (!) to be 0:1 2 Hz for the grouping of Lindblad operators. For example, ifj! ab ! cd j < 0:1 2 Hz, then L ab and L cd are combined into one Lindblad operator L(!). The results are plotted in Fig. 3.17(b). We observe the same sequence of diabatic transitions to high excited states, but they are relatively suppressed. At the end of the anneal there is still have more than half of the population in the ground state. 3.5.6 1=f noise simulation We perform the 1=f noise simulation with the stochastic Hamiltonian Eq. (3.28). We use 10 uctuators in the form ofI p (s). The ipping frequency is sampled from 1= distribution with 2 [ min = 0:01; max = 1] GHz,hbi = 0:46 GHz and b=hbi = 0:2. z (0)= =1:21. The results are shown in Fig. 3.17(c). Compared to closed system evolution and AME simulations, not only are diabatic transitions further suppressed, but there are also population transfers during the annealing. The 1=f noise is also more detrimental to the qubit's lower-level 59 0 1 2 3 4 5 6 0 0.2 0.4 0.6 0.8 1 Eigenstate population 0thpop 1stpop 2ndpop 3rdpop 4thpop (a) 0 1 2 3 4 5 6 0 0.2 0.4 0.6 0.8 1 Eigenstate population 0thpop 1stpop 2ndpop 3rdpop 4thpop (b) 0 1 2 3 4 5 6 0 0.2 0.4 0.6 0.8 1 Eigenstate population 0thpop 1stpop 2ndpop 3rdpop 4thpop (c) Figure 3.17: The eigenstate populations for z (0)= =1:21. (a) Closed-system evolution, (b) Evolution under the adiabatic master equation. (c) 1=f noise: min = 0:01 GHz, max = 1 GHz. Number of uctuators is 100. hbi = 0:46 GHz, b=hbi = 0:2. Results are averaged over 1k trajectories. The initial value of each uctuator is randomly sampled according to the thermal equilibrium distribution p eq = 0:08. Some of the population is lost to higher excited states. populations. In fact, at the end of the anneal, the sum of the rst 5 levels is no longer 1, since some of the population is lost to even higher excited states. 60 3.6 Modeling 1=f noise for superconducting qubit: Transmon 3.6.1 Transmon Circuit Hamiltonian The transmon circuit Hamiltonian H S (s) takes the form of: H S = 4E C ^ n 2 E J + s 1 +d 2 tan 2 x 2 cos x 2 cos ( ^ ' 0 ) ; which can be rewritten as: H S = 4E C ^ n 2 E J + s 1 +d 2 tan 2 x 2 cos x 2 1 2 e i 0 ^ d 1 +e i 0 ^ d y 1 ; where: E J + is the total Josephson junction energy = (E J1 +E J2 ); E J is the Josephson junction energy = (E J1 +E J2 ) 2 ; E C is the mean Capacitive energy; ^ n is the charge operator (with the same dimension as H S ); ^ d 1 is the charge displacement operator by 1 cooper pair (with the same dimension as H S ); d = E J1 E J2 E J1 +E J2 is the x-loop junction asymmetry; x is the time-dependent x (barrier) bias phase; 0 = tan 1 d tan x 2 ; 3.6.2 Transmon: Form of uctuators The transmon circuit Hamiltonian H S (s) is time-independent. Particularly, the rst two energy eigenstatesj0i andj1i do not evolve with time. Therefore, in addition to using the persistent current operator as the uctuator, one can also dene the uctuator to be a logical z =j0ih0jj1ih1j operator as in the longitudinal noise for qubit (Sec. 3.3.1) Denote the form of uctuator by A. The total Stochastic Hamiltonian H(s) can then be formulated 61 as uctuation of the longitudinal noise or persistent current operator (See table. 3.2). In general, H(t) =H S ( x ) + 1 2 X i b i i (t)A; (3.30) where i (t) switches with i =2, with i sampled from 1= distribution in the range of f min ; max g. H S is the time-independent, transmon circuit Hamiltonian. Longitudinal noise Persistent current operator uctuation (A) z I p Table 3.2: Form of uctuation. z =j0ih0jj1ih1j. ; I p =I p ( x ) is the persistent current operator. 3.6.3 Case studies We study the time evolution under the circuit Hamiltonian, with the following parameters (4 decimal gures) and ux schedules: E J = 15:3 GHz; E c = 0:25 GHz; d = 0:03; x = 0:5: (3.31) 0 0.2 0.4 0.6 0.8 1 -10 -5 0 5 10 GHz 0th 1st 2nd 3rd 4th Figure 3.18: Spectrum with x = 0:5. 62 Here we focus on x = =0:5 (the spectrum of the rst 5 levels is plotted in Fig. 3.18). We look into the cases of open-system evolution from adiabatic master equation, and open- system evolution from the uctuator 1=f simulation. 3.6.4 Adiabatic master equation simulation We perform the adiabatic master equation simulations for t f = 10 ns with value of x = = 0:5. Set the temperature to be 20mK and ! c cuto frequency to be 8GHz. Since the Lindblad operator is dened based on energy dierence, We also additionally set the absolute tolerance of energy dierence (!) to be 0:1 2 Hz for the grouping of Lindblad operator. For example, ifj! ab ! cd j < 0:1 2 Hz, then L ab and L cd are combined into one Lindblad operator L(!). We truncate the energy levels of the system to 10. The results of coupling g = 0:008(2) and g = 0:08(2) are plotted in Fig. 3.19(a) and Fig. 3.19(b). As expected, the stronger the coupling to the environment, the faster the coherence is lost. 0 2 4 6 8 10 0 0.2 0.4 0.6 0.8 1 (a) 0 2 4 6 8 10 0 0.2 0.4 0.6 0.8 1 (b) Figure 3.19: The population inj+i for coupling (a) g = 0:008(2). (b) g = 0:08(2). Evolution under (adiabatic) master equation. 63 3.6.5 1=f noise simulation For uctuator 1=f simulation, beneted from the fact that in trajectories approach we sim- ulate the state vector evolution rather than the density matrix's , we can eciently simulate the populations of all the 21 levels, dened by the Transmon circuit Hamiltonian. We per- form the 1=f noise simulation with the stochastic Hamiltonian Eq. (3.30). We consider both form of uctuators: A = z and A = I p ( x ). x = =0:5. We use 20 uctuators in both cases. The ipping frquency is sampled from 1= distribution with 2 [ min ; max ]. The result of A = z is shown in Fig. 3.20 and that of A =I p ( x ) is shown in Fig. 3.21. In both cases we observe non-exponential decay in the population inj+i (o-diagonal elements of density matrix, after averaging over many trajectories) 0 2 4 6 8 10 time (ns) 0.2 0.4 0.6 0.8 1 Figure 3.20: A = z . The population ofj+i for x = =0:5. min = 1 GHz, max = 10 GHz. Number of uctuators is 20. hbi = 0:4600 GHz, b=hbi = 0:2. Results are averaged over 300 trajectories. p eq = 0:08. 64 0 10 20 30 40 50 60 time (ns) 0 0.2 0.4 0.6 0.8 1 1.2 Figure 3.21: A = I p . The population ofj+i for x = =0:5. min = 0:1 GHz, max = 10 GHz. Number of uctuators is 20. hbi = 0:0092 GHz, b=hbi = 0:2. Results are averaged over 300 trajectories. p eq = 0:08. 3.7 Fitting experimental data and reverse engineering of the spectrum The goal of this section is to try to evaluate the eectiveness of 1=f noise simulation in the ttings of Transmon qubits (IBM Quantum) experimental data and how dynamical decoupling can eliminate part of the 1=f noise. There are two directions. Assume that one is given either experimental data of free evolution and dynamical decoupling from a superconducting qubit experiment. One direction is reconstructing the 1=! spectrum by tting the experimental data of T 1 and T 2 times of free evolution and dynamical decoupling. This is what [56] did: nd out the parameters that t with the experimental data from free evolutions and dynamical decoupling, and reconstruct the 1=! spectrum. Recall the form ofS(!) in Eq. (3.8). By nding the correct parameters likeb, M , and m to best t with the experimental data, one can nd out the spectrumS(!). This is what we 65 do in this section. We do not know what the power spectrum of the IBM quantum computer is. Another direction is, given the power spectrum, simulate the qubit free evolutions and evolutions under dynamical decoupling. In this case, we know exactly what S(!) is in a given range of !. With this information, one can deduce the parameters like b, M , and m , and simulate the free evolutions and evolutions under dynamical decoupling with such parameters, etc. Fig. 3.22 illustrates both directions. Reverse Forward Reverse Forward 0 2000 4000 6000 8000 0.4 0.5 0.6 0.7 0.8 0.9 1 (fitting) Experimental data T1, T2 (reverse engineering) Spectrum S(ω) = A ω α γ min , γ max , b Figure 3.22: The processes of determining the spectrum [56] and performing simulations based on the spectrum. 3.7.1 Fitting results with IBM machine We try to t the IBM experimental data with our 1=f noise simulation. The goal is to reduce the sum of L2 norm of the dierence of the simulation data and the experimental data, i.e. P i ks i e i k 2 , and extract the best uctuator parameters. (b; min ; max , etc) We perform the simulations in a rotated frame, where the frame of reference is the original system Hamiltonian. 66 3.7.1.1 Free induction decay: Fluctuation A =I p We try to t the simulations with the IBM experimental data of free induction decay, where the gates along evolution are the identity gates. The uctuation chosen in our simulations is A = I p . We x the standard deviation of uctuator strength to be b=hbi = 0:2. The initial value of each uctuator is randomly sampled according to the thermal equilibrium distribution p eq = 0 (i.e., each uctuator is equally likely to be initialized in the state 1 and state1). The number of uctuators is 20. We rst output the best-t parameters (b; min ; max , etc) to the free induction decay curve with duration of t f = 5405ns, with the initial state beingj+i. The best t parameters of uctuators are b = 0:12(9:2) 10 3 ; min = 0:05 GHz; max = 5 GHz (Fig. 3.23(a)). With these best-t parameters, we can perform simulations with a new initial state. The simulation data of initial statej1i also compares well with the free induction decay experimental data (Fig. 3.23(b)). However, the simulation data of initial state j0i does not compare well experiments (Fig. 3.23(c)). If the initial state of the transmon qubit is the ground state, the experimental data suggests that it stays in the ground state, while the simulation data from 1=f noise suggests it goes to the maximally mixed state. 3.7.1.2 Dynamical decouplings: Fluctuation A = [I p ] 22 We try to t the simulations with the IBM experimental data of dynamical decoupling. We x the standard deviation of uctuator strength to be b=hbi = 0:2. The initial value of each uctuator is randomly sampled according to thermal equilibrium distributionp eq = 0. The number of uctuators is 20. Based on the tting with free evolution, we perform simulations with the XYXY dynamical decoupling sequence. We model the dynamical decoupling as ideal pulses in terms of Pauli two-level x and y matrices. Therefore, our system Hamiltonian and uctuators are also truncated to the rst two energy levels. A = [I p ] 22 . 67 0 1000 2000 3000 4000 5000 6000 0.4 0.5 0.6 0.7 0.8 0.9 1 (a) 0 1000 2000 3000 4000 5000 6000 0.8 0.85 0.9 0.95 1 (b) 0 1000 2000 3000 4000 5000 6000 0.92 0.93 0.94 0.95 0.96 0.97 0.98 0.99 1 (c) Figure 3.23: t f = 5405ns, b = 0:12(9:2) 10 3 ; min = 0:05 GHz; max = 5 GHz, uctuations in I p , 20 uctuators, 200 trajectories. Free induction decay. Rotated frame. We rst nd out the best parameters in (b; min ; max , etc) to the free induction decay curve with duration of t f = 15000ns, with initial statej+i. The best t parameters of uctuators are b = 0:12(9:2) 10 3 GHz; min = 0:05 GHz; max = 5 GHz (Fig. 3.24(a)). With these best-t parameters, we can perform simulations of dynamical decoupling (XYXY sequence between each data point) under 1=f noise, also with initial statej+i. However, as we can see in Fig. 3.24(b), the distance between the experimental data and simulation data diverges at long times. Fast noises, which cannot be removed by dynamical decoupling, are 68 present in Transmon qubit, and thus the experimental dynamical decoupling curve decays faster than pure 1=f noise simulation. Nevertheless, we increase the average uctuator strength by a bit, and can obtain a better t (Fig. 3.24(c)). 0 5000 10000 15000 0.4 0.5 0.6 0.7 0.8 0.9 1 (a) 0 5000 10000 15000 0.65 0.7 0.75 0.8 0.85 0.9 0.95 1 (b) 0 5000 10000 15000 0.65 0.7 0.75 0.8 0.85 0.9 0.95 1 (c) Figure 3.24: Fluctuations in [I p ] 22 , 20 uctuators, 300 trajectories. Rotated frame. min = 0:01 GHz; max = 0:1 GHz. (a), (b): b = 9:292910 4 GHz, (c) Modied averaged uctuator strength: b = 1:05 10 3 GHz. 69 3.8 Conclusion of this chapter In this chapter, based on the previously studied Bloch vector-magnetic eld approach, we have formulated a stochastic uctuator Hamiltonian simulation method that includes a series of uctuator terms with dierent ipping frequencies and Gaussian couplings that altogether constitute the 1=f power spectrum. We, for the rst time, show how to incorporate them into the inherently time-dependent annealing Hamiltonians, and construct an overall time- dependent stochastic Hamiltonian to perform parallelizable simulations. We thoroughly analyzed how noise with 1=f power spectrum aects the quantum annealing, and we also studied temperature eects, loss of coherence (T 1 and T 2 ) in QA under dierent uctuator parameters. We developed a exible approach, allowing the noise uctuators to take any form of multi-level operators and any noise-axis direction depending on the superconducting qubits and circuit decoherence process. We showed how to append these operators to two dierent superconducting circuit Hamiltonians. We t our simulation results with experimental data, and showed how dynamical decoupling can extend the coherence of qubits under 1=f noise. 3.9 Acknowledgments of this chapter Sec. 3.7.1 would be impossible without Vinay Tripathi's experimental data and discussions. The codes in this chapter are available at: https://github.com/USCqserver/1fnoise 70 Chapter 4 Quantum feedback control As mentioned in the conclusion 2.7 of Chapter. 2, an extension of the trajectories method is the study of weak measurement and feedback in annealing system. We devise a quantum feedback error correction method to reverse the eect (error) of thermal excitations during quantum annealing. 4.1 Introduction In QA, the solution of the computational problem is encoded as the ground state of the nal Hamiltonian. In an ideal situation, a quantum state is initiated as the ground state of an initial Hamiltonian, and is evolved to the ground state of the nal Hamiltonian. It is desirable is to maintain the state during the anneal as the ground state of the time- dependent Hamiltonian. However, due to diabatic or thermal transitions, this is not always successful, i.e., the quantum state at the end of the process is not the ground state of the nal Hamiltonian. Quantum error correction (QEC) aims to correct such quantum errors. There were previous extensive studies on quantum error correction in QA [60{62], mostly centered on using error correcting codes to perform majority voting at the end of the anneal. This chapter is about a particular active quantum error correction approach that can used during an anneal: feedback-based error correction for QA. 71 Quantum feedback-based error correction aims to manipulate the system's quantum state or trajectory to evolve the system towards some desired outcome. Depending on the delay be- tween detection and feedback, the eect of the feedback can be Markovian or non-Markovian. We study the following topics about feedback-based error correction on QA. 4.1.1 Feedback and control In active quantum error correction, one may want to know when and what types of corrections to apply to the system qubit(s) during the quantum operation/process. One way to achieve this goal is to perform measurement on the quantum system at a current time-step and use the measurement result as feedback to determine/improve the error correction for future time-steps. 4.1.2 Weak measurement One way to obtain information about quantum system is by weak measurement. Weak mea- surement is a type of quantum measurement that results in an observer obtaining very little information about the system on average, but also disturbs the state very little [63]. Weak measurement is also often referred to as continuous measurement, where (weak) measure- ment is performed continuously in successive tiny time-steps. Weak measurement arises from the system qubit successively interacting with environmental (probe) qubits which are subse- quently measured (See Fig. 4.1). Depending on the basis in which the environmental qubits is measured, one can obtain dierent types (jumpy/diusive) of measurement trajectories for the system qubit being measured. When little information is gathered during the measurement, the quantum state is almost entirely preserved. Nevertheless, taking a suciently large sample of data still permits the same average information to be extracted as a textbook projective measurement that fully collapses the state. By continuously monitoring partial information about a quantum system, 72 Figure 4.1: By successively interacting the system with environmental (probe) qubits which are subsequently measured, we can obtain the weak measurement trajectories. one can control experimental parameters in real time for a variety of purposes, such as error correction, coherence stabilization, state purication. 4.2 Thermal excitation as error: detection In QA, thermal excitation is considered an error since we want to preserve the ground state of the system Hamiltonian. To do that in feedback-based error correction, the system evolution is continuously monitored and feedback is applied conditioned on each weak measurement result. A photo-detection current I(t) of thermal excitation is obtained, 0<t<t f . With the assumption of weak-coupling limit, one can derive a master equation describing the evolution of system qubit as: d dt (t) =i [H S (t) +H LS (t);(t)] +D[A ! (t)](t) +D[A 0 (t)](t) +D[A ! (t)](t): (4.1) Specically, A ! (t)/j" 1 (t)ih" 0 (t)j (4.2) 73 is the Lindblad operator of thermal excitation in the instantaneous energy eigenbasis. D[ ] refers to the decoherence channel (excitation/dephasing/damping). This master equation can be unraveled into a stochastic Schr odinger equation describing the weak measurement. The system evolution under annealing is continuously monitored and feedback is applied conditioned on each continuous measurement result. The continuous measurement results are in the form of photodetection current I(t), 0tt f . The photodetection current I(t) corresponds to the detection of excitation of the system qubit. I(t) = 1 when excitation is detected, = 0 otherwise (See Fig. 4.2). Conditioned on this current, we can obtain an evolution (unravelling in the quantum jump picture) of the conditioned state I (t). 0 1 Figure 4.2: Photodetection current I(t). 4.3 Delay between detection and feedback After the detection of a thermal excitation (I(t) = 1) from ground to excited state, we then want to apply feedback to correct such errors. However, there usually is a delay between detection and feedback. The most general expression for the eect of the feedback on the conditioned state I (t) is [64]: [ _ I (t)] fb =F[t; I [0;t) ] I (t); (4.3) where I [0;t) is the complete photocurrent record from time 0 tot, andF[t; I [0;t) ] is a superop- erator. Assuming that the feedback is of linear functional form, the feedback evolution can be expressed as: [ _ I (t)] fb = Z 1 0 h(s)I(ts)L f I (t)ds; (4.4) 74 whereL f is the Liouville superoperator representing the feedback control. If the delay is xed as , an equation incorporating the eect of the feedback is: [ _ I (t)] fb =I(t)L f I (t): (4.5) Note that Eq. (4.5) is non-Markovian for nite. An analytical solution of Eq. (4.5) is hard to obtain and numerical simulation is needed. 4.4 Feedback superoperatorL f : Lindbladian The feedback superoperatorL f can take the form of a Lindbladian or a Hamiltonian. In the form of a Lindbladian: L f [ ] =F (t) F y (t) 1 2 F y (t)F (t); : (4.6) Physically, for the purpose of QA, this Lindbladian represents a cooling/damping channel in which the qubit is relaxed from an excited state to the ground state. F (t) can take many forms. For example, F (t) = A y ! (t) =j" 0 (t)ih" 1 (t)j, the complex conjugate of the thermal excitation operator [Eq. (4.2)] a time ago. This exactly cancels the error caused by thermal excitation if = 0. With this form of F (t), we can devise a quantum feedback algorithm { see Algorithm. 1. We perform such trajectories feedback simulation in Section 4.5. 75 Algorithm 1 Quantum feedback algorithm with delay procedure Feedback with delay(t f ;t j = 0;) . Initialize t j = 0. while t j +t f do Evolve according to quantum trajectory method. if Excitation A ! (t) occurs at t =t j toj"(t j )i (Indicated by I(t) = 1) then Apply feedback F (t j +) =A y ! (t j ) at t j + toj"(t j +)i. end if end while end procedure 4.5 Case studies We consider the case of a single qubit with annealing Hamiltonian, i.e., H S (t) = A(t) x + B(t) z . TheA(t) =A(1t=t f ) andB(t) =B(t=t f ) schedules are the linear schedules, with A = 1:57 GHz and B = 2:01 GHz (in units ~ = 1). We use F (t) =A y ! (t) =j" 0 (t)ih" 1 (t)j, with dierent delays 0t f . In Fig. 4.3, we simulate the performance of the feedback error protocol for dierent delay and plot the ground state populations along the anneal. For the purpose of QA we want to maximize the nal ground state population, i.e. ground state population at t f . This serves as a metric to measure the eectiveness of the quantum feedback error correction. Note that the nal ground state population need not decrease as delay increases. There is an optimal we can search for. In Fig. 4.4, we can see that there is a nonzero delay where nal ground state population reaches a local maximum. Thus there is an experimentally allowed delay time where the nal ground state population is maximized. 76 0 20 40 60 80 100 0.55 0.6 0.65 0.7 0.75 0.8 0.85 0.9 0.95 1 Instantaneous GS Population no feedback instant 5us 10us 30us Figure 4.3: Evolution of instantaneous ground state population for a total time oft f = 100s and temperature 2:6 GHz, with feedback applied with some dierent delays (0; 5; 10; 30s). Results are obtained from averaging 5k quantum trajectories of the quantum feedback algo- rithm. Error bars are 2 over 5k trajectories. 4.6 Feedback superoperatorL f : Hamiltonian and basis of feedback L f can also take the form of Hamiltonian [64]. L f [ ] =i [H fb (t); ] ; (4.7) where H fb (t) is called the feedback Hamiltonian. Some examples of H fb (t) are • x in the energy eigenbasis H fb (t) = 2 (j" 0 (t)ih" 1 (t)j +j" 1 (t)ih" 0 (t)j) ; (4.8) 77 0 10 20 30 40 50 60 70 80 90 100 0.55 0.6 0.65 0.7 0.75 0.8 0.85 0.9 0.95 1 Final GS population Figure 4.4: Linear Schedule. Final ground state population. Single qubit. Total annealing time t f = 100s. Dashed purple line is the nal ground state population without feed- back. Average over 5k trajectories with Lindbladian feedback. There is a non-small (thus experimentally allowed) 60us where the Final GS population is maximized. • x in the computational basis H fb (t) =H fb = 2 x ; (4.9) • z in the computational basis H fb (t) =H fb = 2 z : (4.10) These are Pauli operators in dierent bases. Note that in the decoherence model of the weak coupling limit, any excitation/relaxation takes place in the energy eigen-basis [13], thus one would expect that the feedback Hamiltonian x in the energy eigenbasis performs better 78 than x in the computational eigenbasis (see simulations in Fig. 4.5). Note also that the signicance of the 2 implies that the feedback Hamiltonian acts instantly (pulse/gate-like). Another advantage of the feedback Hamiltonian formalism is that for delay ! 0, one can obtain the feedback master equation that explains the ensemble average of the eect of feedback. The feedback master equation is d dt (t) =i [H S (t) +H LS (t);(t)] + e iH fb (t) A ! (t)(t)A y ! (t)e iH fb (t) 1 2 n A y ! (t)A ! (t);(t) o +D[A 0 (t)](t) +D[A ! (t)](t) (4.11) This also allows direct master equation simulation, which can be faster than trajectories simulations for a small number of qubits. Fig. 4.5 below compares the two approaches of master equation simulation and the quantum feedback trajectories algorithm, over the three feedback bases mentioned in Eqs. (4.8), (4.9) and (4.10). 79 0 2 4 6 8 10 0.6 0.7 0.8 0.9 1 Algo ME Figure 4.5: Evolution of ground state population under feedback error correction of dierent feedback Hamiltonians ( = 0). The results of quantum jump trajectories (averaged over 5k trajectories) of quantum feedback algorithm agree with the feedback master equation solutions. Error bars are 2 over 5k trajectories. 80 Chapter 5 Reverse Annealing: Iterative reverse annealing We want to apply our open-system simulation techniques the study of open-system eects on a particular QA application called reverse annealing (RA). Reverse annealing has been predicted to have better performance than forward annealing in several systems [16, 19{25]. Reverse annealing can be generally separated into two formalisms: iterated reverse an- nealing (IRA) and adiabatic reverse annealing (ARA). ARA is similar to standard quantum annealing, but with an extra initialization Hamiltonian. IRA does not have such term but starts at s = 1 with a random Initial state that is usually an eigenstate of the problem Hamiltonian. Iterated reverse annealing (IRA) involves (1) annealing backward from a known classical state to a mid-anneal state of quantum superposition, (2) searching for optimum solutions at this mid-anneal point while in the presence of an increased transverse eld (quantum state), and then (3) proceeding forward to a new classical state at the end of the anneal. In this chapter, we study the open-system eects on iterative reverse annealing and perform simulations of IRA in the weak coupling limits to bath, for p-spin model with p = 3. This chapter is based on [65]. In Chapter. 6, we perform the IRA simulations with more realistic experimental settings of D-Wave annealers, for p-spin model with p = 2. In Chapter. 7, we study the open-system eects on adiabatic reverse annealing (ARA), for p-spin model with p = 3. 81 5.1 Introduction In iterative reverse annealing, the initial state is an eigenstate of the nal problem Hamilto- nian and the transverse eld is cycled rather than strictly decreased as in standard (forward) quantum annealing. We present a numerical study of the reverse quantum annealing protocol applied to thep-spin model (p = 3), including pausing, in an open system setting accounting for dephasing in the energy eigenbasis, which results in thermal relaxation. We consider both independent and collective dephasing and demonstrate that in both cases the open system dynamics substantially enhances the performance of reverse annealing. Namely, including dephasing overcomes the failure of purely closed system reverse annealing to converge to the ground state of the p-spin model. We demonstrate that pausing further improves the success probability. The collective dephasing model leads to somewhat better performance than independent dephasing. The protocol we consider corresponds closely to the one imple- mented in the current generation of commercial quantum annealers, and our results help to explain why recent experiments demonstrated enhanced success probabilities under reverse annealing and pausing. 5.2 Iterative reverse annealing: the protocol Standard quantum annealing aims at solving optimization problems by employing quantum uctuations that are slowly decreased to zero, to eciently explore the solution space [1, 3, 4]. A system of N qubits is prepared at t = 0 in the ground state of a transverse eld Hamiltonian, i. e., the uniform superposition over the = 2 N computational basis states fj0i;j1i;:::;j 1ig. The magnitude of transverse eld is then slowly decreased, while the magnitude of the HamiltonianH 0 , encoding the optimization problem, is simultaneously increased. The adiabatic theorem guarantees that if the evolution is long on the timescale set by the inverse of the minimal gap = min t [E 1 (t)E 0 (t)] between the ground state and the rst excited state (we set } = 1 henceforth), then at the end of the anneal t = the system 82 0 10 20 30 40 0.0 0.2 0.4 0.6 0.8 1.0 A( s) B( s) Schedule (GHz) s (a) 0.0 0.2 0.4 0.6 0.8 1.0 0 100 200 300 400 500 s inv t inv s t (ns) Standard annealing Reverse annealing (b) Figure 5.1: (a) Annealing schedules (in units such that ~ = 1) as a function of the annealing fractions(t), chosen to be similar to the schedules of the D-Wave processors. (b) Annealing fraction s(t). The blue solid curve represents standard forward quantum annealing of total annealing time = 500ns. populates the target ground state of H 0 with a probability P 0 that approaches unity [66, 67]. However, any nite sweep rate leads to diabatic Landau-Zener transitions at the avoided crossings [68], thus reducing the success probabilityP 0 of the adiabatic algorithm. Therefore, the output state of quantum annealing is in general a trial solution of the optimization problem, ideally having a large overlap with the target ground state. Reverse annealing instead aims at rening an already available trial solution [16]. For instance, NP-hard optimization tasks are solved using heuristics, whose output is often an approximation of the true global minimum. The algorithm of reverse annealing is the following. 1. At t = 0, the system is prepared in the trial solution state. 2. Quantum uctuations are increased, until a turning point t inv is reached during the dynamics. This ends the reversed part of the dynamics. 83 3. After the turning point, the dynamics follows the standard quantum annealing sched- ule; quantum uctuations are decreased until t = , when the state is eventually measured. Careful choices of the turning point, and of the initial state, can lead to an improvement in the solution. Moreover, this scheme can also be repeated multiple times, each time starting from the output of the previous stage; hence the terminology of iterated reverse annealing. Assuming that each iteration improves the success probabilityP 0 , this procedure can systematically improve the quality of the solution. We consider the following time-dependent Hamiltonian, suitable for IRA but not ARA (which includes another term): H(t) =A s(t) V TF +B s(t) H 0 ; V TF = N X i=1 x i : (5.1) Since we only consider the IRA protocol, henceforth we simply refer to it as reverse annealing. The function of time s(t) in Eq. (5.1) is the annealing fraction (or dimensionless time) and satises 0 s(t) 1 for all t. The two functions A(s) and B(s) determine the annealing schedule, which we choose to match the annealing schedule of the D-Wave processors [see Fig. 5.1(a)]. They satisfy A(0)B(0) and B(1)A(1). The functional form ofs(t) distinguishes between standard (forward) and reverse anneal- ing. In standard quantum annealing the dimensionless time is dened ass(t) =t=, being the annealing time. Thus, s(t) is a monotonic function of t, and is represented in the plane (t;s) by a straight line going from (0;s 0 ) to (;s 1 ), where s 0 =s(0) = 0 and s 1 =s() = 1. This is shown in Fig. 5.1(b) using a blue solid line, for = 500 ns. During standard quantum annealing, quantum uctuations are very large at t = 0, and decrease monotonically until t =. In contrast, in reverse annealing s 0 = s 1 = 1. Starting from s = s 0 , where quantum uctuations are zero, the annealing fraction is rst decreased until it reaches the inversion 84 point,s =s inv , at a timet =t inv . In this rst branch, quantum uctuations are increased. At s =s inv , the annealing fraction is then increased towards s =s 1 , and quantum uctuations are decreased again to zero. In Fig. 5.1(b), we show a typical function s(t) for a reverse annealing of annealing time = 500 ns with an inversion point (t inv = 200 ns;s inv = 0:6), using a red dashed line. In general,s inv andt inv can be chosen independently of each other. However, in this work we choose the following linear relation, in order to have only one free parameter: t inv =(1s inv ) ; s inv 6= 0; 1: (5.2) In this way, we have that s(t) = 8 > > < > > : 1t= for tt inv ; 1s inv s inv t + 2s inv 1 s inv for t>t inv : (5.3) Another possible choice would be to x t inv so that the two slopes are the same, i. e., t inv ==2 for all choices of s inv . This is similar in spirit to what is discussed in Ref. [23]. In what follows, we will focus on the fully-connected ferromagneticp-spin model, a model with a permutationally invariant Hamiltonian and a nontrivial phase diagram, often used as a benchmark for the performance of quantum annealing [69, 70]. 5.3 Ferromagnetic p-spin model The Hamiltonian of the ferromagnetic p-spin model, in dimensionless units, reads H 0 = N 2 1 N N X i=1 z i ! p ; (5.4) with p 2. For even p, there are two degenerate ferromagnetic ground states, whereas for odd p the ferromagnetic ground state is nondegenerate. For p = 2, the Hamiltonian of 85 Eq. (5.1) is subject to a second-order QPT in the thermodynamic limit at the critical point s c (or, equivalently, t c = s c ), separating a para- and a ferromagnetic phase. For p > 2, the QPT is rst-order [71]. The presence of QPTs aects also the nite-size behavior of the system, as the minimal spectral gap , found ats (t =s ), closes asN 1=3 forp = 2 or exponentially in N for p > 2. s (i. e., t ) approaches s c (i. e., t c ) as N!1. First-order QPTs are especially detrimental for quantum annealing, as the annealing time has to grow exponentially with the system size to compensate the closure of the gap at s = s . In the following, we will focus on the case p = 3. We can dene the total spin operator S = (S x ;S y ;S z ) and the dimensionless magnetiza- tion operators m = 1 N S ; S = N X i=1 i ; 2fx;y;zg (5.5) that allow the rewriting of the time-dependent Hamiltonian of Eq. (5.1) as H(t) = A(t) 2 S x B(t)N 2 m p z : (5.6) Since [S 2 ;S z ] = 0 and both are Hermitian, they share an orthonormal basisfjS; S ig such that the eigenvalues of S 2 are S(S + 1) with S2f0; 1=2; 1;:::;N=2g for even N and S2 f1=2; 1;:::;N=2g for oddN, and the eigenvalues ofS z are S 2fS;S + 1;:::;Sg. In the subspace with maximum spinS =N=2, we instead label the basis states asjwijN=2wi, withw2f0; 1;:::;Ng. These are the eigenstates ofm z with eigenvaluesm = 12w=N. The target state is the ferromagnetic ground statej0i, i. e., the eigenstate of m z with eigenvalue m = 1. The Hamiltonian of thep-spin model commutes with S 2 . Hence sectors diering byS do not become coupled under the dynamics generated by the p-spin model Hamiltonian. Since the ferromagnetic ground state and the initial one belong to the subspace with maximum spin S = N=2, the interesting dynamics occurs in this subspace, whose dimension scales 86 0.0 0.2 0.4 0.6 0.8 1.0 0.00 0.25 0.50 η = 0 s Δ = 0.309 τ = 100 ns Success probability s inv m 0 = 0.9 m 0 = 0.8 m 0 = 0.0 m 0 = -1.0 (a) 0.0 0.2 0.4 0.6 0.8 1.0 -1.0 -0.5 0.0 0.5 1.0 Maximum success probability m 0 τ = 100 ns τ = 1000 ns 10 -15 10 -10 10 -5 10 0 -1.0 -0.5 0.0 (b) Figure 5.2: (a) Success probability in unitary reverse annealing as a function of the inversion point s inv , for several values of the magnetization of the initial state. The dashed vertical line indicates s inv =s 0:309. The annealing time is = 100 ns. We sampled the interval s inv 2 (0; 1) using a step size of s = 0:002, and repeated the dynamics for each choice of s inv . (b) Maximum success probability achievable with unitary reverse annealing, as a function of the magnetization of the initial state, for two annealing times: = 100 ns and = 1000 ns. The inset zooms in on the region m2 [1:0; 0:0]. Note the logarithmic scale on the vertical axis. linearly with the number of qubits: = N + 1. This fact enables us to perform numerical calculations with relative large numbers of qubits N. In the following, we will start reverse annealing in each of theN excited statesfj1i;:::;jNig of H 0 in the symmetric sector with S = N=2. The similarity to the ferromagnetic ground state is quantied by the corresponding starting eigenvalue of m z , denoted m 0 . Note that the wth excited state diers from the ferromagnetic ground state by w spin ips. There- fore, the initial state and the target solution dier by a fraction c = N " =N = 1w=N of up-aligned qubits. These parameters are also related to the Hamming distance d H , via d H =NN " =N(1c). 87 5.4 Unitary dynamics In this section we study the closed system case of a system of N = 20 qubits, with p = 3. For our choice of parameters and in terms of the annealing schedules shown in Fig. 5.1(a), the p-spin system has a minimal gap 2:45 GHz at s 0:309. The annealing time is = 100 ns. In Fig. 5.2(a), we report the ground state population P 0 at t = , as a function of the inversion points inv , for several initial states: m 0 = 0:9, 0:8, 0 and1. Recall that the target ground statej0i has m = 1. We focus on the region s inv 2 (0:0; 0:5]. The rightmost part of Fig. 5.2(a) corresponds to cases in which the anneal is reversed too early, i. e., fort inv <t ands inv >s . The system does not cross its quantum critical point, and the success probability is zero. Therefore, no eects on the outcome of the procedure are visible, as the dynamics is slow compared with the minimal inverse level spacing and diabatic transitions are exponentially suppressed. Thus, the system is forced to stay in its initial state, or transition to other excited states. In fact, avoided crossings between pairs of excited eigenstates occur at s>s for this model, and Landau-Zener processes can further excite the p-spin system. On the other hand, if s inv < s the system crosses the minimal gap twice. Here, the success probability benets from Landau-Zener processes, inducing transitions towards the ground state. In this region, we also note some non-adiabatic oscillations of the success probability, due to the nite annealing time. These oscillations are more evident for large m 0 . As expected from the adiabatic theorem, they are suppressed for longer annealing times. For instance, we veried that they are no longer visible for = 1000 ns (not shown). The sharp rise of the success probability for m 0 = 0:9 occurs exactly at s inv = s . For smaller values of m 0 , the success probability rises more smoothly, as the ground state is reached after a preliminary sequence of Landau-Zener transitions between pairs of excited states, whose corresponding avoided crossings occur at s > s . For m 0 = 0:8, a very small rise of the success probability can still be observed around s inv = s . This is due to the fact that 88 during the reverse annealing, the system prepared in the second excited state rst encounters an avoided crossing with the rst excited state, where part of the population is transferred to the latter, and then the avoided crossing with the ground state, where the system populates its ground state. After reversing the dynamics, the two avoided crossings are encountered again (in the reverse order) and part of the population gets excited, thus reducing the success probability P 0 . As expected, reverse annealing is more eective when the initial state is close to the correct ground state. Moreover, as is also clear from Fig. 5.2(a), the inversion times inv must be increasingly close to 0 for decreasingm 0 , in order to obtain a nonzero success probability at t = . This means that almost the entire dynamics is spent in the reverse part of the annealing, and the system is eventually quenched towards s =s 1 for t. Even so, if the initial state is too far in energy from the correct solution, the success probability of reverse annealing is always close to zero, as evident from the curves for m 0 = 0 and m 0 =1 in Fig. 5.2(a). The maximum success probability decreases rapidly as a function of m 0 . This is clearly seen in Fig. 5.2(b), where we report the maximum attainable success probability as a function of m 0 , for annealing times = 100 ns and 1000 ns. Increasing the annealing time reduces non-adiabaticity and results in a lower success probability, compared with that at the end of a faster reverse anneal. As shown in Fig. 5.2(b), which zooms in on the region m 0 2 [1; 0], this decrease can be of several orders of magnitude for poorly chosen trial solutions. The in uence of the annealing time is less pronounced close to m 0 = 1, and more evident for intermediate and lower values of m 0 . This is consistent with the adiabatic theorem, since a longer anneal time guarantees that the system will have a higher probability of remaining close to the initial eigenstate it has the largest overlap with (not necessarily the ground state) [66]. Adopting a conventional quantum annealing procedure, the success probability in the case of = 100 ns would be P 0 = 0:96. This value is larger than any P 0 achievable using reverse 89 0.0 0.2 0.4 0.6 0.8 1.0 0.0 0.2 0.4 0.6 0.8 1.0 η = 1 × 10 -3 s Δ = 0.309 τ = 100 ns Success probability s inv m 0 = 0.9 m 0 = 0.8 m 0 = 0.0 m 0 = -1.0 (a) 0.0 0.2 0.4 0.6 0.8 1.0 0.0 0.2 0.4 0.6 0.8 1.0 η = 1 × 10 -3 s Δ = 0.309 τ = 500 ns Success probability s inv m 0 = 0.9 m 0 = 0.8 m 0 = 0.0 m 0 = -1.0 (b) Figure 5.3: Success probability in reverse annealing as a function of the inversion point s inv , for several values of the magnetization of the initial state. The p-spin system is coupled to a collective dephasing bosonic environment as in Eq. (5.10), and the coupling strength is = 1 10 3 . The dashed vertical line denotes the time s of the avoided crossing between the ground state and the rst excited state. In (a) the annealing time is = 100 ns, in (b) = 500 ns. We sampled the interval s inv 2 (0; 1) using a step s = 0:02. All other parameters are given in the main text. annealing. However, this argument cannot be used to discredit reverse annealing for two main reasons. First, it is clear that, in the analyzed case, we are very close to adiabaticity, where conventional annealing is ecient. Second, as claried in the next Sections, the role of dissipation may strongly aect this scenario. The results of this Section are in agreement with those reported in Ref. [23]. Namely, as is clear from Fig. 5.2(b), upon iteration the IRA protocol will only decrease the success probability under unitary, closed system dynamics, unless the initial state was already chosen as the solution of the optimization problem. 90 5.5 Open system dynamics subject to dephasing-induced relaxation Physical quantum processors always interact with the surrounding environment, which in- duces decoherence and thermal excitation/relaxation, which in turn impacts the performance of quantum annealing [9, 12, 72]. We assume weak coupling between the qubit system and the environment. It can then be shown that the reduced system density matrix evolves ac- cording to a quantum master equation in time-dependent Lindblad form [73], known as the adiabatic master equation: d(t) dt =i (t);H(t) +H LS (t) +D (t) : (5.7) In Eq. (5.7),H LS (t) is a Lamb shift term andD is the dissipator superoperator, which makes the dynamics non-unitary and irreversible. They are expressed in terms of Lindblad opera- tors, inducing dephasing or quantum jumps (pumps and decays) between pairs of adiabatic energy eigenstates [27]: H LS = X ab L y ab (t)L ab (t)S ab (! ab ); (5.8) D[(t)] = X ab ab L ab (t)(t)L y ab (t) 1 2 n L y ab (t)L ab (t);(t) o : (5.9) The Lindblad operators are determined by the instantaneous eigenbasis of the system Hamil- tonian H(t) and the system-bath interaction Hamiltonian H SB . In general, this coupling Hamiltonian involves local operators that break the spin symmetry of the p-spin model, as each qubit is then coupled to its own bath. We study both this independent decoherence model and the collective decoherence model, wherein all the qubits are coupled to a collective 91 bath with the same coupling energyg, in order to preserve the spin symmetry. More speci- cally, we rst consider collective dephasing, for which the system-bath coupling Hamiltonian is H col SB =gS z B; (5.10) where B is a bath operator [e.g., B = P k (a k +a y k ) for an oscillator bath with annihilation operators a k for the kth bosonic mode]. The Lindblad operators are represented in the instantaneous energy eigenbasis of H(t) as [73]: L ab (t) =hE a (t)jS z jE b (t)ijE a (t)ihE b (t)j: (5.11) This represents collective dephasing in the energy eigenbasis, wherein the dephasing process randomizes the relative phase between eigenstates of the system Hamiltonian. Thus, this model does not support phase coherence between energy eigenstates. 1 As a consequence, thermal relaxation tends to equilibrate the system towards its Gibbs state, with a charac- teristic timescale set by the inverse of the bath spectral density at the gap frequency [12]. If instead the qubit system is coupled to independent, identical baths, the system-bath coupling operator becomes H ind SB =g X i z i B i ; (5.12) where, e.g., in the bosonic case B i = P k (a k;i + a y k;i ). Thermal relaxation eects occur here similarly to the collective dephasing case. However, simulations in this case are more 1 It is worth pointing out a caveat. Namely, the collective dephasing model in general supports decoherence free subspaces (DFSs), i. e., subspaces that evolve unitarily despite the coupling to the bath [74, 75]. For instance, the S = 0 subspace (for even N) is a DFS of the p-spin model. However, the p-spin model is unsuitable for performing quantum annealing inside a DFS, since its Hamiltonian consists of operators that preserve the DFS, so that no dynamics would take place if we were to try to encode a computation using states inside the DFS. Instead, to obtain meaningful dynamics (performing a computation) subject to the collective dephasing model, we would need to add Heisenberg exchange terms to the system Hamiltonian [76]. 92 demanding due to the fact that the spin symmetry is broken and that we have N times as many Lindblad operators, i. e., L ab;i (t) =hE a (t)j z i jE b (t)ijE a (t)ihE b (t)j: (5.13) Therefore, in this case we will only investigate reverse annealing starting from the rst excited state in the symmetric subspace with maximum spin, i. e., jw = 1i, for N 2 f3;:::; 8g. Moreover, for the particular cases ofN = 7 andN = 8, we truncate our system to the lowest 29 and 37 eigenstates, respectively, to speed up the numerics. This choice is made since the rst three levels of the maximum spin subspace at s = 1 are spanned by P 2 i=0 7 i = 29 (for N = 7) and P 2 i=0 8 i = 37 (for N = 8) energy eigenstates. We conrm that this is a good approximation by checking that the total population among these levels is close to 1 during the reverse annealing when additional levels are included in the simulation. The adiabatic master equation in Eq. (5.7) is unraveled using a time-dependent Monte Carlo wavefunction (MCWF) approach [27]. The advantage of MCWF is that it allows to work with wavefunctions rather than density matrices, thus saving quadratically in the dimension of the objects we need to store for numerical calculations. The tradeo is that to recover the statistical properties of the density operator we need to average over a large numberK of independent trajectories. For the collective system-bath coupling of Eq. (5.10), the time evolution operator of each trajectory is generated by the eective non-Hermitian Hamiltonian H e (t) =H(t) +H LS (t) i 2 X a6=b ab L y ab (t)L ab (t) 0 i 2 X ab L y aa (t)L bb (t); (5.14) 93 where ab and 0 are the rates for jumps and dephasing, respectively. They are related to the temperature and to the spectral density of the bosonic bath, J(!) =g 2 X k (!! k ) = 2!e !=!c ; (5.15) where! k are the bath eigenfrequencies,! c is a high-frequency cuto and is the dimension- less coupling strength. We x ! c = 1 THz and = 1 10 3 . The working temperature is chosen to be T = 12:1 mK = 1:57 GHz, as in experimental quantum annealing systems [77]. Eq. (5.14) is easily extended to the independent dephasing case by including a summation over i for the Lindblad operators and in the Lamb shift term. The time evolution operator generated by the non-Hermitian Hamiltonian of Eq. (5.14) is not unitary. Therefore, the norm of the wavefunction decays in time. Whenever the squared norm decreases below a randomly extracted threshold r2 [0; 1], a quantum jump occurs, projecting the wavefunction on one of the eigenstates of the Hamiltonian H(t). The adiabatic master equation dynamics is found after averaging over all stochastic trajectories. In this work, we xK = 5000 trajectories and consequently nd relative Monte Carlo errors hOi=hOi of the order of 1:5 % over all observables O. We repeat the simulations we reported in Sec. 5.4 for N = 20,p = 3 and = 100 ns, but now include the role of the environment. We consider both the collective and independent dephasing models. In Fig. 5.3(a), we show the success probability as a function of the inversion time s inv , for the four initial magnetizations m 0 = 0:9, 0:8, 0 and1. Monte Carlo errors are of the order of the point size in all cases and are invisible. As in the unitary case, if the inversion occurs too early (i. e., fort inv t , or, equivalently, s inv s ), the reverse annealing protocol fails to nd the ferromagnetic ground state. In fact, thermal excitations are suppressed, as well as Landau-Zener transitions, due to the large level spacing, compared with the temperature and the inverse of the annealing time. 94 For t inv t (s inv s ), however, the scenario is drastically dierent from the unitary case of Fig. 5.2(a). The rst dierence is that the success probability can be nonzero even if the inversion occurs fors inv &s , especially form 0 = 0:9, where the tail of the curve extends tos inv 0:75. When the instantaneous gap is of the same order of magnitude as the temperature, thermal processes in uence reverse annealing even before crossing the minimal gap. Second, for all m 0 we observe a sudden increase in the success probability arounds inv s , that eventually brings all curves to an almost at region at s inv <s , where the success probability reaches the large value P 0 0:957. The value of the maximum success probability at the plateau is m 0 -independent within Monte Carlo errors. The time at which the success probability starts to increase with respect to the baseline depends on m 0 . Moreover, the at region is wider for larger m 0 , although it has a nite width for all m 0 . These results show that even trial solutions far in Hamming distance from the ferromag- netic ground state can result in a large success probability at the end of a reverse anneal. Moreover, the time window in which inverting the annealing favors the ferromagnetic order- ing is relatively large. We also studied a longer annealing time, = 500 ns, as shown in Fig. 5.3(b). Here, we note that the onset of the success probability plateau shifts towards longer values of s inv , compared with the = 100 ns case. Therefore, the plateau is wider, and the maximum success probability at the plateau is P 0 1 within Monte Carlo errors for all m 0 we consid- ered. This is in contrast with the unitary case of Fig. 5.2(b), where increasing the annealing time had detrimental eects on the algorithm. This evidence supports the notion that the success probability enhancement is due to thermal eects, rather than due to purely unitary quantum dynamics [78]. Moreover, the adiabatic theorem for open quantum system guar- antees convergence to the steady state of the superoperator generator of the dynamics in the large limit [79, 80]. This too helps to explain our observations: the steady state of the Davies-Lindblad generator of the open system dynamics we considered here is the Gibbs 95 distribution of the nal Hamiltonian, which at suciently low temperature relative to the gap is the ferromagnetic ground state. Recall that in our case 2:45 GHz (ats 0:309) and T = 1:57 GHz. Comparing these results with conventional forward annealing in the presence of dissi- pation, we note that the success probability at the plateau is similar to that of standard quantum annealing for = 1 10 3 (P 0 = 0:98). The reason is that collective dephasing favors the ferromagnetic ordering in thep-spin model and its induced relaxation increases the success probability compared to the isolated case, in agreement with previous ndings [78]. We also observe that when the annealing is far from adiabaticity (e.g. = 1 ns), reverse annealing becomes a more favorable approach to forward annealing for the p-spin model. We also compare the collective and independent dephasing models of Eqs. (5.10) and (5.12). Fig. 5.4 shows the simulation results for the two models using the adiabatic master equation of Eq. (5.7) for N 2f3;:::; 8g. As shown in the gure, simulations using the collective dephasing model have larger success probabilities for almost every s inv . This is because, in the independent dephasing model, other states not in the subspace of maximum spin become accessible by thermal excitation or diabatic transition during the reverse anneal. For all of the system sizes we simulated, we had to reverse anneal to a smaller inversion point s inv for the independent dephasing model to achieve the same success probability as the collective dephasing model. Moreover, the maximum success probability achievable is always smaller for the independent dephasing model. The success probabilities from both models, however, are very similar as s inv ! 0, i. e., in the quench limit of the direct part of the evolution. Figure 5.5 shows how the maximum success probability (over s inv ) of both bath models depends on the number of qubits. As N increases, the maximum success probability of the independent dephasing model decreases more rapidly than that of the collective dephasing model. While we can infer that if we modeled independent dephasing for N = 20 we would not observe as large success probabilities as in Fig. 5.3, we stress that reverse annealing in 96 the independent dephasing model still yields a signicantly larger success probability (for the same N values) than the unitary dynamics case described in Section 5.4. 5.6 Open system dynamics with a pause Quantum annealing in the presence of a low temperature bath can benet from pauses in- serted at certain times during the dynamics [25, 81]. During a pause,s(t) = constant, and the system evolves with a time independent Hamiltonian, subject to dephasing. When a pause is inserted some time after s , the environment favors a redistribution of the repopulation according to the Gibbs state at the pause point; at suciently low temperature (relative to the gap at this point), this can result in a repopulation of the instantaneous ground state. In this section, we show that pauses at the inversion point can further improve the performance of reverse annealing of the p-spin model. We repeat the simulations for N = 20, p = 3 and = 100 ns, using the collective dephasing model. A pause of duration t p = is inserted at t = t inv , so that the total annealing time, including the pause, is 0 = +t p = 200 ns. In Fig. 5.6(a), we report the success probability as a function of the inversion point, for starting magnetizations m 0 = 0:9, 0:8, 0 and1. We compare the paused case with the unpaused case, for which = 100 ns. As can be seen in the gure, if the dynamics is reversed too early (s inv s ), the success probability at the end of the anneal vanishes. The level spacing is large compared with the temperature. The relaxation rate is small and the pause is too short to have impact on the dynamics. However, the presence of a pause signicantly changes the outcome of the annealing around s inv s . In fact, when a pause is inserted at s inv & s , the success probability reachesP 0 1 for a wide range of inversion points and for allm 0 , within Monte Carlo errors. Here, the ground state is completely repopulated by thermal relaxation. This is in contrast with conventional quantum annealing, where the success probability exhibits a peak as a 97 function of the pausing time, when the pause is inserted about 20 % later than s , and then rapidly returns to its baseline value [25, 81]. In contrast, for s inv < s , the eect of the pause is negligible; the solid (with pause) and dotted (no pause) lines in Fig. 5.6(a) overlap in this region. We repeated our analysis for a pause duration t p = 400 ns, with total annealing time 0 = 500 ns. As shown in Fig. 5.6(b), the longer pause duration aects the results only marginally. Comparing with Fig. 5.6(a), we note that the qualitative behavior of the curves is the same in the two cases. The pause duration aects mostly the region s inv & s . A longer pause enhances thermal relaxation, thus the success probability starts to increase from its baseline earlier than for shortert p . This results in a wider plateau where the success probability is large, compared with Fig. 5.6(a). Finally, we compare the collective and independent dephasing models while including pausing, starting from the rst excited state of the maximal spin sector. The results are shown in Fig. 5.7, for a pause of duration t p = 100 ns inserted at the inversion point. The collective dephasing model continues to exhibit higher success probabilities than the independent dephasing model, as in the case discussed in the previous section, but the results of the two models coincide when s inv <s . Thus, relaxation to the ground state during the pause improves performance for both dephasing models. Note that, as N increases, the maximum success probability of the collective dephasing model is achieved at s inv > s , while it is achieved at s inv <s in the independent dephasing model. This is in agreement with the N = 20 result shown in Fig. 5.3(a). 5.7 Conclusion of this chapter Earlier work revealed an intriguing tension between experimental results demonstrating a substantial enhancement in success probabilities for random spin glass instances under re- verse annealing compared to standard (forward) annealing [25], and theoretical results nding 98 that reverse annealing adversely aects performance for thep-spin model, in a closed system setting [23]. In this work we resolved this tension by performing a numerical study of reverse annealing of thep-spin model in an open system setting, where we included dephasing in the instantaneous energy eigenbasis. We found that the associated thermal relaxation results in signicant increase in the success probabilities, as long as the inversion point of the reverse annealing protocol is chosen to be close to the avoided crossing point, or before it. Pausing at the inversion point further improves performance. Since closed-system, unitary dynamics predicts that reverse annealing fails, yet its open system analogue succeeds, it follows that thermal relaxation is the mechanism responsible for the success. Reverse annealing is thus an example of a family of protocols that strictly benet from thermal eects [82, 83]. It may be worth noting that quantum eects are likely to play an important role in thermal relaxation because the success probability is very small for s inv close to 1, i.e. when the Hamiltonian stays classical during the anneal, even with a pause. Whether this can lead to any quantum speedups is an interesting problem worthy of future investigations. 5.8 Acknowledgments of this chapter This chapter was originally published in [65]. The simulation data of collective coupling model is provided by Gianluca Passarelli. The research of KY, DL, and HN is based upon work (partially) supported by the Oce of the Director of National Intelligence (ODNI), Intelligence Advanced Research Projects Activity (IARPA), via the U.S. Army Research Oce contract W911NF-17-C-0050. The views and conclusions contained herein are those of the authors and should not be interpreted as necessarily representing the ocial policies or endorsements, either expressed or implied, of the ODNI, IARPA, or the U.S. Government. The U.S. Government is authorized to reproduce and distribute reprints for Governmental purposes notwithstanding any copyright annotation thereon. Computation for some of the 99 work described in this chapter was supported by the University of Southern California Center for High-Performance Computing and Communications (hpcc.usc.edu). 100 N N N N N N Figure 5.4: Success probability in reverse annealing as a function of the inversion point s inv , forN2f3;:::; 8g. The initial state is the rst excited state of the maximum spin subspace (m 0 = 1 2=N). The dashed vertical line denotes the time s of the avoided crossing between the ground state and the rst excited state. The annealing time is = 100 ns. We sampled the interval s inv 2 (0; 1) using a step size of s = 0:005. 101 Figure 5.5: Maximum success probability achievable with reverse annealing as a function of the number of qubits N, using the collective and the independent dephasing models of Eqs. (5.10) and (5.12), respectively. The annealing time is = 100 ns. 0.0 0.2 0.4 0.6 0.8 1.0 0.0 0.2 0.4 0.6 0.8 1.0 η = 1 × 10 -3 s Δ = 0.309 τ = 100 ns t p = 100 ns Success probability s inv m 0 = 0.9 m 0 = 0.8 m 0 = 0.0 m 0 = -1.0 (a) 0.0 0.2 0.4 0.6 0.8 1.0 0.0 0.2 0.4 0.6 0.8 1.0 η = 1 × 10 -3 s Δ = 0.309 τ = 100 ns t p = 400 ns Success probability s inv m 0 = 0.9 m 0 = 0.8 m 0 = 0.0 m 0 = -1.0 (b) Figure 5.6: Success probability in paused reverse annealing for the collective dephasing model, as a function of the inversion points inv , for several values of the magnetization of the initial state. The coupling strength is = 1 10 3 . The annealing time is = 100 ns. In (a) a pause of durationt p = 100 ns is inserted at the inversion point, while in (b)t p = 400 ns. All other parameters are given in the main text. Dotted lines with empty symbols refer to open system reverse annealing of time = 100 ns and no pauses. The dashed vertical line denotes s inv =s . The interval s inv 2 (0; 1) is sampled using a step size s = 0:02. 102 N N N N N N Figure 5.7: Comparison of the success probability in reverse annealing for the collective and independent dephasing models, as a function of the inversion point s inv , for N2f3;:::; 8g. The initial state is the rst excited state of the maximum spin subspace (m 0 = 1 2=N). The dashed vertical line denotes the time s of the avoided crossing between the ground state and the rst excited state. The annealing time is = 100 ns, and a pause of duration t p = 100 ns is inserted at the inversion point. 103 Chapter 6 Experimental iterative reverse annealing In this chapter, we continue to study the open-system eects on iterative reverse annealing. We perform the IRA simulations with more realistic experimental settings of D-Wave an- nealers, forp-spin model withp = 2. We explore both quantum and semiclassical simulation methods, and also compare our simulation results to the experiments. 6.1 Introduction Iterated reverse annealing (IRA) is a variant of quantum annealing, in which the system is prepared in a trial solution state, reverse-annealed to an inversion point, and then forward- annealed. It may also be iterated with the last output state as the input of the next iteration. We perform experiments of reverse annealing on the D-Wave 2000Q device, with a focus on the p-spin problem with p = 2. We examine the dependence of the performance on param- eters such as problem size, annealing time and inversion points, as well as pausing and the number of iterations. The performance is evaluated in terms of total and partial success probabilities, the latter being the probabilities for each of two degenerate ground states. To explain the experimental results, we perform open-system simulations with decoherence in the instantaneous energy eigenbasis. The experimental total success probabilities match the simulation data well. We nd a dierence between the independent and collective dephas- ing models, which is explained by non-conservation of spin symmetry of classical inputs. 104 Simulations show that iterated reverse annealing fails in a closed system, suggesting the im- portance of thermal relaxation in the convergence to the correct solution. We also perform semi-classical simulations of spin-vector Monte Carlo to supplement our understanding of the experimental data. Recently, D-Wave devices that actually realize quantum annealing on a hardware have been developed. Not only can the D-Wave device be used to solve optimization problems by following quantum annealing physically, but it is also able to perform physics experiments, for example, spin-glass phase transitions [84], the Kosterlitz-Thouless phase transition [20, 85], the Kibble-Zurek mechanism [86], the Griths-McCoy singularity [87], spin ice [88], and the Shastry-Southernland model [89]. As in the case of simulated annealing [90], which is a classical counterpart of quantum annealing, one of the simple ways to enhance the performance of quantum annealing is to control the \annealing schedule". For conventional quantum annealing [1, 4], we choose the ground state of the transverse eld term as the initial state and dene forward annealing as the temporal change of the system from the transverse eld Hamiltonian to the target Hamiltonian. Traditionally, this transition process is supposed to be adiabatic. Diabatic transition processes have also been proposed to overcome the diculties of adiabatic com- puting, see e.g., [26, 91, 92]. On the other hand, reverse annealing is dened as the process of choosing an appropriate classical state as the initial state, starting from the target Hamil- tonian and mixing it with the transverse eld Hamiltonian to return to the original target Hamiltonian [16, 22, 23]. In chapter. 5 (Ref. [65]) we tested the eectiveness of reverse annealing on the p-spin models withp = 3. Numerical calculations using the time-dependent Lindblad equation [13, 27] showed that the presence of coupling with the boson eld signicantly changes the per- formance of reverse annealing. This result shows that the interaction with the environment can lead to an improvement in performance. Since qubits and the environment should always interact with each other on the hardware, it is likely that this result re ects real physical 105 processes in the D-Wave device. However, it is non-trivial if the coupling model to the environment of chapter. 5 is good enough to explain experimental data. In chapter. 5, the p-spin model with p = 3 is discussed because it has a rst-order phase transition in the thermodynamic limit [70, 93, 94] and thus is a dicult problem for simple quantum annealing. In the present chapter, we instead choose the p-spin model with p = 2. This is because the p-spin model with p = 2 can be more easily embedded in the hardware graph and thus may avoid possible errors by minor embedding of three-body interactions, and also larger system sizes can be implemented up to about 40. Another reason is the existence of ground-state degeneracy, which leads to additional information as will be described below. Our goal is not to study if a high success probability can be achieved but to understand the mechanism of physical processes working in the device. In this chapter, we perform reverse annealing on the p-spin model with p = 2 exper- imentally and numerically. We compare the success probability with numerical results to verify how the performance of reverse annealing on D-Wave 2000Q is enhanced by coupling to baths of boson eld. The structure of this chapter is as follows. We rst dene the p-spin model with p = 2 in Sec. 6.2 and describe the sampling details in this Chapter. In Sec. 6.3, we discuss the sampling results from the D-Wave device. We show and examine the simulation results based on closed system Schr odinger equation and quantum adiabatic master equation with two system-bath coupling models in Sec. 6.4. In Sec. 6.5, we perform semi-classical simulations of spin-vector Monte Carlo to further understand the experimental data. The chapter is concluded in Sec. 6.6. 6.2 Problem denition and protocols of experiment In this section we describe the problem Hamiltonian and the protocol of reverse annealing we used both in our experiments on the D-Wave device and in numerical simulations. 106 6.2.1 p-spin model with p = 2 The total Hamiltonian of quantum annealing is composed of the driver Hamiltonian H D and the target HamiltonianH 0 , the latter encoding the combinatorial optimization problem represented as an Ising model, as a function of dimensionless time 0s(t) 1, H(s) = A(s) 2 H D + B(s) 2 H 0 : (6.1) Here, A(s) andB(s) are device-dependent annealing schedules. The D-Wave device used in the present experiment has A(s) and B(s) as shown in Fig. 6.1(a). (a) (b) 0.0 1.0 s t [μs] FA RA+Pause t inv+t p s inv t inv 2t inv+t p 0 2 4 6 8 10 12 14 0 0.2 0.4 0.6 0.8 1 Energy [GHz] s A(s) B(s) Figure 6.1: (a) Annealing schedule of D-Wave 2000Q with the \DW 2000Q 6\ solver. (b) Forward (FA) and reverse annealing (RA) protocols are dened in Eq. (6.4)) and Eq. (6.5)), respectively. The latter incorporates a pause of duration t p . The parameter values for this plot are t a = 1:2 [s] for forward annealing, and = 1 [s] and s inv = 0:2 for reverse annealing. These symbols are dened below Eq. (6.5). The nal values of the schedule functions are A(s = 1) 0 and B(s = 1) > 0 so that ideally the ground state of H 0 is realized as the nal state, irrespective of the type of annealing protocol, i.e., the time dependence ofs(t) as discussed in the next subsection. The driver Hamiltonian H D is usually chosen as H D = N X i=1 x i (6.2) 107 where N is the number of qubits and x i is the x component of the Pauli matrix acting on the ith qubit. We study the p-spin model, H 0 =N 1 N N X i=1 z i ! p (6.3) with p = 2, in contrast to p = 3 used in chapter. 5. The main reason for the choice of p = 2 is that we can study the eects of ground-state degeneracy in p = 2 in contrast to the non- degenerate ground state of p = 3. Another technical reason is that two-body interactions can be directly embedded on the D-Wave device, which would reduce possible minor errors coming from reduction of three-body interactions to two-body representations. 6.2.2 Reverse annealing protocol The traditional protocol of forward annealing has s(t) = t t a t2 [0;t a ]; (6.4) where t a is the total annealing time (see Fig. 6.1(b), the red full line). The initial state of forward annealing is the ground state of the driver Hamiltonian H D . In reverse annealing, we start from s = 1 and decrease s to an intermediate value s inv , pause, then increase s to nish the process at s = 1. The explicit time dependence s(t) realized on the D-Wave device is s(t) = 8 > > > > > > < > > > > > > : 1 t 0tt inv 1 t inv =s inv t inv tt inv +t p 2s inv 1 + t t inv +t p <t 2t inv +t p ; (6.5) 108 where 1= is the annealing rate,t inv is the inversion time dened ast inv =(1s inv ), andt p is the length of an intermediate pause. Our reverse annealing begins from s(t = 0) = 1 and ends at s(t =t a ) = 1, where t a = 2t inv +t p = 2(1s inv ) +t p ; (6.6) passing through the inversion points inv ; see the blue dashed line in Fig. 6.1(b). To distinguish between andt a , we henceforth refer to the former as the \annealing time" and the latter as the \total annealing time". The initial state of reverse annealing is a classical state, usually a candidate solution to the combinatorial optimization problem, i.e, a state that is supposed to be close to the solution. We can iteratively repeat the process of reverse annealing with the nal state of a process fed into the system as the initial condition of the next iteration. We call this protocol iterated reverse annealing and the denote the iteration number asr. In our experiments for single-run reverse annealing, we constructed 15 instances (embeddings on the Chimera graph of the device) and generated 10 random gauges for each set of parameter values, and 1000 annealing runs were repeated for a given random gauge. Thus, the total number of runs for each set of parameter values is 150,000. In iterated reverse annealing, we constructed 15 instances and generated 90 random gauges for each r. Thus, the total number of samples at each r and each s inv is 1350. Notice that the present protocol of reverse annealing on the D-Wave device is dierent from the one adopted in Chapter. 5 ( Ref. [65]) for numerical computation. In particular, whens inv is close to 0, the second part to increases of their protocol has a very sharp increase of s as a function of t, which seems to have lead to a signicant drop of success probability fors inv 0. This is not the case in the present protocol, where the time derivative of s(t) is a constant1= (or 0 in pausing) irrespective of s inv . 109 6.3 Experimental results In this section we report on the performance of reverse annealing for the p = 2p-spin model on the D-Wave 2000Q device. 6.3.1 Dependence on initial conditions Figure 6.2 shows the experimental success probability, i.e., the ground-state probability of the nal state forN = 20 and no pause. The initial condition is specied by the initial value of the normalized total magnetization, m 0 = 1 N N X i=1 h 0 j z i j 0 i; (6.7) wherej 0 i denotes the initial wave function (a classical state). Since the doubly degenerate ground state of the problem Hamiltonian H 0 has1 as the normalized magnetization, a value of m 0 closer to 1 (or1) represents an initial condition exhibiting higher overlap with the ground state. In Fig. 6.2 we present results for a completely unbiased initial condition (m 0 = 0) and initial conditions strongly biased toward the all-up state z i = 18i (m 0 = 0:8; 0:9). The state with m 0 = 0 is the highest excited state, and the state with m 0 = 0:9 is the rst excited state. 6.3.1.1 m 0 = 0 Figure 6.2(a) shows the probability that the system reaches the all-up state (denoted \up") and the all-down state ( z i =1; 8i, denoted \down"). We call these the partial success probabilities. When the initial condition is m 0 = 0, the up and down probabilities are equal to within the error bars for all s inv , as expected because the initial state is unbiased. Note that both success probabilities are close to 0:5 for s inv < 0:4 whereas they are zero for s inv > 0:5. The vertical dashed line ats inv 0:36 indicates the minimum-gap point, i.e., the point where the 110 (a) (b) Figure 6.2: Success probability for dierent initial conditions m 0 (magnetization) in reverse annealing on the D-Wave 2000Q device as a function of the inversion points inv . The dashed line is the minimum-gap point at s 0:36 for N = 20. When s inv < s the reverse direction of the anneal goes through and past the minimum gap point, and then crosses it again during the forward anneal. When s inv >s there is no crossing of the minimum gap point. (a) Partial success probability for m 0 = 0; 0:8, and 0:9. (b) Total success probability for the same set of m 0 as in (a). Here and in all subsequent gures the labels \up" and \down" mean the populations of the all-up state and the all-down state, respectively. Error bars denote standard deviations. energy gap between the ground and the second excited state becomes minimum, which we denote bys (the rst excited state becomes degenerate with the ground state at the end of the anneal, hence it is the second excited state that is relevant). It is likely that the system remains close to its initial state for s inv > 0:5, where the transverse eld and the associated quantum uctuations are small, and the true nal ground state with m 0 =1 is hard to reach, since the initial state with m 0 = 0 has small overlap with the latter. The situation is dierent when s inv < 0:4, i.e., when the system traverses the minimum-gap region. In this case the system enters the paramagnetic (quantum-disordered) phase dominated by the driverH D , so that the initial condition is eectively erased. This renders the process similar to conventional forward quantum annealing, by which the two true ground states are reached with nearly equal probability 0:5. Recall that the phase transition around s is of second order in the present problem, and therefore the system nds the true nal ground state 111 relatively easily by forward annealing because of the mild, polynomial, closing of the energy gap, as we discuss in more detail later (see Fig. 6.5(c)). 6.3.1.2 m 0 = 0:8; 0:9: up/down symmetry breaking for 0:3.s inv . 0:5 For the initial states with m 0 = 0:8 and m 0 = 0:9, the experimental results shown in Fig. 6.2(a) reveal signicant dierences between the probabilities of the nal all-up and all- down states. Of course, the initial conditionsm 0 = 0:8; 0:9 have much larger overlap with the all-up state, and this fact alone suggest a mechanism for breaking the symmetry between the probabilities of the nal state being the all-up or all-down state. However, it is remarkable that the success probability for the all-up state is almost 1 in a region 0:3.s inv . 0:5 both to the left and to the right of s 0:36 (we discuss the additional asymmetry for s inv . 0:3 below). As we explain in Sec. 6.4.2, this behavior is inconsistent with a model of an open quantum system that is weakly coupled to its environment and is thus described by the adiabatic master equation [13]. However, it is consistent with a simple semiclassical model captured by the spin-vector Monte Carlo algorithm [95, 96], as explained in Sec. 6.5.2. 6.3.1.3 s inv & 0:5>s : freezing Note that for all three initial conditions the success probability eventually vanishes fors inv & 0:5 > s . The reason is that in all cases the system is initialized in an excited state, and remains in an excited state at the end of the anneal, since there is no mechanism for thermal relaxation to a lower energy state when s inv is large. This is a manifestation of the phenomenon of freezing [13, 97], i.e., the extreme slowdown of relaxation due to the fact that the system-bath interaction (nearly) commutes with the system Hamiltonian when the transverse eld is very small. In addition, the annealing timescale is manifestly too slow for downward diabatic transitions. This is true even in the absence of pausing, which suggests a discontinuity in the derivative of the schedule depicted in Fig. 6.1(b). Despite this, the 112 reversal of the anneal direction is apparently too slow in practice to have a non-adiabatic eect, or if diabatic transitions do occur, then they exclusively populate higher excited states. 6.3.1.4 s inv . 0:3<s : spin-bath polarization It is also noteworthy from Fig. 6.2(a) that the initial condition results in dierent probabilities for the all-up and all-down states even in the paramagnetic region s inv . 0:3 < s , where quantum uctuations are large and the all-up and all-down states are expected to have the same probability in equilibrium. The system \remembers" the initial condition to a certain extent, an anomaly that may be attributable to spin-bath polarization [98]. Namely, the persistent current owing in the qubit body during the anneal produces a magnetic eld that can partially align or polarize an ensemble of environmental spins local to the qubit wiring, with a much slower relaxation time than the anneal duration. Given the polarized initial condition m 0 = 0:8 or 0:9, this polarized spin-bath will be aligned with the all-up state even after the system crosses into the paramagnetic phase, thus preventing the system from equilibrating and explaining the observed memory eect. 6.3.1.5 Total success probability Figure 6.2(b) shows the total success probability, the sum of the nal up and down proba- bilities. The resulting curve stays almost at and close to 1 fors inv < 0:4. This constant 1 is a mixture of the two eects manifest in Fig. 6.2(a), the genuine eects of reverse annealing around s inv 0:4 for m 0 = 0:9 and 0.8 and the eectively-forward-annealing-like behavior for m 0 = 0. We shall see in Sec. 6.4.2.2 that all the salient features seen in Fig. 6.2(b), such as the shift to the left with decreasing m 0 , are captured by our open system simulations. The results seen in Fig. 6.2(b) are also consistent with the trend found numerically in Ref. [65] for the p = 3 p-spin model except for the very small values of success probability near s inv = 0 seen in the latter's numerical simulations. This may be explained by the very fast forward annealing in this parameter region in Ref. [65], which keeps the system 113 almost unchanged from the quantum-disordered state at s =s inv . Aside from this subtlety, our experimental data are consistent with the result of Ref. [65], supporting the latter's expectation that the system-environment interaction (through spin-boson coupling) prompts relaxation of the system toward the ground state around the minimum-gap region. 6.3.2 Dependence on annealing time As seen from Eq. (6.6), the total annealing time t a is linearly dependent on the annealing time , and proportional to when t p = 0. Figure 6.3 shows the success probability for dierent with N = 20, m 0 = 0:9, and no pausing. The overall trend is similar to Fig. 6.2. As shown in Fig. 6.3(a), an increase of leads to slightly higher all-up success probabilities for s inv > s but the other way around for s inv < s . This can be explained in terms of an increased relaxation to the ferromagnetic ground state near the minimum energy gap for larger and an enhanced relaxation to the paramagnetic ground state for s inv <s . (a) (b) Figure 6.3: Success probability for dierent annealing times in reverse annealing on the D-Wave 2000Q device as a function ofs inv . (a) Partial success probability for = 5; 10, and 20[s]. (b) Total success probability for the same set of as in (a). The dashed line is the minimum-gap point at N = 20, s 0:36. The total success probability stays close to 1 for s inv < s and benets slightly from increasing for s inv >s (Fig. 6.3(b)). 114 6.3.3 Eects of pausing Figure 6.4 shows success probabilities for pausing time t p 2f0; 2; 10g [s] as a function of s inv with N = 20, m 0 = 0:9 and = 5 [s]. (a) (b) Figure 6.4: Success probability for dierent pausing timest p in reverse annealing as a function of s inv . (a) Partial success probability for pause t p = 0; 2, and 10 [s]. (b) Total success probability for the same set of t p as in (a). The dashed line is the minimum-gap point at N = 20, s 0:36. The overall trend is similar to Fig. 6.3, i.e., a longer pause leads to an increased relaxation, but the eects are more signicant in the present case. Pausing at 0:5<s inv < 0:6 greatly improves the success probability from nearly 0 (t p = 0 [s]) to nearly 1 (t p = 2 [s] and 10 [s]), implying that pausing in a relatively narrow region slightly past the minimum gap point is most eective. This is in line with previous ndings [25, 96, 99]. Another signicant dierence between pausing and not pausing is that the partial success probabilities shown in Fig. 6.4(a) are at for t p > 0, showing that the eect responsible for the anomaly disappears under pausing. This is consistent with the spin-bath polarization eect discussed in Sec. 6.3.1, in that the pause provides the time needed for this polarization to relax to equilibrium, on a timescale of a few s. 115 (a) (b) 0 5 10 15 20 25 2 2.5 3 3.5 0 10 20 0.36 0.38 0.4 0.42 0.44 (c) Figure 6.5: Success probability for dierent system sizesN in reverse annealing as a function of s inv . (a) Partial success probability for system size N = 4; 8; 12; 16; 20, and 24. (b) Total success probability for the same set of N as in (a). (c) The position of the minimum gap for 4N 22, along with the value of s inv at which the success probability = :5 for each N. Inset: the minimum gap for each N. 6.3.4 Size dependence Figures 6.5(a) and (b), respectively, show the partial and total success probability for dif- ferent system sizes N as a function of s inv with = 5 [s]. The initial state for each N is the rst excited state with m 0 < 1:0 (closest to the ground state for which m 0 = 1:0), i.e., a state with one spin ipped. Although we observe a slight non-monotonicity with N in the regions inv s in Fig. 6.5(a), the overall trend is clear, namely, the drop-o to zero success 116 probability occurs at smaller s inv as N is increased. This is consistent with the reduction in s as a function of N, shown in Fig. 6.5(c), which is tracked by the value of s inv at which the success probability equals 0:5. Additionally, the asymmetry between the all-up and all-down states for s inv < s is enhanced with increasing N. This is consistent with the decreasing gap seen in Fig. 6.5(c), which should lead to slower relaxation in thermal equilibrium: at inverse temperature , the probability of an excitation is exp() times the probability of relaxation between the same pair of states, separated by a gap . I.e., the asymmetry induced by the spin-bath polarization eect is expected to take longer to dissipate for larger N due to the slower relaxation of the system. Finally, the overlap of data forN = 20 and 24 suggests that large-N eects have already converged at these sizes, i.e., that N 20 is suciently large to infer the behavior at larger N. Also, even the smallest systems withN = 4 and 8 already share qualitative features with larger systems. 6.3.5 Eects of iteration We next study the eects of iteration on reverse annealing, i.e., how the success probability depends on the number of iterationsr. Figure 6.6 shows the result forr = 1; 10; 25 and 50 at s inv withN = 20,m 0 = 0:9,t p = 0 [s], and = 5 [s]. In Fig. 6.6, the success probabilities improve in the region s inv s as the number of iterations r increases, regardless of the initial state m 0 . The relaxation to the low-energy state due to coupling to the environment should be successively induced by iterations, and the occupation of the ground state should increase. The results in Fig. 6.6 conrm this expectation. Comparing Figs. 6.6(a) and (b) with m 0 = 0:9 and Figs. 6.6(c) and (d) with m 0 = 0, we can see that the initial state with m 0 = 0:9 has a larger improvement in success probability with fewer iterationsr in the region wheres inv s . The initial state withm 0 = 0 deviates greatly from the ground state, and there are many excited states between it and the ground 117 (a) 0.0 0.2 0.4 0.6 0.8 1.0 0.0 0.2 0.4 0.6 0.8 1.0 Success Probability s inv N = 20 τ = 5 [μs] tp = 0 [μs] m0 = 0.9 up down up down up down up down r = 1 r = 10 r = 25 r = 50 0.0 0.2 0.4 0.6 0.8 1.0 0.0 0.2 0.4 0.6 0.8 1.0 Success Probability s inv N = 20 τ = 5 [μs] tp = 0 [μs] m0 = 0.9 r = 50 r = 25 r = 10 r = 1 (b) 0.0 0.2 0.4 0.6 0.8 1.0 0.0 0.2 0.4 0.6 0.8 1.0 Success Probability s inv N = 20 τ = 5 [μs] tp = 0 [μs] m0 = 0.0 up down up down up down up down r = 1 r = 10 r = 25 r = 50 (d) (c) 0.0 0.2 0.4 0.6 0.8 1.0 0.0 0.2 0.4 0.6 0.8 1.0 Success Probability s inv N = 20 τ = 5 [μs] tp = 0 [μs] m0 = 0.0 r = 50 r = 25 r = 10 r = 1 Figure 6.6: Success probability for dierent number of iterationsr in iterated reverse anneal- ing as a function of s inv . (a) and (b) are partial and total success probabilities for number of iterations r = 1; 10; 25, and 50 with the initial state m 0 = 0:9 (the rst excited state), respectively. (c) and (d) are partial and total success probabilities with the initial state m 0 = 0 (the highest excited state). The dashed line is the minimum-gap point at N = 20, s 0:36. state. Therefore, when the system state transitions from the initial state of m 0 = 0 to other states by relaxation due to coupling to the environment, those excited states become candidates. Many iterations are expected to be necessary to obtain a signicant improvement in the success probability. On the other hand, there are only a few excited states between the initial state m 0 . 1:0 and the ground state. Therefore, it is expected that only a few iterations will be sucient to make the transition to the ground state. Our results support this picture of improved performance under iterated reverse annealing. 118 6.4 Numerical simulations We next present closed and open-systems simulations and compare the results with the data presented in the previous section. We choose relatively small system sizes N = 4 and 8 to facilitate numerical computations. We note that while we can thus simulate the initial conditionm 0 = 0 (equal number of spins up and down), we are limited tom 0 = 0:5 forN = 4 andm 0 = 0:75 forN = 8 (a single spin down). However, our experimental data do not show a strong size dependence, as discussed in the previous section, so the qualitative comparison our numerical results will enable should suce to draw relevant physical conclusions. For simplicity, we focus on reverse annealing without pausing (i.e.,t p = 0) in this section. We expect that in general pausing will lead to an overall success probability in the results of open-system simulation. 6.4.1 Closed system model An analysis of the closed system case, while being unrealistic due to the absence of ther- mal eects, is instrumental in isolating the eect and importance of diabatic transitions in explaining our experimental results. The time evolution of reverse annealing in a closed system is described as follows. The state after a single iteration (cycle) is j (2t inv )i =U(2t inv ; 0)j (0)i ; (6.8) wherej (0)i is the initial state and U(2t inv ; 0) =T exp i Z 2t inv 0 H(t 0 )dt 0 (6.9) 119 is the unitary time-evolution operator andT denotes forward time ordering. At the end of r cycles, the nal statej (2rt inv )i can be expressed as: j (2rt inv )i =U(2rt inv ; 0)j (0)i ; (6.10) with U(2rt inv ; 0) = r1 Y i=0 U(2(i + 1)t inv ; 2it inv ): (6.11) Here, the nal state of the rth cycle is the initial state of the next cycle. This condition is shared by the experiments in the previous section. The solution states are doubly degenerate, and the total success probability at the end of r cycles is p(r) =jh (2rt inv )jupij 2 +jh (2rt inv )jjdownij 2 ; (6.12) wherejupi =j"i N andjdowni =j#i N . For higher eciency without loss of accuracy, we rotate the state vector into the instantaneous energy eigenbasis representation at each time step in our numerical simulations. 6.4.1.1 Dependence on system size and annealing time We initialize the state to have a single spin down, i.e., as the computational basis state j0001i (j0ij"i,j1ij#i (m 0 = 0:5)) for N = 4 andj0000001i (m 0 = 0:75) for N = 8, respectively. Note that in our simulations these are not exact eigenstates of H(1), due to a small residual transverse eld ats = 1, as in the D-Wave annealing schedule shown in Fig. 6.1. Note further that our chosen computational basis initial states are not eigenstates of S 2 = P 2fx;y;zg (S ) 2 , whereS = 1 2 P N i=1 i . Therefore the computational basis states do not lie in the subspace of any xed value of S, where S represents the total spin quantum number of N qubits. When the initial statej (0)i is a computational basis state, there is a simple upper bound on the success probability achievable in closed system reverse annealing: the 120 population of the initial state in the maximum-spin sector (see App. B.2). This upper bound explains why the following closed system results have relatively low success probabilities. We plot in Fig. 6.7 the simulation results for the total success probability for various annealing rates, subject to the D-Wave annealing schedule shown in Fig. 6.1, in the same units. We see that the success probability after a single cycle is non-negligible only when is small enough ( < 1 [ns]), in which case diabatic transitions to states with lower energies take place for s in <s . However, the success probability remains small even in this case, in any case much smaller than in our experimental results where thermal relaxation plays the dominant role. 0 0.2 0.4 0.6 0.8 1 0 0.02 0.04 0.06 0.08 0.1 0.12 0.14 (a) 0 0.2 0.4 0.6 0.8 1 0 0.01 0.02 0.03 0.04 0.05 0.06 0.07 (b) Figure 6.7: Total success probability for closed system withr = 1 (single cycle) and a single- spin down as the initial state. (a)N = 4, (b)N = 8. Here and in the subsequent gures the nanoscecond timescale is set by the energy scale of the D-Wave annealing schedule shown in Fig. 6.1. 6.4.1.2 Dependence on initial state and number of iterations We next choose the initial statej (0)i with two spins downj0011i (m 0 = 0) for N = 4 and j00000011i (m 0 = 0:5) forN = 8, respectively, i.e., the second excited states ofH 0 for these respective system sizes. 121 0 0.2 0.4 0.6 0.8 1 0 0.02 0.04 0.06 0.08 0.1 0.12 (a) 0 0.2 0.4 0.6 0.8 1 0 0.002 0.004 0.006 0.008 0.01 0.012 (b) 0 0.2 0.4 0.6 0.8 1 0 0.005 0.01 0.015 0.02 0.025 0.03 0.035 (c) 0 0.2 0.4 0.6 0.8 1 0 0.002 0.004 0.006 0.008 0.01 0.012 0.014 (d) Figure 6.8: Total success probability with two spins down as the initial state. (a) N = 4 and r = 1, (b) N = 8 and r = 1, (c) N = 4 and r = 2, (d) N = 8 and r = 2. Figure 6.8 (top row) shows that as expected, forN = 4 andr = 1, the success probability is overall lower than the case when the initial state is the rst excited state, Fig. 6.7. Only for the highest annealing rate (smallest, most diabatic) is the maximum success probability similar to the case of the rst excited state. For N = 8, the success probability is much smaller, no matter how diabatic the process is. This indirectly conrms the dominant role played by thermal eects in our experimental results. 122 We plot in Fig. 6.8 (bottom row) the results forr = 2 cycles, where we see that the success probability decreases compared to the r = 1 case. This is consistent with the conclusions of Ref. [23]. Recalling the experimental results reported in the previous section, we may therefore safely conclude that, as expected, the closed-system picture is far from the experimental reality of the D-Wave device. 6.4.2 Open quantum system model For an open system, the state after the rst cycle is (2t inv ) =V (2t inv ; 0)(0); (6.13) where V (t; 0) =T exp Z t 0 L(t 0 )dt 0 : (6.14) L(t) is the time-dependent Liouville superoperator. We use the adiabatic master equation (AME) [13], whereL(t) at time t satises d(t) dt =L(t)(t) (6.15) =i (t);H(t) +H LS (t) +D (t) : (6.16) Here,H LS (t) is the Lamb shift term andD is the dissipator. Again, the time dependence of the integrand in Eq. (7.11) is incorporated into s(t). In the weak coupling limit,D can be expressed in a diagonal form with Lindblad operators L i;! (t): 123 D[(t)] = X i X ! i (!) L i;! (t)(t)L y i;! (t) 1 2 n L y i;! (t)L i;! (t);(t) o ; (6.17) where the summation runs over the qubit index i and the Bohr frequencies ! [all possible dierences of the time-dependent eigenvalues ofH(t)]. This dissipator expresses decoherence in the energy eigenbasis [12]: quantum jumps occur only between the eigenstates ofH(t) [27]. At the end of r cyles, the nal state (2rt inv ) is expressed as: (2rt inv ) =V (2rt inv ; 0)(0); (6.18) with V (2rt inv ; 0) = r1 Y i=0 V (2(i + 1)t inv ; 2it inv ): (6.19) We consider two dierent models of system-bath coupling: independent and collective dephasing [100]. In rst case, we assume that the qubit system is coupled to independent, identical bosonic baths, with the bath and interaction Hamiltonians being H B = N X i=1 1 X k=1 ! k b y k;i b k;i ; (6.20a) H ind SB =g N X i=1 z i X k b y k;i +b k;i ; (6.20b) where b y k;i and b k;i are, respectively, the raising and lowering operators for the kth oscillator mode with natural frequency ! k . The rates appearing in Eq. (6.17) are given by i (!) = 2g 2 !e j!j=!c 1e ! ; (6.21) 124 arising from an Ohmic spectral density, and satisfying the Kubo-Martin Schwinger condi- tion [101, 102] ( i (!) = e ! i (!), = 1=T being the inverse temperature). We use g 2 = 10 3 , the cuto frequency ! c = 1 [THz], and the D-Wave device operating tem- perature T = 12:1 [mK] = 1:57 [GHz]. With the independent dephasing assumption, the Lindblad operators are: L i;! (t) = X b a=! j" a (t)ih" a (t)j z i j" b (t)ih" b (t)j; (6.22) corresponding to dephasing in the instantaneous eigenbasisfj" a (t)ig of H(t). Similarly to the closed system case, we rotate the density matrix and Lindblad operators into the instantaneous energy eigenbasis at each time step in our numerical simulations. This keeps the matrices sparse, without loss of accuracy. For N = 8 we truncate the system size to the lowest n = 18 = 1 + 1 + 8 + 8 (i.e. total number of degenerate the ground states and rst excited states at s = 1) levels out of 256. We also consider the collective dephasing model, where all the qubits are coupled to a collective bath with the same coupling strength g, which preserves the spin symmetry. In this case, the interaction Hamiltonian becomes H col SB =gS z B; (6.23) where S z = X i z i ; B = X k b y k +b k : (6.24) With this assumption, we can group together the Lindblad operators corresponding to dif- ferent qubit i into a single one: L ! (t) = X b a=! j" a (t)ih" a (t)jS z j" b (t)ih" b (t)j: (6.25) 125 The resulting number of Lindblad operators is a factor of N smaller than that of the inde- pendent system-bath coupling model. The total success probability at the end of the r cycles is p(r) =hupj(2rt inv )jupi +hdownj(2rt inv )jdowni : (6.26) Any relaxation during the reverse annealing dynamics to the global instantaneous ground state or the instantaneous rst excited state of H(t) is benecial, the latter since it be- comes degenerate with the ground statesfjupi;jdownig ofH 0 at the end of the anneal (see App. B.1). 6.4.2.1 Dependence on annealing time Figure 6.9 (top row) shows the simulation results for success probability as a function ofs inv with N = 4, various , and the initial statej0001i (m 0 = 0:5), using the independent and collective dephasing models. The independent dephasing model leads to a maximum success probability of 1 for large and smalls inv . However, the collective dephasing model leads to a maximum success probability of 1=4. The reason is the same as in closed system simulations. For the initial statej0001i, only 1=4 of its population is in the subspace of maximum quantum spin number. However, the global (instantaneous) ground state and the rst excited state of the annealing Hamiltonian H(s) both belong to the maximum-spin subspace. Since the collective dephasing model preserves the spin symmetry and the dynamics is restricted to each subspace, at most 1=4 of the initial population can be relaxed to the instantaneous ground state and the rst excited state, and reach the correct solutions at the end of the anneal (see App. B.2 for more details). It is also noteworthy that coherent oscillations are visible for = 1 [ns] (compare with Fig. 6.7(a)), but not for the larger values of we have simulated. Recall that the experimental timescale in Sec. 6.3 is on the order of a [s]. 126 0 0.2 0.4 0.6 0.8 1 0 0.2 0.4 0.6 0.8 1 (a) 0 0.2 0.4 0.6 0.8 1 0 0.05 0.1 0.15 0.2 0.25 (b) 0 0.2 0.4 0.6 0.8 1 0 0.2 0.4 0.6 0.8 1 (c) 0 0.2 0.4 0.6 0.8 1 0 0.2 0.4 0.6 0.8 1 (d) Figure 6.9: Top row: Success probability as a function of s inv with (a) independent and (b) collective dephasing (SB denotes system-bath coupling). Initial state: j0001i (m 0 = 0:5) and r = 1, a single cycle. Bottom row: Success probability with dierent . Initial state: j0011i (m 0 = 0). (a) r = 1, (b) r = 2. Comparing the simulation results in Fig. 6.9 and the experimental results in Fig. 6.3(b), we observe that they share the same main features. Namely, the success probability in- creases with; and, as increases the maximum success probability can be achieved with a larger inversion point s inv . Most notably, the total success probability also drops to zero for suciently large s inv in our simulations. This is the freeze-out eect that is well captured 127 by the adiabatic master equation (as rst reported in Ref. [13]), since thermal relaxation is suppressed when the transverse eld magnitude becomes so small that the system and system-bath Hamiltonians eectively commute. Note that the experimental results in Sec. 6.3 have a maximum success probability as high as 1, while our closed system simulations and open-system simulations with the collec- tive dephasing model have success probabilities upper bounded by some constants < 1, as already discussed. The high total success probability observed in our experiments is evidence that, as expected, the dynamics in the D-Wave device do not preserve spin symmetry. Col- lective system-bath coupling cannot explain the experimental results, leaving independent system-bath coupling as the only candidate consistent with the experiments according to our simulations. 6.4.2.2 Dependence on initial condition and the number of cycles For the initial statej0011i (m 0 = 0) and with r = 1, the independent dephasing model still gives a maximum success probability of 1 as seen in Fig. 6.9(c). The collective dephasing model (not shown) has a maximum success probability of 1=6, since the initial state has 1= 4 2 = 1=6 of its population in the maximum-spin subspace. Comparing Fig. 6.9(a) and (c), the dependence of the success probability on s inv is sim- ilar to that of the initial statej0001i. The = 1 [ns] coherent oscillations visible for the latter are more attenuated forj0011i, but this was also the case for the closed system simu- lations (contrast Fig. 6.7(a) and Fig. 6.8(a) at = 1 [ns]). The main feature distinguishing Fig. 6.9(a) and (c) is the shift to the left of the 100 [ns] curves for the initial statej0011i, i.e., as m 0 is reduced from 0:5 to 0. This is consistent with our experimental results, as can be seen in Fig. 6.2(b). In Fig. 6.9(d), we consider the dependence on the number of cyclesr, using the indepen- dent system-bath model. For r = 2, we see that the results are similar to those of a single cycle. However, we do see a small improvement in the sense of a slight shift to the right of 128 the curves with 100 [ns] compared withr = 1, which is consistent with the experimental result shown in Fig. 6.6(d). We note that this improvement was not observed in the closed system case, as can be seen by contrasting Figs. 6.8(a) and (c). Interestingly, there is also a small signature of coherent oscillations for 100 [ns]. Finally, we note that compared to the closed systems case, the results depend much less on the initial condition. 6.4.2.3 Size dependence and partial success probability In Fig. 6.10, we show the success probability for two dierent system sizes N = 4 and 8. The initial state has a single spin down. We observe that the results do not depend much on the system size, with the total success probabilities slightly larger (shifted to the right) for N = 4. This is consistent with the experimental results forN = 4 and 8 shown in Fig. 6.5(b). 0 0.2 0.4 0.6 0.8 1 0 0.2 0.4 0.6 0.8 1 Figure 6.10: Success probability as a function of s inv for independent dephasing. The initial state has one spin down and is 1[s]. The main plot shows the results for N = 4 and 8. The inset shows the partial success probabilities for N = 8. While the total success probabilities of our open-system simulations with independent dephasing are generally consistent with the experimental total success probabilities, our open-system simulations always produce symmetric partial success probabilities as shown in Fig. 6.10, in stark contrast with the experiments, where the all-up state is strongly favored 129 over the all-down state in the region around the minimum gap (see all the top panels in the gures in Sec. 6.3). Clearly, this re ects a signicant failure of the AME model. One reason for this discrepancy is that our simulations do not include a mechanism such as the polarized spin-bath which we believe explains the experimentally observed asymmetry for small s inv (see Sec. 6.3.1.4). Moreover, our adiabatic master equation simulations result in a symmetric nal all-up and all-down population since any relaxation event during the anneal (at anys) to either the global instantaneous ground statej 0 (s)i or the instantaneous rst excited statej 1 (s)i of H(s), which become degenerate at s = 1, eventually populates j 0 (s = 1)i = (j"i N +j#i N )= p 2 orj 1 (s = 1)i = (j"i N j#i N )= p 2. These two states have equal all-up and all-down populations. Similar observations are also found in Ref. [103]. 6.5 Semi-classical simulations In an eort to understand the asymmetric partial success probability observed in our exper- iments, we performed semi-classical simulations of the same problem using the spin-vector Monte Carlo (SVMC) [95] algorithm and a new variant with transverse-eld-dependent up- dates (SVMC-TF) [96]. For the p-spin problem, we replace the Hamiltonian of Eq. (7.1) by a classical Hamiltonian: H(s) = A(s) 2 N X i sin i ! B(s)N 2 1 N N X i cos i ! p : (6.27) Each qubit i is replaced by a classical O(2) spin ~ M i = (sin i ; 0; cos i ), i 2 [0;]. For the purpose of reverse annealing, we again need to specify thet dependence ofs(t). The concept of timet is here replaced by the number of Monte Carlo sweeps. we replace by a specied number of total sweeps. The total number of sweeps is then 2(1s inv ), in analogy to the total annealing time in Eq. (6.6). To simulate the eect of thermal hopping through this semi-classical landscape with in- verse temperature, we perform at each time step a spin update according to the Metropolis 130 rule. In SVMC, a random angle 0 i 2 [0;] is picked for each spin i. Updates of the spin angles i to 0 i are accepted according to the standard Metropolis rule associated with the change in energy (E) of the classical Hamiltonian. For the p-spin problem, E cannot be expressed in a simple form as in case of the Ising problem Hamiltonian [95]. In SVMC-TF, the random angle 0 i = i + i (s) is picked in a restricted range i (s)2 min 1; A(s) B(s) ; min 1; A(s) B(s) : (6.28) The goal of SVMC-TF is to restrict the angle update for A(s) < B(s), and imitate the freeze-out eect discussed in Sec. 6.3.1.3. The full SVMC and SVMC-TF algorithms for reverse annealing are summarized in App. B.3, including the expression for E. 6.5.1 Total success probability We report on SVMC and SVMC-TF simulations for N = 4 and 8. We choose the initial condition with a single spin down. In terms of angles, for N = 4 the initial angles are f0; 0; 0;g. We again use the temperature T = 12:1 [mK]. The classical analogue of the annealing time is chosen to be = 10 3 and 10 4 sweeps. In Fig. 6.11 we display the simulation results for the total success probability using SVMC and SVMC-TF. The number of samples is 10 4 . SVMC gives high total success probabilities even with large inversion pointss inv . A larges inv value means that during the whole reverse annealing process the ratio A(s)=B(s) is small. In the D-Wave device, it is expected that for smallA(s)=B(s) the dynamics freezes, which makes it dicult to reach the ground state when the initial state is excited. With the number of sweeps increased from = 10 3 to 10 4 , we observe that even for very large inversion points s inv the total success probability of SVMC can be as high as 1. This is because the angle updates in SVMC are completely random and thus with a sucient number of sweeps, it is possible for the state to ip to the correct solution. 131 0 0.2 0.4 0.6 0.8 1 0 0.2 0.4 0.6 0.8 1 (a) 0 0.2 0.4 0.6 0.8 1 0 0.2 0.4 0.6 0.8 1 (b) Figure 6.11: Total success probability for SVMC and SVMC-TF. (a) = 10 3 sweeps, (b) = 10 4 sweeps. Error bars are 2 over 10 4 samples. However, in SVMC-TF, the range of angle update is restricted for A(s)=B(s)< 1. The restricted angle updates (freeze-out eect) prevent the state from ipping to the correct solutions. Therefore the total success probability for large inversion points s inv is basically zero in SVMC-TF, regardless of the number of sweeps. This is also what we observe in the experimental data and in the adiabatic master equation simulations, as discussed in Sec. 6.4.2.1. 6.5.2 Partial success probability We compare the partial success probability obtained from SVMC and SVMC-TF in Figs. 6.12(a) and (b), respectively, for N = 4 and = 10 3 . The SVMC-TF result for N = 8 is plotted in Fig. 6.12(c). Comparing the experimental data in Fig. 6.5(a) and Fig. 6.12, we observe a satisfactory agreement with SVMC-TF, in particular the correct trend for the unequal up and down partial success probabilities. Also, for a limited region of s inv , the partial success probabilities of the all up states can be as high as 1. Classical SVMC-TF simulations also explain the D-Wave results of random systems in Ref. [96]. 132 0 0.2 0.4 0.6 0.8 1 0 0.2 0.4 0.6 0.8 1 (a) 0 0.2 0.4 0.6 0.8 1 0 0.2 0.4 0.6 0.8 1 (b) 0 0.2 0.4 0.6 0.8 1 0 0.2 0.4 0.6 0.8 1 (c) Figure 6.12: Partial success probability by (a) SVMC, N = 4and (b) SVMC-TF, N = 4. = 10 3 sweeps. (c) SVMC-TF,N = 8. = 10 3 sweeps. Error bars are 2 over 10 4 samples. For N = 8, we observe a deviation from the experimental data shown in Fig. 6.5(a) for s inv . 0:3, where there exists a small but clear dierence in the probabilities of all-up and all-down states, whereas the SVMC-TF data do not show such a trend. As discussed in Sec. 6.3.1.4, we attribute the dierence for s inv . 0:3 to the spin-bath polarization eect, which is not modeled in our SVMC-TF simulations. In Fig. 6.13 we display SVMC-TF reverse annealing simulation results of partial success probabilities for N = 16; 32 with = 10 3 ; 10 4 sweeps. For both sizes shown, the regime 133 of high partial success probability for the all-up state is shifted slightly to higher s inv for = 10 4 sweeps then for = 10 3 . This is consistent with the trend in the experimental results observed in Fig. 6.3(a). However, the trend with system size is inconsistent with the experimental results shown in Fig. 6.5(a): the results for N = 16 and N = 32 are virtually indistinguishable (apart from statistical uctuations), which the experimental data shows that N = 16 is not yet large enough for convergence. This (small) failure of the SVMC- TF model may hint at an interesting way in which to identify a \quantum signature" in experimental quantum annealing [47, 103]. 6.6 Conclusion of this chapter The previous chapter (Ref. [65]) has reported that the coupling with the environment im- proves the performance of reverse annealing for the p-spin model withp = 3. In the present study, we have performed reverse annealing on a real device D-Wave 2000Q for the p-spin model with p = 2 and have conrmed that the reverse annealing process yields high success probabilities as in the previous chapter (Ref. [65]). We perform closed system simulations, and open-system simulations of quantum adia- batic master equation, which take into account the eect of decoherence in the instantaneous energy eigenbasis. The numerical simulations reveal that thermal relaxation indeed leads to high total success probabilities, which also match well with our experimental results. Ex- perimental data of reverse annealing also in turn helps us understand the noise mechanism in D-Wave machines. Based on the experimental data, we make careful analysis on the simulation methods and the correct choice of system-bath coupling model. We argue that reverse annealing would fail in a closed system setting or open system simulation with the assumption of collective dephasing. To further support the experimental results of asymmetrical partial success probabilities, we perform semi-classical simulations of spin-vector Monte Carlo. Our results show that 134 0 0.2 0.4 0.6 0.8 1 0 0.2 0.4 0.6 0.8 1 (a) 0 0.2 0.4 0.6 0.8 1 0 0.2 0.4 0.6 0.8 1 (b) 0 0.2 0.4 0.6 0.8 1 0 0.2 0.4 0.6 0.8 1 (c) 0 0.2 0.4 0.6 0.8 1 0 0.2 0.4 0.6 0.8 1 (d) Figure 6.13: Partial success probability for SVMC-TF. (a) N = 16, = 10 3 sweeps, (b) N = 16, = 10 4 sweeps, (c) N = 32, = 10 3 sweeps, (b) N = 32, = 10 4 sweeps. Error bars are 2 over 10 4 samples. the semi-classical thermal eect also benets the performance of reverse annealing in real- world quantum annealing devices. We hope that this realization would be a step forward in the understanding and development of devices where thermal eects lead to quantum acceleration. 135 6.7 Acknowledgments of this chapter All experimental data is provided by Yuki Bando. The experimental section is written by Yuki Bando. We thank Masayuki Ohzeki for useful comments and Sigma-i Co., Ltd. for providing machine time on D-Wave 2000Q. The research is based partially upon work sup- ported by the Oce of the Director of National Intelligence (ODNI), Intelligence Advanced Research Projects Activity (IARPA) and the Defense Advanced Research Projects Agency (DARPA), via the U.S. Army Research Oce contract W911NF-17-C-0050. The views and conclusions contained herein are those of the authors and should not be interpreted as nec- essarily representing the ocial policies or endorsements, either expressed or implied, of the ODNI, IARPA, DARPA, or the U.S. Government. The U.S. Government is authorized to reproduce and distribute reprints for Governmental purposes notwithstanding any copyright annotation thereon. Computation for some of the work described in this chapter was sup- ported by the University of Southern California Center for High-Performance Computing and Communications (hpcc.usc.edu). 136 Chapter 7 Adiabatic reverse annealing We now move on to the study of another variant of reverse annealing called adiabatic reverse annealing, originally introduced in [16, 22, 23], in an open system setting. We focus on the p-spin Hamiltonian with p = 3. 7.1 Introduction and terminology The annealing Hamiltonian for adiabatic reverse annealing (ARA) is [23, 26]: H(s;) =sH 0 + (1s)(1)H init + (1s)V TF : (7.1) The time-dependent parameterss(t) and(t) both change from 0 to 1 as timet proceeds from 0 to , where is the total anneal time. We focus on the case of a linear function s(t) =t=. The notation of Eq. (7.1) is as follows. H 0 =N 1 N N X i=1 z i ! p (7.2) 137 is the target problem Hamiltonian, with p = 3 and N the number of qubits. H init = N X i=1 i z i (7.3) is the initial Hamiltonian, i =1. V TF = N X i=1 x i (7.4) is the transverse eld Hamiltonian. We want to perform (adiabatic) reverse annealing with various initial statesj (s = 0)i =j 1 ; 2 ;:::; N i (the ground state of H init ). We further assume that =s. With this assumption, Eq. (7.1) is simplied as: H(s) =sH 0 + (1s) 2 H init + (1s)sV TF : (7.5) In the special case that(t) = 1, Eq. (7.1) is the quantum annealing (QA) Hamiltonian. 7.2 Preliminary I: Adiabatic reverse annealing (ARA) Hamiltonian spectrum In this section, we want to explore how the spectrum changes by tuning certain parameters of Eq. (7.5). This is important in the calculations of diabatic transition rates, excitation and relaxation rates, and thus the understanding of open-system behavior of adiabatic reverse annealing (ARA). The solution state of the p-spin problem withN spins isj"i N . Consider the initial state with the number of spin-downs N # = 2. We plot the spectrum for N = 10 in Fig. 7.1. 138 0 0.2 0.4 0.6 0.8 1 -10 -8 -6 -4 -2 Figure 7.1: The 16 lowest lying energies of 10 qubits, with H init having N # = 2. 7.2.1 Controlling the parameters: The energy gap between instantaneous eigenenergies i (s) and j (s) is ij (s) = i (s) j (s). The corresponding minimum energy gap is ij = min s ij (s). In adiabatic reverse annealing, we are interested in the minimum energy gap between the ground state 0 (s) and the 1st excited state 1 (s). It is denoted by (= min ). = 10 = min s 10 (s) = min s 1 (s) 0 (s): (7.6) Consider the initial state withN " = 1. We want to explore how the gap properties change with in Eq. (7.5). We calculate and plot in Fig. 7.3 the value of 10 (s), for N = 8 and =f1; ; 5g. In the inset we plot the dependence of on . In general decreases as decreases. Also, the shape of the gap is sharper with smaller . This allows diabatic transitions. 139 0 0.2 0.4 0.6 0.8 1 0 1 2 3 4 5 0 2 4 6 0.5 1 1.5 Figure 7.2: (s) for =f1; ; 5g in units where h = 1. 7.2.2 Controlling the parameters: Initial states i or (H init ) We want to investigate how the gap properties change with dierent initial states. In this case we focus on = 1. 0 0.2 0.4 0.6 0.8 1 0 1 2 3 4 5 Figure 7.3: (s) for dierent i . = 1 in units where h = 1. 140 From Fig. 7.3, we see that in general decreases as the Hamming distance between the initial state and the target state increases. 7.2.3 Controlling the parameters: N We want to investigate how the minimum gap changes with system size N. In this case we focus on the initial state of N # = 1 and = 1. The value of for each N is plotted in Fig. 7.4. We also plot the s value where the minimum gap is found, i.e., s min = argmin 10 (s): (7.7) This is important for the calculation of relaxation rates for reverse annealing. 0 5 10 15 20 0.34 0.36 0.38 0.4 0.42 0.44 0 10 20 0.025 0.03 0.035 Figure 7.4: forN =f3; ; 20g in units whereh = 1. Inset: the value ofs min for eachN. 7.3 Preliminary II: Time evolution Letj (s = 0)i (j (t = 0)i) denote the initial state. The anneal time is . The nal state at the end of the adiabatic reverse annealing (ARA) and the success probability can be expressed as follows. 141 7.3.1 Closed system The state after the ARA in a closed system is j ()i =U(; 0)j (0)i ; (7.8) where U(; 0) =T exp i Z 0 H(t 0 )dt 0 (7.9) is the unitary operator andT denotes the forward time-ordering. The time-dependence of the Hamiltonian (in the integrand) is given by the s(t) function. 7.3.2 Open system: Lindbladian dynamics For an open system, the state after ARA is () =V (; 0)(0); (7.10) where V (; 0) =T exp Z 0 L(t 0 )dt 0 : (7.11) L(t) is the time-dependent Liouville superoperator. For example, in adiabatic master equa- tion,L(t) at time t takes the form of d(t) dt =L(t)(t) (7.12) =i (t);H(t) +H LS (t) +D (t) : (7.13) H LS (t) is a Lamb shift term and D is the dissipator superoperator. Again, the time- dependence of the integrand in Eq. (7.11) depends on the s(t) function. 142 7.4 Preliminary III: Success metrics and classes of problems We use two success metrics for ARA: success probability and time to solution (TTS). 7.4.1 Success probability The success probability refers to the probability of the nal state being the target solution state. For the example of p-spin (p = 3) model, the success probability at the end of ARA can be expressed as: Successprobability p g () = 8 > > < > > : j N h"j ()ij 2 ; in a closed system. N h"j()j"i N ; in an open system. (7.14) 7.4.2 Time to solution (TTS) The time to solution (TTS) is dened as [104]: TTS(;p d ) = log(1p d ) logp e () : (7.15) p d is the desired threshold probability and p e () = 1p g (). 7.4.3 Classes of problems We want to explore four general classes of problems: 1. Adiabatic reverse annealing in an open system (ARA Open ), 2. Adiabatic reverse annealing in a closed system (ARA Closed ), 3. Standard quantum annealing in an open system (QA Open ), and 4. Standard quantum annealing in a closed system (QA Closed ). The comparison between ARA Closed and QA Closed has been made in [23] and it was shown that ARA Closed can have an advantage over QA Closed . Since there are many studies (e.g. [78] and [65]) showing that dissipation can give quantum enhancement, we want to investigate further if and under what conditions ARA Open performs 143 better than ARA Closed . On an equal setting, we want also to compare if ARA Open has any computational advantage over QA Open , in order to support the use of the formal protocol in a realistic experimental setting. Any advantage/enhancement is dened according to the two success metrics mentioned before. We have already studied how the spectrum depends on , H init and N. We want to evaluate how the performance of ARA and QA depends on these parameters. For each set of such parameters, we performed simulations in a range of anneal times . As noted in [26], diabatic transition to higher excited states may be a shortcut of the route to the nal solution. Therefore, we include in our studies anneal times shorter than the one set by adiabatic condition [3]. We also explore the cases with small and large Hamming distance between the initial state and the target state, which result in a very small and sharp gap and thus diabatic transitions to higher excited states. 7.5 Illustrative examples We perform numerical experiments for both closed and open system. We use the full Hamil- tonian without sector decomposition. The coupling is = 1e 3 and the cuto frequency is ! c = 1THz. We use an independent dephasing bath H ind SB = g P i z i B i and a collec- tive dephasing bath H col SB = g P i z i B, and perform simulations of the adiabatic master equation Eq. (7.12). We begin by studying three examples in the following subsections (Sec. 7.5.1- 7.5.3) and illustrate the dynamics of ground state population during the anneal. 7.5.1 Results of 8-qubit simulation: = 40ns. = 1. We simulate the closed and open-system evolution with dierent initial states. In this exam- ple we have N = 8, = 40ns. = 1. The results are plotted in Fig. 7.5. The solid lines are open-system simulation results while the dotted lines are closed-system simulation results. 144 0 10 20 30 40 0.88 0.9 0.92 0.94 0.96 0.98 1 Ground state population Figure 7.5: Open-system and close-system adiabatic reverse annealing (ARA) simulation results. = 40ns. = 1. Conclusion: From Fig. 7.5, we see that open-system eect (thermal excitation) is detri- mental in this case, i.e., for all initial states specied, the solid lines are below the dotted lines along the evolution and at the end of the evolution. 7.5.2 Results of 4-qubit simulation: = 4ns. = 1. In this example we have N = 4, = 4ns, and = 1. The results are plotted in Fig. 7.6. The solid lines are open-system simulation results while the dotted lines are closed-system simulation results. Conclusion: From Fig. 7.6, we see that open-system eect makes no dierence, i.e. for all initial states specied, the solid lines almost overlap with the dotted lines along the evolution and at the end of the evolution. The total anneal time = 4ns is too short for the environment (bath) to come into eect. The environment causes little dierence to the success probabilities. This short anneal time also leads to observable diabatic transitions to higher excited states. 145 0 1 2 3 4 0 0.2 0.4 0.6 0.8 1 Ground state population Figure 7.6: Open-system and close-system adiabatic reverse annealing (ARA) simulation results. = 4ns. = 1. 7.5.3 Results of 4-qubit simulation: = 400ns. = 0:3. In this example we have N = 4, = 400ns. = 0:3. The results are plotted in Fig. 7.7. Conclusion: From Fig. 7.7, we observe sudden diabatic transitions, since = 0:3 results in very small gaps for most of theH init . Meanwhile, a relatively long anneal time of = 400ns allows open-system relaxation mechanism to increase the (instantaneous and nal success probability) ground state population for some of the initial states. 7.6 ARA: Random initial state c = 0:5. As in previous studies, the success of ARA depends on whether our initial state is close to the correct solution or not. In this section, we study the most general case where the initial state is a random guess (c = 0:5, where c = N " =N) and want to see if ARA Open can give an enhancement over ARA Closed . Instead of looking at ground state population along the 146 0 100 200 300 400 0 0.2 0.4 0.6 0.8 1 Ground state population Figure 7.7: Open-system and close-system adiabatic reverse annealing (ARA) simulation results. = 400ns. = 0:3. anneal, we now focus only on the results of ARA at the end of anneal. A summary of the results as a function of N and is shown in Fig. 7.8(d). First of all, we observe that success probabilities go up monotonically with anneal time for the case of ARA Closed . As expected, the longer the anneal time , the more adiabatic the evolution and the higher the success probability. For ARA Open , however, we observe that there can exist an intermediate anneal time that gives a locally maximal success probability, and this is true for both SB coupling models. Such non-monotonic dependence of success probabilities in terms of anneal times is due to the environmental induced relaxation after diabatic transitions to higher excited states. Secondly, at long anneal times, for the case of ARA Closed the ground state is reached with with certainty (success probability = 1); while for the case of ARA Open , the Gibbs state is reached. Therefore, in all situations ARA Closed always outperforms ARA Open for long enough anneal times . However, as shown in later section 7.8, the steady state reached by the two dierent system bath coupling model is dierent. 147 10 0 10 2 10 4 0 0.2 0.4 0.6 0.8 1 (a) 10 0 10 2 10 4 0 0.2 0.4 0.6 0.8 1 (b) 10 0 10 2 10 4 0 0.2 0.4 0.6 0.8 1 (c) 10 1 10 2 10 3 500 1000 1500 2000 2500 3000 (d) Figure 7.8: ARA performance for (a) closed system, (b) open system with independent system-bath coupling, (c) open system with collective system-bath coupling, (d) comparison of the 3 situations in terms of TTS (N = 10). To see if relaxation can provide any quantum enhancement, we do observe that for some N ARA Open performs better than ARA Closed in the cases of intermediate anneal times . We demonstrate in Fig. 7.8(d) for the case ofN = 10 thatARA performs better in an open system settings in terms of shorter TTS. For both ARA Closed and ARA Open we observe optimal TTS. To avoid the misinterpretation of TTS, we avoid the calculation of TTS() where the success probabilityp g () is very close to 0 and 1. Lastly, we can see that in general the performance 148 of both ARA Closed and ARA Open worsens asN increases, as expected due to the closure of the gap with increasing system size. 7.6.1 QA 10 0 10 2 10 4 0 0.2 0.4 0.6 0.8 1 (a) 10 0 10 2 10 4 0 0.2 0.4 0.6 0.8 1 (b) Figure 7.9: Standard QA in an open system with two dierent system-bath coupling model. We show in Fig. 7.9 the case of quantum annealing in open system ARA Open . Note that the initial state of QA is the equal superpositionj+i N . Comparing Fig. 7.9 with Fig. 7.8(d), we can see easily thatQA Open achieves higher success probabilities for allN and. Therefore, for random initial states c = 0:5, ARA Open cannot outperform QA Open . 7.7 ARA performance in terms of and H init We now want to study the dependence of ARA performance in terms of H init and . For this purpose, we x the number of qubits N in our studies in this section. 149 0 1 2 3 4 5 0 0.5 1 1.5 2 Figure 7.10: The value of minimum gap as a function of H init and . N = 4. 7.7.1 Min. Gap vs H init vs . N = 4. We rst plot in Fig. 7.10 the value of minimum gap of ARA (N = 4) as a function of initial state H init and . Except for initial states that are close to the correct solutions, in general the larger the transverse eld the larger the minimum gap min . 7.7.2 Dierent Anneal times + Dierent = f1; 2; 5g. N = 4. ARA Open We simulate the performance ofARA Open as a function of initial stateH init and , and plot the results in Fig. 7.11 . We focus on the assumption of Independent system-bath coupling. In general the closer the initial state to the correct solution, the higher the success probabilities of ARA Open . Except for initial states that are close to the correct solutions, in general we want to increase the transverse eld to enhance ARA performance because it results in a 150 10 -1 10 0 10 1 10 2 10 3 0 0.2 0.4 0.6 0.8 1 (a) 10 -1 10 0 10 1 10 2 10 3 0 0.2 0.4 0.6 0.8 1 (b) 10 -1 10 0 10 1 10 2 10 3 0 0.2 0.4 0.6 0.8 1 (c) Figure 7.11: ARA open results (a) = 1, (b) = 2, (c) = 5. Independent system-bath coupling. larger gap, in agreement with the observation in Fig. 7.10. We also plot this dependence on in Fig. 7.12 in terms of success probability and TTS. TTS is also shorter for larger and we observe an optimal TTS for all three values of . 151 10 -1 10 0 10 1 10 2 10 3 0 0.2 0.4 0.6 0.8 1 (a) 10 -1 10 0 10 1 10 2 10 3 10 0 10 1 10 2 10 3 (b) Figure 7.12: ARA open results, with initial state beingj##""i. (a) success probability, (b) TTS. Independent system-bath coupling. 7.7.3 Dierent Anneal times + Dierent < 1. N = 4. ARA Closed and ARA Open We study the case of < 1, which would results in a smaller gap (Fig. 7.10) except for the case where the initial state is equal to the correct solution. We plot in Fig. 7.13 and Fig. 7.14 the simulated performance of ARA Closed and ARA Open respectively as a function of initial state H init and . Again we focus on the assumption of independent system-bath coupling. Comparing the results of ARA Closed (Fig. 7.13) and ARA Open (Fig. 7.14), we again observe the non-monotonic dependence of success probabilities in terms of anneal times exclusively for ARA Open . We also see noticeable improvements of success probabilities due to open- system eect at intermediate anneal times, especially if the initial states are far away from the correct solution. For both ARA Closed and ARA Open , the larger the higher the success probabilities, except for the case where the initial state is equal to the correct solution. 152 10 0 10 2 10 4 0 0.2 0.4 0.6 0.8 1 (a) 10 -1 10 0 10 1 10 2 10 3 0 0.2 0.4 0.6 0.8 1 (b) Figure 7.13: ARA Closed (a) = 0:1, (b) = 0:5. N = 4. 10 0 10 2 10 4 0 0.2 0.4 0.6 0.8 1 (a) 10 -1 10 0 10 1 10 2 10 3 0 0.2 0.4 0.6 0.8 1 (b) Figure 7.14: ARA Open . (a) = 0:1, (b) = 0:5. N = 4. Independent system-bath coupling. 7.7.4 Dierent Anneal times. = 0:1, N # = 2. N = f8; 10g. ARA Open . We study the extreme case of = 0:1 for larger N. Note that for such small = 0:1, the transverse eld is very weak and the Hamiltonian is dominated by the problem Hamiltonian H 0 andH init . The resulting gap is very small and sharp and there is a large degree of diabatic transition to higher excited states. 153 We plot in Fig. 7.15 the TTS dependence of anneal times , for N =f8; 10g both with N # = 2. Due to the signicant diabatic transition to higher excited states, ARA Open tends to perform better than ARA Closed because of the relaxation to the ground state. 10 0 10 1 10 2 10 3 10 5 10 6 (a) 10 0 10 1 10 2 10 3 10 5 10 6 (b) Figure 7.15: (a) N = 8, (b) N = 10. = 0:1 and N # = 2. Independent system-bath coupling. 7.8 Independent vs collective SB coupling We investigate the dierence between the independent and collective system-bath coupling assumption. 7.8.1 Long anneal times As observed in Fig. 7.8(d), for long anneal times the success probabilities given by simula- tion with the collective SB coupling assumption are higher than those given by independent SB coupling. In fact, the steady state solution (s = 1;!1) of the Lindbladian master equation in the weak coupling limit is the Gibbs state set by the problem Hamiltonian. (s = 1;!1) = e H S (1) Z = e H 0 Z : (7.16) 154 whereZ = Tr (exp (H 0 )). Recall that the permutation symmetry allows us to write the Hamiltonian in terms of only a set of total spin operators of two subsystems, ^ S x;z 1 1 2 jN " j X i=1 ^ x;z i ; ^ S x;z 2 1 2 N X i=jN " j+1 ^ x;z i : (7.17) For collective SB coupling model, the dynamics occurs entirely in the two subspaces with maximum eigenvalues of S 2 1 and S 2 2 . The Lindbladian is dened according to the projection of H(s) to such subspaces. For collective SB coupling model, the steady state is given by: 0 (s = 1;!1) = e H 0 0 Z 0 (7.18) whereZ 0 = Tr (exp (H 0 0 )) andH 0 0 =PH 0 P . P = P j ih j andj i =j1;w 1 i j2;w 2 i is the basis that has the maximum eigenvalues of S 2 1 and S 2 2 . A simple proof below shows that for long anneal times the success probability given by collective SB coupling model is always higher than the one given by independent SB coupling model. Proof. Since H 0 0 is the projection of H 0 in the subspace with the largest values of ( ^ S 1 ) 2 and ( ^ S 2 ) 2 ,Z 0 <Z. The success (ground state) probabilityhgj 0 (s = 1;!1)jgi = e eg Z 0 > e eg Z =hgj(s = 1;!1)jgi, where e g is the ground state energy. 7.8.2 Intermediate anneal times In general, since the assumption of independent SB coupling allows decoherence to many more eigenstates of H(s) outside the maximal spin subspaces, it gives a lower success prob- ability than the independent SB model for almost all anneal times . For the example of c = 0:5 and = 10ns, we see in Fig. 7.16 that the dierence of the success probability set 155 by the two models diverges as system sizeN increases. This is expected since the dimension reduction due to the restriction to maximal spin subspaces also scales exponentially in N. 2 4 6 8 10 0 0.2 0.4 0.6 0.8 1 Figure 7.16: Dependence of success probability as a function of N due to the two dierent SB coupling model. c = 0:5, = 1 and = 10ns. 7.9 Conclusion of this chapter We studied the performance of ARA as a function of transverse eld , initial state H init and system size N. In general, except for the cases where the initial state is equal to the correct solution, ARA performs better with larger . We observe non-monotonic dependence of success probabilities in terms of the anneal times for ARA Open . For long anneal times, ARA Closed reaches the solution with certainty, while ARA Open reaches the Gibbs state. For ARA Open , the dierence due to independent system-bath coupling assumption and collective system-bath coupling assumption is studied 156 for both long anneal times and intermediate anneal times . In general, the former case always results in a lower success probability (or longer TTS). ARA Open tends for outperform ARA Closed in the following conditions: 1. Small . 2. Initial state far away from the correct solution. 3. Large N. 4. Intermediate anneal times . The rst three conditions result in a small and sharp gap for diabatic transitions for in- termediate anneal times . Because of the environmentally-induced relaxation back to the ground state after diabatic transitions to higher excited states, ARA Open performs better than ARA Closed with the above conditions, in agreement with the prospects of quantum en- hancement for diabatic annealing [26]. 157 Chapter 8 Conclusion In the rst part of this dissertation, we studied time-dependent quantum trajectory tech- niques which can be applied to three dierent areas: 1) simulation of the adiabatic master equation, 2) simulation of slow noise with 1=f power spectrum and a time-dependent system Hamiltonian, 3) simulation of weak measurements and feedback error correction protocols in quantum annealing. For the trajectories simulation of the adiabatic master equation, we derived the time- dependent quantum trajectory theory and a corresponding numerical method, and performed simulations of several examples of problem gadgets with D-Wave annealers settings and rel- atively large problem sizesN. We also showed that this technique can provide up to a factor N advantage over directly solving the master equation. The cost of running many trajecto- ries, which is required to recover the master equation evolution, can be further minimized by running the trajectories in parallel. The quantum trajectories method provides insight into individual quantum jump trajectories and their statistics, thus shedding light on open system quantum adiabatic evolution beyond the master equation. For the simulation of 1=f noise in QA, based on the previously studied Bloch vector-magnetic eld approach, we formulated a stochastic uctuator Hamiltonian simulation method that includes a series of uctuator terms with dierent ipping frequencies and Gaussian couplings that altogether constitute the 1=f power spectrum. We, for the rst time, showed how to incorporate them into the inherently time-dependent annealing Hamiltonians, and construct an overall time-dependent 158 stochastic Hamiltonian for parallelizable simulations. We studied in depth how noise with 1=f power spectrum aects quantum annealing, and we also studied temperature eects, loss of coherence (T 1 and T 2 ) in QA under dierent uctuator parameters. We allowed the noise uctuators to take any form of multi-level operators and any noise-axis direction de- pending on the superconducting qubits and circuit decoherence process. We showed how to append these operators to two dierent superconducting circuit Hamiltonians. We t our simulation results with experimental data, and showed how dynamical decoupling can ex- tend the coherence of qubits under 1=f noise. For the simulation of weak measurement and feedback error correction protocols in quantum annealing, we derived the feedback master equation provided that the feedback is Markovian (feedback delay ! 0) and further gave the timescale condition for feedback Markovianity. For realistic feedback which is subjected to non-negligible feedback delay and restrictions of the form of the feedback Hamiltonian, we studied the eectiveness of feedback error correction under such limitations and showed how to optimize the feedback delay time depending on the annealing schedule. In the second part of this dissertation we studied the open-system behavior of reverse annealing and showed how it can benet from environment-induced relaxation. We presented a numerical study of the iterative reverse annealing (IRA) protocol applied to the p-spin model (p = 3), including pausing, in an open system setting accounting for dephasing in the energy eigenbasis which results in thermal relaxation. We considered both independent and collective dephasing and demonstrated that in both cases the open system dynamics can substantially enhance the performance of reverse annealing. We also studied the p- spin problem with p = 2 with D-Wave annealer settings and examined the dependence of the performance on parameters such as problem size, anneal time, inversion points, etc. We performed open-system simulations that match the experimental data well. We found that only the independent dephasing model gives a correct description of reverse annealing in noisy D-Wave annealers, due to the spin symmetry breaking of classical inputs. Semi- classical simulations of spin-vector Monte Carlo also supplement our understanding of the 159 experimental data well. For adiabatic reverse annealing (ARA), we performed simulations of the p-spin model (p = 3) with an extra initialization Hamiltonian term. We found the conditions that result in a small and sharp gap for diabatic transitions with intermediate anneal times. Under such conditions, ARA performs better in an open-system setting, due to the environmentally-induced relaxation back to the ground state after diabatic transitions to higher excited states. There is still much to understand about the topics we have studied, for example, the simulation of slow noise in superconducting qubits. First, we do not know how to simulate slow noise with 1=f spectrum, with 6= 1, in particular using the stochastic Schr oodinger equation approach of time-dependent annealing and superconducting qubit Hamiltonian. The steady state of a qubit after Landau-Zener transition under 1=f telegraph noise is non- trivial [105] and depends on the coupling strength and switching rate of the uctuators. It would be interesting to come up with a similar observation in the case of quantum annealing and many qubits. As for weak measurements and feedback error correction, there are also questions worth further research, such as: how to perform the protocol for more than 1 qubit, whether we should apply feedback based on the global excited states or the local excited states of each qubit, and how to adapt the protocol beyond the weak-coupling limit. 160 References 1. Kadowaki, T. & Nishimori, H. Quantum annealing in the transverse Ising model. Phys. Rev. E 58, 5355{5363 (5 1998). 2. Farhi, E., Goldstone, J., Gutmann, S., Lapan, J., Lundgren, A. & Preda, D. A quantum adiabatic evolution algorithm applied to random instances of an NP-complete problem. Science 292, 472{475 (2001). 3. Albash, T. & Lidar, D. A. Adiabatic quantum computation. Reviews of Modern Physics 90, 015002 (2018). 4. Santoro, G. E., Marto n ak, R., Tosatti, E. & Car, R. Theory of Quantum Annealing of an Ising Spin Glass. Science 295, 2427{2430 (2002). 5. Das, A. & Chakrabarti, B. K. Colloquium: Quantum annealing and analog quantum computation. Rev. Mod. Phys. 80, 1061{1081. doi:10.1103/RevModPhys.80.1061 (3 2008). 6. Morita, S. & Nishimori, H. Mathematical foundation of quantum annealing. Journal of Mathematical Physics 49, 125210. doi:10.1063/1.2995837 (2008). 7. Hauke, P., Katzgraber, H. G., Lechner, W., Nishimori, H. & Oliver, W. D. Perspectives of quantum annealing: Methods and implementations. Rep. Prog. Phys. 83, 054401 (2020). 8. Novikov, S., Hinkey, R., Disseler, S., Basham, J. I., Albash, T., Risinger, A., et al. Exploring more-coherent quantum annealing in 2018 IEEE International Conference on Rebooting Computing (ICRC) (2018), 1{7. 9. Amin, M. H. S., Averin, D. V. & Nestero, J. A. Decoherence in adiabatic quantum computation. Phys. Rev. A 79, 022107 (2009). 10. Smirnov, A. Y. & Amin, M. H. Theory of open quantum dynamics with hybrid noise. New Journal of Physics 20, 103037 (2018). 11. Davies, E. & Spohn, H. Open quantum systems with time-dependent Hamiltonians and their linear response. Journal of Statistical Physics 19, 511{523 (1978). 12. Albash, T. & Lidar, D. A. Decoherence in adiabatic quantum computation. Physical Review A 91, 062320 (2015). 13. Albash, T., Boixo, S., Lidar, D. A. & Zanardi, P. Quantum adiabatic Markovian master equations. New J. of Phys. 14, 123016. doi:10.1088/1367-2630/14/12/123016 (2012). 14. Amin, M. H. & Averin, D. V. Macroscopic resonant tunneling in the presence of low frequency noise. Physical review letters 100, 197001 (2008). 15. Campos Venuti, L., Albash, T., Marvian, M., Lidar, D. & Zanardi, P. Relaxation vs. adiabatic quantum steady state preparation: which wins? arXiv e-prints, arXiv{1612 (2016). 161 16. Perdomo-Ortiz, A., Venegas-Andraca, S. E. & Aspuru-Guzik, A. A study of heuristic guesses for adiabatic quantum computation. Quantum Information Processing 10, 33{ 52. doi:10.1007/s11128-010-0168-z (2011). 17. Albash, T. & Lidar, D. A. Adiabatic Quantum Computing. arXiv:1611.04471 (2016). 18. Albash, T. & Kowalsky, M. Diagonal catalysts in quantum adiabatic optimization. Physical Review A 103, 022608 (2021). 19. Chancellor, N. Modernizing quantum annealing using local searches. New Journal of Physics 19, 023024. doi:10.1088/1367-2630/aa59c4 (2017). 20. King, A. D., Carrasquilla, J., Raymond, J., Ozdan, I., Andriyash, E., Berkley, A., et al. Observation of topological phenomena in a programmable lattice of 1,800 qubits. Nature 560, 456{460. doi:10.1038/s41586-018-0410-x (2018). 21. Ottaviani, D. & Amendola, A. Low Rank Non-Negative Matrix Factorization with D-Wave 2000Q. arXiv:1808.08721 (2018). 22. Ohkuwa, M., Nishimori, H. & Lidar, D. A. Reverse annealing for the fully connected p-spin model. Phys. Rev. A 98, 022314. doi:10.1103/PhysRevA.98.022314 (2 2018). 23. Yamashiro, Y., Ohkuwa, M., Nishimori, H. & Lidar, D. A. Dynamics of reverse an- nealing for the fully-connectedp-spin model. arXiv preprint arXiv:1906.10889 (2019). 24. Venturelli, D. & Kondratyev, A. Reverse quantum annealing approach to portfolio optimization problems. Quantum Mach. Intell. 1, 17{30. doi:10.1007/s42484-019- 00001-w (2019). 25. Marshall, J., Venturelli, D., Hen, I. & Rieel, E. G. Power of pausing: Advancing understanding of thermalization in experimental quantum annealers. Physical Review Applied 11, 044083 (2019). 26. Crosson, E. J. & Lidar, D. A. Prospects for Quantum Enhancement with Diabatic Quantum Annealing. arXiv:2008.09913 (2020). 27. Yip, K. W., Albash, T. & Lidar, D. A. Quantum trajectories for time-dependent adiabatic master equations. Physical Review A 97, 022116 (2018). 28. Wiseman, H. & Milburn, G. Quantum Measurement and Control (Cambridge Univer- sity Press, 2010). 29. Breuer, H.-P. & Petruccione, F. The Theory of Open Quantum Systems (Oxford Uni- versity Press, 2002). 30. Carmichael, H. J. Statistical Methods in Quantum Optics 2: Non-Classical Fields (Springer Science & Business Media, 2009). 31. Dum, R., Parkins, A., Zoller, P. & Gardiner, C. Monte Carlo simulation of master equations in quantum optics for vacuum, thermal, and squeezed reservoirs. Phys. Rev. A 46, 4382 (1992). 32. Mlmer, K., Castin, Y. & Dalibard, J. Monte Carlo wave-function method in quantum optics. JOSA B 10, 524{538. doi:10.1364/josab.10.000524 (1993). 33. Daley, A. J. Quantum trajectories and open many-body quantum systems. Advances in Physics 63, 77{149. doi:10.1080/00018732.2014.933502 (2014). 34. Brun, T. A. A simple model of quantum trajectories. ajp 70, 719{737. doi:10.1119/ 1.1475328 (2002). 35. Imamoglu, A. Stochastic wave-function approach to non-Markovian systems. Physical Review A 50, 3650{3653 (1994). 162 36. Breuer, H.-P. Genuine quantum trajectories for non-Markovian processes. Physical Review A 70, 012106 (2004). 37. Caiaa, M., Smirne, A. & Bassi, A. Stochastic unraveling of positive quantum dynam- ics. Physical Review A 95, 062101{ (). 38. Davies, E. & Spohn, H. Open quantum systems with time-dependent Hamiltonians and their linear response. Journal of Statistical Physics 19, 511{523 (1978). 39. Heim, B., Rnnow, T. F., Isakov, S. V. & Troyer, M. Quantum versus classical an- nealing of Ising spin glasses. Science 348, 215{217 (2015). 40. Albash, T. & Lidar, D. A. Evidence for a Limited Quantum Speedup on a Quantum Annealer. arXiv:1705.07452 (2017). 41. Laine, E.-M., Piilo, J. & Breuer, H.-P. Measure for the non-Markovianity of quantum processes. Physical Review A 81, 062115 (2010). 42. Gardiner, C. & Zoller, P. Quantum Noise: A Handbook of Markovian and Non-Markovian Quantum Stochastic Methods with Applications to Quantum Optics (Springer, 2004). 43. Breuer, H.-P. & Piilo, J. Stochastic jump processes for non-Markovian quantum dy- namics. EPL (Europhysics Letters) 85, 50004 (2009). 44. Piilo, J., Maniscalco, S., H ark onen, K. & Suominen, K.-A. Non-Markovian quantum jumps. Physical review letters 100, 180402 (2008). 45. Piilo, J., H ark onen, K., Maniscalco, S. & Suominen, K.-A. Open system dynamics with non-Markovian quantum jumps. Physical Review A 79, 062112 (2009). 46. H ark onen, K. Jump probabilities in the non-Markovian quantum jump method. Jour- nal of Physics A: Mathematical and Theoretical 43, 065302 (2010). 47. Boixo, S., Smelyanskiy, V. N., Shabani, A., Isakov, S. V., Dykman, M., Denchev, V. S., et al. Computational multiqubit tunnelling in programmable quantum annealers. Nature communications 7, 1{7 (2016). 48. Boixo, S., Rnnow, T. F., Isakov, S. V., Wang, Z., Wecker, D., Lidar, D. A., et al. Evidence for quantum annealing with more than one hundred qubits. Nature physics 10, 218{224 (2014). 49. Gao, Y., Neuhauser, D., Baer, R. & Rabani, E. Sublinear scaling for time-dependent stochastic density functional theory. The Journal of chemical physics 142, 034106 (2015). 50. Paladino, E., Galperin, Y., Falci, G. & Altshuler, B. 1/f noise: Implications for solid- state quantum information. Reviews of Modern Physics 86, 361 (2014). 51. Bergli, J., Galperin, Y. M. & Altshuler, B. Decoherence in qubits due to low-frequency noise. New Journal of Physics 11, 025002 (2009). 52. Paladino, E., Faoro, L. & Falci, G. Decoherence due to discrete noise in josephson qubits. Advances in Solid State Physics, edited by B. Kramer,(Springer Berlin Heidel- berg, Berlin, Heidelberg, 2003), 747{762. 53. Faoro, L. & Viola, L. Dynamical suppression of 1/f noise processes in qubit systems. Physical review letters 92, 117905 (2004). 54. Radclie, J. M. Some properties of coherent spin states. Journal of Physics A: General Physics 4, 313{323. doi:10.1088/0305-4470/4/3/009 (1971). 55. Klauder, J. R. Path integrals and stationary-phase approximations. Physical Review D 19, 2349 (1979). 163 56. Bylander, J., Gustavsson, S., Yan, F., Yoshihara, F., Harrabi, K., Fitch, G., et al. Noise spectroscopy through dynamical decoupling with a superconducting ux qubit. Nature Physics 7, 565{570 (2011). 57. Quintana, C., Chen, Y., Sank, D., Petukhov, A., White, T., Kafri, D., et al. Observa- tion of classical-quantum crossover of 1/f ux noise and its paramagnetic temperature dependence. Physical review letters 118, 057702 (2017). 58. Gillespie, D. T. et al. Exact stochastic simulation of coupled chemical reactions. J. phys. Chem 81, 2340{2361. doi:10.1021/j100540a008 (1977). 59. Khezri, M., Grover, J. A., Basham, J. I., Disseler, S. M., Chen, H., Novikov, S., et al. Anneal-path correction in ux qubits. npj Quantum Information 7, 1{8 (2021). 60. Young, K. C., Sarovar, M. & Blume-Kohout, R. Error suppression and error correction in adiabatic quantum computation: Techniques and challenges. Physical Review X 3, 041013 (2013). 61. Lidar, D. A. Towards fault tolerant adiabatic quantum computation. Physical Review Letters 100, 160506 (2008). 62. Mishra, A., Albash, T. & Lidar, D. A. Performance of two dierent quantum annealing correction codes. Quantum Information Processing 15, 609{636 (2016). 63. Brun, T. A. A simple model of quantum trajectories. American Journal of Physics 70, 719{737 (2002). 64. Wiseman, H. M. & Milburn, G. J. Quantum measurement and control (Cambridge university press, 2009). 65. Passarelli, G., Yip, K.-W., Lidar, D. A., Nishimori, H. & Lucignano, P. Reverse quan- tum annealing of the p-spin model with relaxation. Physical Review A 101, 022331 (2020). 66. Jansen, S., Ruskai, M.-B. & Seiler, R. Bounds for the adiabatic approximation with applications to quantum computation. J. Math. Phys. 48, 102111 (2007). 67. Lidar, D. A., Rezakhani, A. T. & Hamma, A. Adiabatic approximation with exponen- tial accuracy for many-body systems and quantum computation. J. Math. Phys. 50, 102106 (2009). 68. Joye, A. Proof of the Landau{Zener formula. Asymptotic Analysis 9, 209 (1994). 69. Gross, D. & Mezard, M. The simplest spin glass. Nuclear Physics B 240, 431{452. doi:https://doi.org/10.1016/0550-3213(84)90237-2 (1984). 70. Derrida, B. Random-energy model: An exactly solvable model of disordered systems. Physical Review B 24, 2613 (1981). 71. Bapst, V. & Semerjian, G. On quantum mean-eld models and their quantum anneal- ing. Journal of Statistical Mechanics: Theory and Experiment 2012, P06007 (2012). 72. Childs, A. M., Farhi, E. & Preskill, J. Robustness of adiabatic quantum computation. Physical Review A 65, 012322 (2001). 73. Albash, T., Boixo, S., Lidar, D. A. & Zanardi, P. Quantum adiabatic Markovian master equations. New Journal of Physics 14, 123016 (2012). 74. Zanardi, P. & Rasetti, M. Noiseless quantum codes. Physical Review Letters 79, 3306 (1997). 75. Lidar, D. A., Chuang, I. L. & Whaley, K. B. Decoherence-free subspaces for quantum computation. Physical Review Letters 81, 2594 (1998). 164 76. Kempe, J., Bacon, D., Lidar, D. A. & Whaley, K. B. Theory of decoherence-free fault-tolerant universal quantum computation. Physical Review A 63, 042307 (2001). 77. D-Wave System Documentation. Technical Description of the D-Wave QPU. https: //docs.dwavesys.com/docs/latest/doc_qpu.html. Accessed: 2019-07-24. 78. Passarelli, G., De Filippis, G., Cataudella, V. & Lucignano, P. Dissipative environment may improve the quantum annealing performances of the ferromagnetic p-spin model. Physical Review A 97, 022319 (2018). 79. Venuti, L. C., Albash, T., Lidar, D. A. & Zanardi, P. Adiabaticity in open quantum systems. Physical Review A 93, 032118 (2016). 80. Venuti, L. C. & Lidar, D. A. Error reduction in quantum annealing using boundary cancellation: Only the end matters. Physical Review A 98, 022315 (2018). 81. Passarelli, G., Cataudella, V. & Lucignano, P. Improving quantum annealing of the ferromagnetic p-spin model through pausing. Physical Review B 100, 024302 (2019). 82. Verstraete, F., Wolf, M. M. & Cirac, J. I. Quantum computation and quantum-state engineering driven by dissipation. Nature physics 5, 633{636 (2009). 83. Venuti, L. C., Albash, T., Marvian, M., Lidar, D. & Zanardi, P. Relaxation versus adiabatic quantum steady-state preparation. Physical Review A 95, 042302 (2017). 84. Harris, R., Sato, Y., Berkley, A. J., Reis, M., Altomare, F., Amin, M. H., et al. Phase transitions in a programmable quantum spin glass simulator. Science 361, 162{165. doi:10.1126/science.aat2025 (2018). 85. King, A. D., Raymond, J., Lanting, T., Isakov, S. V., Mohseni, M., Poulin-Lamarre, G., et al. Scaling advantage in quantum simulation of geometrically frustrated magnets. arXiv:1911.03446 (2019). 86. Gardas, B., Dziarmaga, J., Zurek, W. H. & Zwolak, M. Defects in Quantum Comput- ers. Sci. Rep. 8, 4539. doi:10.1038/s41598-018-22763-2 (2018). 87. Nishimura, K., Nishimori, H. & Katzgraber, H. G. Griths-McCoy singularity on the diluted Chimera graph: Monte Carlo simulations and experiments on quantum hardware. Phys. Rev. A 102. doi:10.1103/physreva.102.042403 (2020). 88. King, A. D., Nisoli, C., Dahl, E. D., Poulin-Lamarre, G. & Lopez-Bezanilla, A. Quan- tum Articial Spin Ice. arXiv:2007.10555 (2020). 89. Kairys, P., King, A. D., Ozdan, I., Boothby, K., Raymond, J., Banerjee, A., et al. Sim- ulating the Shastry-Sutherland Ising Model using Quantum Annealing. arXiv:2003.01019 (2020). 90. Kirkpatrick, S., Gelatt, C. D. & Vecchi, M. P. Optimization by simulated annealing. science 220, 671{680 (1983). 91. Brady, L. T., Baldwin, C. L., Bapat, A., Kharkov, Y. & Gorshkov, A. V. Optimal protocols in quantum annealing and qaoa problems. arXiv preprint arXiv:2003.08952 (2020). 92. Mbeng, G. B., Fazio, R. & Santoro, G. Quantum Annealing: a journey through Digital- ization, Control, and hybrid Quantum Variational schemes. arXiv:1906.08948 (2019). 93. Gross, D. J. & Mezard, M. The simplest spin glass. Nucl. Phys. B 240, 431{452. doi:10.1016/0550-3213(84)90237-2 (1984). 94. Bapst, V. & Semerjian, G. On quantum mean-eld models and their quantum anneal- ing. J. Stat. Mech.: Theory Exp. 2012, P06007. doi:10.1088/1742-5468/2012/06/ p06007 (2012). 165 95. Shin, S. W., Smith, G., Smolin, J. A. & Vazirani, U. How" quantum" is the D-Wave machine? arXiv preprint arXiv:1401.7087 (2014). 96. Albash, T. & Marshall, J. Comparing relaxation mechanisms in quantum and classical transverse-eld annealing. Physical Review Applied 15, 014029 (2021). 97. Amin, M. H. Searching for quantum speedup in quasistatic quantum annealers. Phys- ical Review A 92, 052323 (2015). 98. D-Wave System Documentation { Spin-Bath Polarization Eect 2021. 99. Chen, H. & Lidar, D. A. Why and When Pausing is Benecial in Quantum Annealing. Physical Review Applied 14, 014100{. doi:10.1103/PhysRevApplied.14.014100 (2020). 100. Lidar, D. A. & Whaley, K. B. in Irreversible Quantum Dynamics 83 (Springer, Berlin, 2003). 101. Kubo, R. Statistical-mechanical theory of irreversible processes. I. General theory and simple applications to magnetic and conduction problems. J. Phys. Soc. Jpn. 12, 570{ 586. doi:10.1143/JPSJ.12.570 (1957). 102. Martin, P. C. & Schwinger, J. Theory of many-particle systems. I. Phys. Rev. 115, 1342. doi:10.1103/PhysRev.115.1342 (1959). 103. Albash, T., Vinci, W., Mishra, A., Warburton, P. A. & Lidar, D. A. Consistency tests of classical and quantum models for a quantum annealer. Phys. Rev. A 91, 042314. doi:10.1103/PhysRevA.91.042314 (4 2015). 104. Rnnow, T. F., Wang, Z., Job, J., Boixo, S., Isakov, S. V., Wecker, D., et al. Den- ing and detecting quantum speedup. Science 345, 420{424. doi:10.1126/science. 1252319 (2014). 105. Vestgaarden, J., Bergli, J. & Galperin, Y. Nonlinearly driven Landau-Zener transition in a qubit with telegraph noise. Physical Review B 77, 014514 (2008). 166 Appendices A Quantum trajectories for time-dependent adiabatic master equations A.1 Error estimates Let us assume that the system is in a pure state at time t, i.e., S (t) =j (t)ih (t)j, and let us consider a single time-step. In a single trajectory, the evolution ofj (t)i involves two possibilities: no jump or a jump. The ensemble average of trajectories after one nite time-step mainly involves two kinds of errors: the error associated with the norm of the state vector (Sec. A.1.1), and the error associated with the probability elements in a nite time-step (Sec. A.1.2). A.1.1 Error associated with the norm squared of a no-jump trajectory The Schr odinger equation of the eective Hamiltonian in the case of a no-jump trajectory is given by dj (t)i dt =iH e (t)j (t)i : (A.1) The resulting state after one time-step t is: j ~ (t + t)i =V e (t + t;t)j (t)i ; (A.2) 167 where V e (t + t;t) =T exp i Z t+t t H e (t 0 )dt 0 (A.3) is a non-unitary contractive evolution operator (kV e (t + t;t)k 1 for every induced norm) andj ~ (t + t)i is unnormalized since H e (t) is not hermitian. We can expand the time- ordered exponential for V e (t + t;t) as T exp i Z t+t t H e (t 0 )dt 0 (A.4) =1i Z t+t t dt 0 H e (t 0 ) + (i) 2 1 2! T Z t+t t dt 0 H e (t 0 ) 2 + (i) 3 1 3! T Z t+t t dt 0 H e (t 0 ) 3 ! For any H e (t) that is C K on an open interval containing [0;t f ] (H (k) e (t) is continuous and bounded for 1kK): Z t+t t dt 0 H e (t 0 ) = K1 X k=0 H (k) e (t) (k + 1)! (t) k+1 +O((t) K+1 ); (A.5) as t! 0. Assume H e (t) is C K and K 2. Let _ H e (t) H (1) e (t) denote the time derivative of H e (t). We have: Z t+t t dt 0 H e (t 0 ) =H e (t)t + 1 2! _ H e (t)(t) 2 +O((t) 3 ); (A.6a) T Z t+t t dt 0 H e (t 0 ) 2 = 2 Z t+t t dt 0 H e (t 0 ) Z t 0 t dt 00 H e (t 00 ) =H e (t) 2 (t) 2 + 1 2 H e (t) _ H e (t)(t) 3 + 1 2 _ H e (t)H e (t)(t) 3 +O((t) 4 ) ; (A.6b) T Z t+t t dt 0 H e (t 0 ) 3 =H e (t) 3 (t) 3 +O((t) 4 ): (A.6c) Therefore, we have, to second order in t: V e (t+t;t) =1i H e (t)t + 1 2! _ H e (t)(t) 2 +(i) 2 1 2! H e (t) 2 (t) 2 +O((t) 3 ); (A.7) 168 so that: V y e (t + t;t)V e (t + t;t) =1 +it H y e (t)H e (t) + 1 2! (t) 2 h i _ H y e (t)i _ H e (t) H y2 e (t) +H 2 e (t) i +O((t) 3 ) : (A.8) We approximate N(t; t)kj ~ (t + t)ik 2 by: N(t; t) =h ~ (t + t)j ~ (t + t)ih (t)j 1 +it H y e (t)H e (t) j (t)iE(t; t): (A.9) The approximation error is given by: (t; t)N(t; t)E(t; t) (A.10a) = 1 2 (t) 2 h (t)j h i _ H y e (t)i _ H e (t) H y2 e (t) +H 2 e (t) i j (t)i +O((t) 3 ) (A.10b) (t) 2 2 ki( _ H y e (t) _ H e (t)) (H y2 e (t) +H 2 e (t))k +kO((t) 3 )k ; (A.10c) wherekk is the operator norm (largest singular value). Eq. (A.10c) gives the relation between the error of the norm square approximation and the time-step. It is also helpful to consider the sources of error here as they will be used later. This error mainly arises from two sources during the truncations of Taylor expansion: 1. The truncation of the Taylor expansion of the integral [Eq. (A.6a)] to keep only H e (t)t. This turns Eq. (A.4) into exp (iH e (t)t). For this to hold, we require 1 2 _ H e (t)(t) 2 +O((t) 3 ) kH e (t)tk; (A.11) 169 implying: t 2 kH e (t)k k _ H e (t)k ; (A.12) assuming _ H e (t)6= 0; if it is then the condition becomes t K!kH e (t)k kH (K) e (t)k 1=K for the lowest value of K such that H (K) e (t)6= 0, where the superscript denotes the K-th derivative. 2. Keeping only the rst order term in Eq. (A.7) afterwards, i.e., 1 iH e (t). This requires, in addition to Eq. (A.12): 1 2 H 2 e (t)(t) 2 +O((t) 3 ) kH e (t)tk; (A.13) implying, for all t such that the denominators do not vanish: t min kH e (t)k k _ H e (t)k ; 1 kH e (t)k ; (A.14) where we ignored the factor of 2. In conclusion, the norm after one time-step is related to the jump probabilities as: N(t; t) =h (t)j 1 +i(t) H y e (t)H e j (t)i +(t; t) (A.15a) = 1 t X i hA y i (t)A i (t)i +(t; t) (A.15b) = 1 p(t) +(t; t); (A.15c) 170 where we used Eq. (2.10) and dened the (approximate) jump probabilities as: p(t) = X i p i (t) (A.16a) p i (t) = thA y i (t)A i (t)i : (A.16b) This explains Eqs. (2.19) and (2.20). A.1.2 Error associated with probability elements We have dened the p i (t) and p(t) in a xed time-step in Eq. (A.16). Note that even in the time-independent case, where both the Hamiltonian and Lindblad operators are time- independent, 1 p(t) is never exactly equal to the norm squared of the state vector after one time-step, and p i (t) is not exactly the jump probability inside the time-step. They are only approximations. For a nite time-step t, with p 0 the probability of having no jump inside the interval [t;t + t] andp 1 the probability of having one jump inside the interval [t;t + t], we have: p 0 =e R t+t t (t 0 )dt 0 = 1p 1 : (A.17) Note that as t! 0,p 0 = 1(t)t+o(t),p 1 '(t)t. We shall focus on the case where the time-step is suciently small such that the probability of two or more jumps occurring within a single time-step is negligible. 1 First, we can expand the exponential as p 0 = 1 Z t+t t (t 0 )dt 0 + 1 2 Z t+t t (t 0 )dt 0 2 1 E.g., the quantum jump can be described by a Poisson process with a state-dependent inhomogeneous jump rate, with two or more jumps as successive one and no-jump processes, and P n2 p n = o(t) as t! 0. 171 During the time window from t to t + t,j (t 0 )i [required to calculate (t 0 )] is obtained by solving the Schr odinger equation with the eective Hamiltonians [Eq. (A.1)] and renormal- izing the solution during the integration. As we show in Sec. A.3, in any nite interval the norm squared of the unnormalized state vector [Eq. (A.9)] is equal to p 0 . This is the reason why we can use the waiting time distribution (or Gillespie algorithm [58]) as our simulation method in the main text. Second, since Z t+t t dt 0 (t 0 ) = 1 X k=0 (k) (t) (k + 1)! (t) k+1 ; (A.18) we have = 1(t)t +e p ; (A.19) where the error e p associated with the probability elements in a xed time-step is: e p = 1 2 d dt (t)(t) 2 + 1 2 2 (t)(t) 2 +O((t) 3 ): (A.20) This should be much smaller than the rst order term (t)t. Therefore, we need t (t) 2 (t) d dt (t) : (A.21) In the time-independent case, this reduces to t 1 : (A.22) 172 A.2 Proof of equivalence between the master equation and trajectories formulations Our goal in this section is to show how the master equation, Eq. (2.9), can be recovered from the quantum trajectories formulation, and to nd a bound on the time-step t. This generalizes the proof for the time-independent case found in [32]. A.2.1 To jump or not to jump The probability elements p(t) and p i (t) are important for determining whether a jump occurs and if a jump does occur, which jump type occurs. In order to determine if a jump occurs or not, we draw a random number , uniformly distributed between 0 and 1. If p(t)<, which is almost always the case since p(t) is very small, no jump occurs. In the case of no jump,j (t)i evolves according to the eective Schr odinger equation, Eq. (A.1). At time t + t we simply renormalize the solution of Eq. (A.2): j (t + t)i = 1 q h ~ (t + t)j ~ (t + t)i j ~ (t + t)i (A.23a) (A.15c) = 1 p 1 p(t) + j ~ (t + t)i : (A.23b) If p(t) > , the state undergoes an abrupt jump and we choose the new wavefunction among the dierent states A i j (t)i and renormalize: j (t + t)i = A i (t)j (t)i q h (t)jA y i (t)A i (t)j (t)i (A.24a) (A.16b) = s t p i (t) A i (t)j (t)i : (A.24b) 173 Which type of jumps occurs is determined according to the probability i (t) = p i (t) p(t) = h (t)jA y i (t)A i (t)j (t)i t P i h (t)jA y i (t)A i (t)j (t)i t (A.25a) = hA y i (t)A i (t)i (t) ; (A.25b) where (t) = X i hA y i (t)A i (t)i (A.26) is the time-dependent jump rate. A.2.2 Averaging over trajectories LetH e (t) beC K withK 2. We rst express the mean value S (t) as a sum over the non- Hermitian evolution [with probability 1 p(t)] and the jump trajectories [with probability p(t)], so that as t! 0 we have: S (t + t) = [1 p(t)] j ~ (t + t)i p 1 p(t) + h ~ (t + t)j p 1 p(t) + + p(t) X i i (t) s t p i (t) A i (t)j (t)i s t p i (t) h (t)jA y i (t) (A.27a) (A.25b) = 1 p(t) 1 p(t) + j ~ (t + t)ih ~ (t + t)j + t X i A i (t) S (t)A y i (t) (A.27b) Combining Eq. (A.2) with Eq. (A.7) we have: j ~ (t + t)i = 1i H e (t)t + 1 2! _ H e (t)(t) 2 + (i) 2 1 2! H e (t) 2 (t) 2 +O((t) 3 ) j (t)i : (A.28) 174 Recall that in Eq. (A.14) we gave conditions allowing us to neglect the O((t) 2 ) terms. Thus: S (t + t) = 1 p(t) 1 p(t) + 1iH e (t)t +O((t) 2 ) S (t) 1 +iH y e (t)t +O((t) 2 ) (A.29) + t X i A i (t) S (t)A y i (t): where we have replaced S (t) by S (t) after averaging over many trajectories. Rearranging this expression into a form that exposes the terms that will become the master equation, the expression for the averaged state at t + t becomes: S (t + t) = S (t) +it S (t)H y e (t)H e (t) S (t) + t X i A i (t) S (t)A y i (t) (A.30a) 1 p(t) + h S (t) +i S (t) H y e (t)H e (t) S (t)t i +O((t) 2 ) : (A.30b) Note that is O((t) 2 ) as t! 0 [Eq. (A.10c)], and p(t) = t P i h (t)jA y i (t)A i (t)j (t)i is O(t) as t! 0, so that: 1 p(t) + = O((t) 2 ) 1O(t) +O((t) 2 ) (A.31a) =O((t) 2 ): (A.31b) Therefore line (A.30b) can be absorbed into O((t) 2 ), and we are left with: S (t + t) S (t) t =i H e (t) S (t) S (t)H y e (t) + X i A i (t) S (t)A y i (t) +O(t) ; (A.32) which becomes the master equation, Eq. (2.9), in the t! 0 limit. 175 A.2.3 Upper bound on t The above proof takes t! 0. We would like to know how small the time-step t should be in order for the approximations made to be valid. In Eq. (A.29), we expanded the time- ordered exponential, and kept only the rst order terms. This is equivalent to the criteria in Eqs. (A.11) and (A.13), summarized as a single condition in Eq. (A.14). As shown in Sec. A.1.1, this also automatically makes the error in the norm squared approximation small. We also need to satisfy Eq. (A.21), in order to accurately approximate the probability elements. Taken together, therefore: t min kH e (t)k k _ H e (t)k ; 1 kH e (t)k ; (t) 2 (t) _ (t) : (A.33) In practice, choosing a constant time-step that satises Eq. (A.33) in the whole timespan [0;t f ] is sucient, though one might prefer to implement an adaptive time-step tailored to the instantaneous value of the R.H.S. A.3 On the validity of waiting times (quantum time-dependent operators) Here we show the validity of using the waiting time distribution in the case of time-dependent operators. The argument presented here is based on Ref. [29] and we extend it to the time- dependent case. Let us denote byj (t)i andj ~ (t)i the normalized and unnormalized state vectors respec- tively, and let us assume they are equal at time t. This can happen when t = 0 or any time immediately after each jump. Lett + t +, where can be as large as is possible until the next jump occurs, and V e (t + ;t) =T exp " i Z t + t H e (t 0 )dt 0 # ; (A.34) 176 Then: j ~ (t +)i =V e (t + ;t)j (t)i ; (A.35a) j (t + )i = V e (t + ;t)j (t)i kV e (t + ;t)j (t)ik : (A.35b) Then, starting from t, for any future t + >t, we have d dt + j ~ (t + )i 2 = d dt + T exp " i Z t + t H e (t 0 )dt 0 # j (t)i 2 (A.36a) = d dt + h (t)jV y e (t + ;t)V e (t + ;t)j (t)i (A.36b) =h (t)jV y e (t + ;t)(+i)H y e (t + )V e (t + ;t)j (t)i +h (t)jV y e (t + ;t)(i)H e (t + )V e (t + ;t)j (t)i (A.36c) =h (t)jV y e (t + ;t) X i A y i (t + )A i (t + ) ! V e (t + ;t)j (t)i (A.36d) = V e (t + ;t)j (t)i 2 X i h (t + )jA y i (t + )A i (t + )j (t + )i ; (A.36e) where in the last equality we used Eq. (A.35b). Let N(t + )kj ~ (t + )ik 2 , as in Eq. (A.9). We have d dt + N(t + ) =N(t + )(t + ); (A.37) where is the time-dependent jump rate [Eq. (A.26)]. The solution to this dierential equation with the initial condition N(t) = 1 is N(t + ) = exp Z t + t (t 0 )dt 0 ! =p 0 (t + ); (A.38) wherep 0 is the probability of not having any jump inside the interval [Eq. (A.17)], which we have now shown to be equal to the the norm squared of the unnormalized state vector for any nite interval [t;t +]. 177 No commutators of operators at dierent times appears in the derivation. The use of the waiting time distribution is therefore valid for time-dependent operators as long as the correlation matrix is positive. A.4 Eect of the annealing schedule on the 8-qubit chain example We provide the functional form for the (D-Wave 1) annealing schedule functions A(t) and B(t) used in Sec. 2.5.1 in Fig. 8.1(a). We compare this non-linear schedule to a linear annealing schedule with an overall energy scale chosen to closely match the energy scale at which the rst annealing schedule curves cross. The dynamics associated with the linear annealing schedule [shown in Fig. 8.1(b)] is somewhat dierent than the non-linear one. This can be attributed to the dierences between the two schedules. In order to ensure that the two sets of schedules intersect at the same energy scale, this requires the linear schedule to start from a signicantly smaller energy scale. This results in some important eects on the dynamics. First, this energy scale is not large enough to ensure that the thermal state at s = 0 has negligible weight on excited states. Since we use the ground state as the initial state of the simulation, which is now suciently dierent from the thermal state at s = 0, the dissipative dynamics causes visible changes in the state immediately, which can be see as both the dip in Fig 8.1(a) near s = 0 and the large number of excitations at the beginning of the anneal in Fig. 8.1(c). Second, even after this initial dip, the lower energy scale associated with the linear schedule means that thermal depopulation of the ground state occurs generally sooner relative to the non-linear schedule studied in the main text, although at a lower rate because of the linear form of the schedule. Furthermore, because the transverse eld remains signicant compared to the Ising term for longer in the anneal than in the linear case, repopulation of the ground state due to thermal relaxation occurs for a longer period of time. However, ultimately, the ground state probability at the end of the anneal is not signicantly dierent than that shown in Fig. 2.3. The jump statistics are shown in Fig. 8.1(c), and closely resemble the non-linear case shown in the inset of Fig. 2.4. 178 0 0.2 0.4 0.6 0.8 1 0 5 10 15 20 25 30 35 Energy (GHz) (a) 0 0.2 0.4 0.6 0.8 1 0.4 0.5 0.6 0.7 0.8 0.9 1 Instantaneous GS Population QT ME 10 2 10 3 10 4 0.4 0.5 0.6 (b) (c) Figure 8.1: (a): Dotted lines: Linear schedule; Solid lines: D-Wave schedule used in Fig. 2.3; (b): Same as in Fig. 2.3 (the evolution of the population in the instantaneous ground state for the 8-qubit problem), but with a linear schedule. Inset: the convergence of the ground state population towards the AME results. (c): A histogram of the net number of jumps to the instantaneous ground state (GS). It shares the same pattern as in the non-linear schedule case (the inset of Fig. 2.4), but the net number of jumps out of the ground state is smaller due to the smaller energy scale in the linear schedule. A.5 Proof that eigenstates of H S are modied only under jumps Recall that the Lindblad operators are dened in Eq. (2.4) of the main text as: L ;! (t) = X a;b !;" b (t)"a(t) h" a (t)jA j" b (t)ij" a (t)ih" b (t)j: (A.39) After inverting Eq. (2.6) we have A i;! (t) = X u i; (!)L ;! (t) (A.40) for the new Lindblad operators corresponding to a diagonalized matrix of decay rates. Assume for simplicity that the matrix is already diagonal, so that u is just an identity transformation and i is a relabeling of ; then A i;! (t) = X a;b !;" b (t)"a(t) h" a (t)jA j" b (t)ij" a (t)ih" b (t)j A y i;! (t) = X a;b !;" b (t)"a(t) h" b (t)jA j" a (t)ij" b (t)ih" a (t)j: (A.41) 179 Consider the drift term 1 2 X i h A y i (t)A i (t)hA y i (t)A i (t)i i j (t)idt : Since A i (t) comes from the redenition where the index i includes the Bohr frequencies [Eq. (2.8)]: p 0 i (!)A i;! (t)!A i (t), the drift term becomes the following after we reintroduce the Bohr frequencies: 1 2 X i X ! 0 i (!) h A y i;! (t)A i;! (t)h (t)jA y i;! (t)A i;! (t)j (t)i i j (t)idt: (A.42) A.5.1 Whenj (t)i is an eigenstate of H S (t) If thej (t)i is an eigenstate of H S (t) (denoted asj" b (t)i here), A i;! (t)j (t)i =A i;! (t)j" b (t)i = X a !;" b (t)"a(t) h" a (t)jA j" b (t)ij" a (t)i h (t)jA y i;! (t) =h" b (t)jA y i;! (t) = X a !;" b (t)"a(t) h" b (t)jA j" a (t)ih" a (t)j: (A.43) Concentrate on the term of the parenthesis inside each summation of Eq. (A.42), i.e. h A y i;! (t)A i;! (t)h" b (t)jA y i;! (t)A i;! (t)j" b (t)i i j" b (t)i : (A.44) The rst term is: A y i;! (t)A i;! (t)j" b (t)i = X a 0 ;b 0 !;" b 0(t)" a 0(t) h" b 0(t)jA j" a 0(t)ij" b 0(t)ih" a 0(t)j ! X a !;" b (t)"a(t) h" a (t)jA j" b (t)ij" a (t)i = X a;b 0 !;" b 0(t)"a(t) !;" b (t)"a(t) h" b 0(t)jA j" a (t)ih" a (t)jA j" b (t)ij" b 0(t)i = X a;b 0 0;" b 0(t)" b (t) !;" b (t)"a(t) h" b 0(t)jA j" a (t)ih" a (t)jA j" b (t)ij" b 0(t)i; (A.45) 180 where the sum overb 0 denotes the sum overj" b 0i sharing the same energy asj" b i. The second term is: h" b (t)jA y i;! (t)A i;! (t)j" b (t)ij" b (t)i = X a 0 !;" b (t)" a 0(t) h" b (t)jA j" a 0(t)ih" a 0(t)j ! X a !;" b (t)"a(t) h" a (t)jA j" b (t)ij" a (t)i ! j" b (t)i = X a !;" b (t)"a(t) h" b (t)jA j" a (t)ih" a (t)jA j" b (t)i ! j" b (t)i : (A.46) Subtracting Eq. (A.46) from Eq. (A.45) yields the drift term [Eq. (A.44)], which is not zero, but a linear combination of degenerate eigenstates with the same energy " b (t). Before the jump happens the environment leads to the redistribution ofj" b (t)i to other states in the same energy manifold. (The Lamb shiftH LS (t) = P i;! S i (!)A y i;! (t)A i;! (t) also yields the same eect.) Since they all share the same energy, this does not aect the overlap with the ground state. If the evolution by H S (t) is adiabatic, such a linear combination will stay in the same energy manifold and this explains the square-pulse like behavior in the overlapping with the ground state in Fig. 2.4 of the main text. A.5.2 No degeneracy in " b (t) If there are no degenerate states with energy " b (t), A y i;! (t)A i;! (t)j" b (t)i = X a !;" b (t)"a(t) h" b (t)jA j" a (t)ih" a (t)jA j" b (t)ij" b (t)i (A.47) This cancels withh" b (t)jA y i;! (t)A i;! (t)j" b (t)ij" b (t)i [Eq. (A.46)] and the drift term [Eq. (A.44)] becomes zero. 181 A.6 Derivation of Eq. (2.31) When the state isj 1 (t)i, its jump rate 1!0 (t) toj 0 (t)i comprises the summation of Lindblad terms responsible for the 1! 0 transition: 1!0 (t) = X 2f1!0g hA y (t)A (t)i : (A.48) The summation is over the number of qubits n. Since each qubit is coupled to its own environment with an independent noise source, the matrix in Eq. (2.5) is already diagonal. From Eq. (2.4) we know that L ;! 10 (t) = X a;b ! 10 ;" b (t)"a(t) h" a (t)j z j" b (t)ij" a (t)ih" b (t)j: (A.49) Assume that Bohr frequency ! 10 (t) is due only to the 1! 0 transition (even if it is not, the other terms would be annihilated by the matrix elementh 1 (t)j:::j 1 (t)i). The Lindblad operators A (t) have the form: A (t) = p (! 10 )h 0 (t)j z j 1 (t)ij 0 (t)ih 1 (t)j: (A.50) Therefore: 1!0 (t) = X h 1 (t)jA y (t)A (t)j 1 (t)i (A.51) = n X =1 (! 10 )jh 0 (t)j z j 1 (t)ij 2 : Here (! 10 ) is evaluated with respect to the Ohmic spectral density. 182 B Experimental iterative reverse annealing B.1 Spectrum of the p-spin problem with p = 2 and scaling of the minimum gap with N The spectrum of p-spin problem and the ordering of energies of dierent spin sectors can be found in [94]. For the particular case of n = 4;p = 2 and the D-Wave 2000Q annealing schedule, we plot the spectrum in Fig. 8.2. From Fig. 8.2 we can see that the instantaneous rst excited state 1 (s) and the instantaneous ground state 0 (s) converge at s = 1. The energy gap between instantaneous eigenenergies i (s) and j (s) is ij (s) = i (s) j (s). The corresponding minimum energy gap is ij = min s ij (s). For p = 2, we are interested instead in the minimum energy gap between the ground state 0 (s) and the second excited state 2 (s), since the 1st excited state 1 (s) and the ground state converge at s = 1. It is denoted by . = 20 = min s 20 (s) = min s 2 (s) 0 (s): (B.1) We calculate and plot in Figure 6.5(c) the value of , for N =f4; ; 22g. We also show the s where the minimum gap is found, i.e. s min = argmin 20 (s): (B.2) B.2 Bound on success probabilities, for closed system and open- system dynamics with collective system-bath interaction If the dynamics preserve spin symmetry (for example in closed system Schr odinger equation and open-system simulation with the collective system-bath coupling assumption), there exists a natural upper bound on the maximum success probabilities achievable in reverse 183 0 0.2 0.4 0.6 0.8 1 -30 -20 -10 0 10 20 Figure 8.2: Full spectrum of 4 qubits. annealing. The upper bound is the population of the initial state in the maximum spin sector. For the example ofN = 4, there are two degenerate ground states (j0000i (j0ij"i,j1i j#i) andj1111i) of the problem Hamiltonian (H 0 ) and they both belong to the maximum spin subspace of S = S max = N 2 = 2. While the computational basis state with one spin down, for example,j0001i is a rst excited state ofH 0 , it is not a basis of the maximum spin subspace S = 2. A superposition of the computational basis states with one spin down, i.e. j0001i+j0010i+j0100i+j1000i p 4 , however belong to the maximum spin subspace S = 2, and this is the case for [65]. Unfortunately, the D-Wave machine does not allow such an initialization. In general, suppose that the initial computational basis state isj (t = 0)i comp , with a particular (normalized) magnetizationm 0 (Eq. 6.7). It is represented as a linear combination of states, each of which has a xed value of total spin momentum S and (normalized) magnetization m 0 , j (t = 0)i comp = Smax X S=S min a S jS;m S =m 0 i: (B.3) 184 From the theory of addition of angular momentum, we know that S min = N=2bN=2c and S max = N=2. The total spin momentum (integer or half-integer) S 2fS min ;S min + 1;:::;S max g. jS;m S =m 0 i is a simultaneous eigenstate of S 2 = P 2fx;y;zg S and S z = 1 2 P N i=1 z i (or H p ), with eigenvalues: S 2 jS;m S =m 0 i =S(S + 1)jS;m S =m 0 i ; (B.4) S z jS;m S =m 0 i = N 2 m 0 jS;m S =m 0 i (B.5) In Eq. (B.3),ja S j 2 is the initial state's population in that spin subspace. For the p-spin Hamiltonian, in a closed system the state in each spin subspace develops independently due to the preservation of spin symmetry of the Hamiltonian [23]. At the end of a single cycle of reverse annealing (r = 1), we have from Eq. (B.3) that U(2t inv ; 0)jS;m S =m 0 i = S X M=S c M jS;m S = 2M=Ni; (B.6) where c M is the relative phase developed in the other bases of the spin S subspace at the end of the time evolution. The (unnormalized) magnetization M = ( N 2 m S )2fS;S + 1;:::;Sg. Therefore, the nal state of 1-cycle reverse annealing (denoted byj i n ) can be expressed as: j i n =U(2t inv ; 0)j (t = 0)i comp = Smax X S=S min a S S X M=S c M jS;m S = 2M=Ni ! : (B.7) We know that the all-up state (j0i N =j"i N =jupi) and all-down state (j1i N =j#i N = jdowni) are both solution states ofH 0 , and moreover they lie in the maximum spin subspace 185 S = S max . Particularly,jupi =jS max ; 1i andjdowni =jS max ;1i. Therefore, projecting j i n onto thejupi gives: hupj i n = Smax X S=S min a S S X M=S c M hupjS;m S = 2M=Ni ! (B.8) =a Smax c Smax : (B.9) Similarly,hdownj i n =a Smax c Smax . The total success probability is thus bounded by: jhupj i n j 2 +jhdownj i n j 2 =ja Smax j 2 (jc Smax j 2 +jc Smax j 2 ) ja Smax j 2 ; (B.10) where equality is met with (jc Smax j 2 +jc Smax j 2 ) = 1. The upper boundja Smax j 2 is the population of the initial state in the maximum spin sub- space. For the initial state ofj0001i, we have a Smax = 1 p 4 sincejS =S max = 2;m S = 0:5i = j0001i+j0010i+j0100i+j1000i p 4 . Therefore, the total success probability is bounded byja Smax j 2 = 1 4 . For the other examples in the main text, the initial state ofj0011i has a Smax = 1 q ( 4 2 ) = 1 p 6 ; while the initial state ofj00000001i has a Smax = 1 p 8 . By expressing the above equations in terms of density matrix, we can also generalize the same conclusion to open-system simulation of collective system-bath coupling, where the population of each spin-sector is preserved during dynamics. 186 B.3 SVMC and SVMC-TF for reverse annealing We outline the algorithm of SVMC and SVMC-TF with application to reverse annealing. Algorithm 2 SVMC and SVMC-TF (p-spin reverse annealing) 1: H(s) = A(s) 2 ( P N i sin i ) B(s)N 2 1 N P N i cos i p , with a specied reverse annealing s =t dependence s(t). 2: procedure SVMC 3: for k = 1 to K (number of samples) do 4: Initial computational basis state:j0i i ! k i;t : 0,j1i i ! k i;t :. 5: for t = 1 to T (total number of sweeps) do 6: for i = 1 to n (number of qubits) do 7: Randomly choose a new angle 0 k i;t 2 [0;]. 8: Calculate E. 9: if E 0 then 10: k i;t ! 0 k i;t 11: else if p< exp(E) where p2 [0; 1] is drawn with uniform probability. then 12: k i;t ! 0 k i;t 13: end if 14: end for 15: end for 16: Take the mean of K samples: i;t = ( P K k=1 0 k i;t )=K. 17: end for 18: end procedure 19: return i;t . 20: procedure SVMC-TF 21: for k = 1 to K (number of samples) do 22: Initial computational basis state:j0i i ! k i;t : 0,j1i i ! k i;t :. 23: for t = 1 to T (total number of sweeps) do 24: for i = 1 to n (number of qubits) do 25: 0 k i;t = k i;t + i (s(t=T )), 26: a random i (s(t=T )2 [ min 1; A(s(t=T) B(s(t=T) ; min 1; A(s(t=T) B(s(t=T) ] 27: Calculate E. 28: if E 0 then 29: k i;t ! 0 k i;t 30: else if p< exp(E) where p2 [0; 1] is drawn with uniform probability. then 31: k i;t ! 0 k i;t 32: end if 33: end for 34: end for 35: Take the mean of K samples: i;t = ( P K k=1 0 k i;t )=K. 36: end for 37: end procedure 38: return i;t . 39: procedure Projection onto computational basis 40: if 0 k i;t =2 then 41: j k (t)i i =j0i 42: else if =2 k i;t then 43: j k (t)i i =j1i 44: end if 45: end procedure 46: Remark: E, for example, due to the update of qubit 2, can be expressed in the following form: E 2 = 0 B B @ A(s) 2 0 B B @ N X i=1 i6=2 sin i + sin 0 2 1 C C A B(s)N 2 0 B B @ 1 N 0 B B @ N X i=1 i6=2 cos i + cos 0 2 1 C C A 1 C C A 2 1 C C A 0 @ A(s) 2 N X i=1 sin i ! B(s)N 2 1 N N X i=1 cos i !! 2 1 A : 187
Abstract (if available)
Abstract
In this dissertation, we explore how quantum annealing (QA) and its applications behave in an open system setting. We give derivations and numerical recipes for effective parallel simulation methods for time-dependent open dynamics of quantum annealing devices. We consider the weak-coupling limit to an environment and also the case of 1/f noise. The stochastic time-dependent quantum trajectories technique can be utilized in studying weak measurements and feedback error correction in quantum annealing. Then we focus on open-system descriptions of reverse annealing (RA), which is a promising variant and application of quantum annealing. We show that, with various simulation tools, reverse annealing can benefit from the interaction between the annealer and its environment.
Linked assets
University of Southern California Dissertations and Theses
Conceptually similar
PDF
The theory and practice of benchmarking quantum annealers
PDF
Error correction and quantumness testing of quantum annealing devices
PDF
Theory and simulation of Hamiltonian open quantum systems
PDF
Topics in quantum information and the theory of open quantum systems
PDF
Tunneling, cascades, and semiclassical methods in analog quantum optimization
PDF
Open quantum systems and error correction
PDF
Minor embedding for experimental investigation of a quantum annealer
PDF
Topics in quantum information -- Continuous quantum measurements and quantum walks
PDF
Error suppression in quantum annealing
PDF
Advancing the state of the art in quantum many-body physics simulations: Permutation Matrix Representation Quantum Monte Carlo and its Applications
PDF
Dissipation as a resource for quantum information processing: harnessing the power of open quantum systems
PDF
Coherence generation, incompatibility, and symmetry in quantum processes
PDF
Quantum information-theoretic aspects of chaos, localization, and scrambling
PDF
Protecting Hamiltonian-based quantum computation using error suppression and error correction
PDF
Applications and error correction for adiabatic quantum optimization
PDF
Applications of quantum error-correcting codes to quantum information processing
PDF
Towards robust dynamical decoupling and high fidelity adiabatic quantum computation
PDF
Dynamic topology reconfiguration of Boltzmann machines on quantum annealers
PDF
Trainability, dynamics, and applications of quantum neural networks
PDF
Imposing classical symmetries on quantum operators with applications to optimization
Asset Metadata
Creator
Yip, Ka Wa
(author)
Core Title
Open-system modeling of quantum annealing: theory and applications
School
College of Letters, Arts and Sciences
Degree
Doctor of Philosophy
Degree Program
Physics
Publication Date
04/15/2021
Defense Date
03/12/2021
Publisher
University of Southern California
(original),
University of Southern California. Libraries
(digital)
Tag
1/f noise,OAI-PMH Harvest,quantum annealing,quantum feedback,quantum open system,reverse annealing,weak measurement
Language
English
Contributor
Electronically uploaded by the author
(provenance)
Advisor
Lidar, Daniel (
committee chair
), Brun, Todd (
committee member
), Haas, Stephan (
committee member
), Spedalieri, Federico (
committee member
), Zanardi, Paolo (
committee member
)
Creator Email
kawayip@usc.edu
Permanent Link (DOI)
https://doi.org/10.25549/usctheses-c89-443939
Unique identifier
UC11668949
Identifier
etd-YipKaWa-9471.pdf (filename),usctheses-c89-443939 (legacy record id)
Legacy Identifier
etd-YipKaWa-9471.pdf
Dmrecord
443939
Document Type
Dissertation
Rights
Yip, Ka Wa
Type
texts
Source
University of Southern California
(contributing entity),
University of Southern California Dissertations and Theses
(collection)
Access Conditions
The author retains rights to his/her dissertation, thesis or other graduate work according to U.S. copyright law. Electronic access is being provided by the USC Libraries in agreement with the a...
Repository Name
University of Southern California Digital Library
Repository Location
USC Digital Library, University of Southern California, University Park Campus MC 2810, 3434 South Grand Avenue, 2nd Floor, Los Angeles, California 90089-2810, USA
Tags
1/f noise
quantum annealing
quantum feedback
quantum open system
reverse annealing
weak measurement