Real-time simulation of flux qubits used for quantum annealing

TL;DR

This paper uses numerical simulations to analyze flux qubit dynamics during quantum annealing, validating the qubit model and environmental effects by comparing with experimental data from a D-Wave system.

Contribution

It demonstrates that flux qubits can be accurately modeled during annealing and highlights the importance of environmental coupling for matching experimental results.

Findings

Qubit representation effectively describes flux dynamics.

Tunable coupling does not impair annealing performance.

Environmental effects improve simulation-experiment agreement.

Abstract

The real-time flux dynamics of up to three superconducting quantum interference devices (SQUIDs) are studied by numerically solving the time-dependent Schr\"odinger equation. The numerical results are used to scrutinize the mapping of the flux degrees of freedom onto two-level systems (the qubits) as well as the performance of the intermediate SQUID as a tunable coupling element. It is shown that the qubit representation yields a good description of the flux dynamics during quantum annealing and the presence of the tunable coupling element does not have negative effects on the overall performance. Additionally, data obtained from a simulation of the dynamics of two-level systems during quantum annealing are compared to experimental data produced by the D-Wave 2000Q quantum annealer. The effects of finite temperature are incorporated in the simulation by coupling the qubit system to a bath of two-level systems. It is shown that an environment modeled as non-interacting two-level systems coupled to the qubits can produce data which matches the experimental data much better than the simulation data of the qubits without coupling to an environment and better than data obtained from a simulation of an environment modeled as interacting two-level systems coupling to the qubits.

Figures

Tables

Table 1. Table 1: Values of the device parameters appearing in the Hamiltonian Eq. ( 8 ) and used in our numerical work.

Parameter	Value
$E_C$	$4.68\mathrm{GHz}$
$E_L$	$3.48\mathrm{THz}$
$E_{C_J}$	$133.02\mathrm{GHz}$
$E_{L_J}$	$54.54\mathrm{THz}$
$E_J$	$7.80\mathrm{THz}$
$M$	$15.97\mathrm{pH}$

Table 2. Table 2: Percentage for finding the ground state (GS) on D-Wave’s DW_2000Q_2 chip for three different problem instances and four different annealing times. The minimal and final gap (in GHz) are denoted by δ ​ E 𝛿 𝐸 \delta E and Δ p subscript Δ 𝑝 \Delta_{p} , respectively.

Case	$J$	$h_1$	$h_2$	$\delta E$	${\Delta }_p$	GS probability in %
						$1\mu$s	$20\mu$s	$100\mu$s	$1\mathrm{ms}$
(a)	$-1$	$0$	$0.05$	$1.206$	$3.519$	$63.0$	$65.6$	$67.1$	$69.7$
(b)	$-1$	$0.96$	$0.94$	$0.627$	$1.407$	$51.4$	$52.9$	$53.6$	$55.6$
(c)	$0.1$	$0.3$	$-0.3$	$5.481$	$14.07$	$92.5$	$96.2$	$97.6$	$98.5$

Table 3. Table 3: Parameter values of the problems shown in Fig.  8 . (a) Success probability for the qubit model ( t a = 5 subscript 𝑡 𝑎 5 t_{a}=5 ns). (b) Success probability for the full model ( t a = 5 subscript 𝑡 𝑎 5 t_{a}=5 ns). (c) Percentage of successful runs on D-Wave’s DW_2000Q_2 and DW_2000Q_2_1 chips ( t a = 20 ​ μ subscript 𝑡 𝑎 20 𝜇 t_{a}=20\mu s).

$h_1$	$h_2$	$J$	Minimal gap	Success probability in %
			$\delta E$ in GHz	(a)	(b)	(c)
$0.2$	$0.2$	$0.2$	$7.958914$	$99.9$	$99.5$	$99.7$
$0.2$	$-0.2$	$0$	$6.524809$	$99.8$	$99.6$	$95.8$
$0.3$	$-0.3$	$0.1$	$6.509859$	$99.9$	$99.6$	$96.2$
$0.1$	$-0.1$	$-0.1$	$4.817172$	$96.5$	$96.3$	$94.0$
$0.9$	$0.7$	$-1$	$4.660374$	$96.3$	$96.2$	$93.9$
$0.99$	$-0.8$	$1$	$4.367788$	0$95.8$0	0$95.2$0	0$94.1$0
$0.1$	$0.1$	$0$	$3.750846$	$93.0$	$92.8$	$80.5$
$0.3$	$0.3$	$-0.2$	$3.740396$	$96.2$	$96.2$	$81.4$
$0.07$	$-0.07$	$0$	$2.786031$	$82.1$	$81.6$	$66.8$
$0.07$	$0.07$	$0$	$2.786031$	$82.1$	$82.1$	$66.9$
$0.9$	$-0.8$	$1$	$2.766581$	$83.8$	$83.9$	$80.0$
$0.02$	$0.08$	$0.05$	$2.542547$	$77.7$	$77.6$	$71.2$
$0.99$	$-1$	$0.94$	$2.134413$	$84.8$	$57.9$	$48.3$
$0.05$	$0.05$	$0$	$2.092326$	$69.3$	$69.4$	$56.8$
$0$	$-0.05$	$0.05$	$1.585987$	$60.3$	$60.3$	$55.6$
$0$	$0.05$	$1$	$1.433807$	$69.0$	$67.6$	$64.4$
$0$	$0.05$	$-1$	$1.433807$	$69.0$	$69.1$	$65.3$
$0.01$	$0.04$	$0.025$	$1.419405$	$55.6$	$55.7$	$50.4$
$0.99$	$-1$	$0.96$	$1.366784$	$74.7$	$42.3$	$38.1$
$0.99$	$1$	$-0.96$	$1.366784$	$74.7$	$74.4$	$46.9$
$0.02$	$-0.02$	$-0.02$	$1.305954$	$51.1$	$50.6$	$43.8$
$0.02$	$0.02$	$0.02$	$1.305954$	$51.1$	$51.2$	$43.9$
$0.95$	$-0.99$	$0.98$	$1.145772$	$47.5$	$63.4$	$55.7$
$0.95$	$0.99$	$-0.98$	$1.145772$	$47.5$	$47.3$	$52.3$
$0.99$	$0.96$	$-1$	$1.018001$	$48.4$	$48.0$	$54.8$
$0.02$	$-0.02$	$0$	$0.939407$	$42.7$	$42.2$	$36.2$
$0$	$0.03$	$1$	$0.917871$	$61.8$	$60.8$	$59.0$
$0.96$	$-0.94$	$1$	$0.742309$	$50.7$	$57.5$	$55.1$
$0.96$	$0.94$	$-1$	$0.742309$	$50.7$	$50.0$	$54.0$
$0.98$	$-0.96$	$1$	$0.740480$	$45.4$	$56.3$	$53.4$
$0.98$	$0.96$	$-1$	$0.740480$	$45.4$	$44.7$	$51.0$
$0.01$	$-0.01$	$-0.01$	$0.716800$	$37.8$	$37.3$	$33.2$
$0.01$	$0.01$	$0.01$	$0.716800$	$37.8$	$38.0$	$33.9$
$0$	$0.02$	$-1$	$0.640932$	$57.9$	$57.8$	$57.0$
$0.01$	$0.009$	$0.002$	$0.543047$	$33.9$	$34.2$	$31.4$
$0.99$	$-1$	$0.98$	$0.510495$	$61.1$	$27.0$	$27.5$
$0.01$	$0.01$	$0$	$0.505201$	$33.5$	$33.8$	$31.7$
$0.99$	$-0.98$	$1$	$0.406202$	$35.3$	$49.3$	$47.2$
$0.99$	$0.98$	$-1$	$0.406202$	$35.3$	$34.6$	$45.6$
$0$	$-0.01$	$0.01$	$0.392828$	$33.0$	$33.2$	$31.6$
$0.005$	$0.005$	$0.005$	$0.388173$	$31.2$	$31.5$	$29.8$
$0.005$	$-0.005$	$-0.005$	$0.388173$	$31.2$	$30.7$	$29.1$
$0$	$0.01$	$1$	$0.343987$	$54.0$	$53.4$	$52.7$
$0$	$0.01$	$-1$	$0.343987$	$54.0$	$53.8$	$53.1$
$0.007$	$0$	$-0.01$	$0.302380$	$31.8$	$31.3$	$30.4$
$0.007$	$0$	$0.01$	$0.302380$	$31.8$	$32.0$	$30.5$
$1$	$0$	$-0.005$	$0.269309$	$54.0$	$53.1$	$52.8$
$1$	$0$	$0.005$	$0.269309$	$54.0$	$54.6$	$52.8$
$0.005$	$0.001$	$0.01$	$0.268263$	$31.4$	$31.6$	$30.2$
$0.005$	$-0.001$	$-0.01$	$0.268263$	$31.4$	$30.9$	$30.0$
$0.005$	$0$	$0.01$	$0.228548$	$30.9$	$31.2$	$29.6$
$0$	$0.005$	$0.5$	$0.193353$	$52.0$	$51.8$	$51.3$
$0$	$0.005$	$-0.5$	$0.193353$	$52.0$	$51.9$	$52.3$
$0.005$	$-0.001$	$0.01$	$0.187420$	$30.5$	$30.7$	$29.3$
$0.005$	$0.001$	$-0.01$	$0.187420$	$30.5$	$30.0$	$28.4$
$0$	$0.005$	$1$	$0.183074$	$52.0$	$51.5$	$50.7$
$0.005$	$0$	$-1$	$0.183074$	$52.0$	$51.8$	$50.7$
$0.003$	$0$	$0.01$	$0.145233$	$30.1$	$30.3$	$29.2$
$0$	$0.003$	$-1$	$0.114519$	$51.2$	$51.0$	$51.0$
$0$	$0.003$	$1$	$0.114519$	$51.2$	$50.8$	$50.1$

Equations

$$H(s)=A(s)H_{\mathrm{init}}+B(s)H_{\mathrm{final}},$$

$$\underset{s\in [0,1]}{\max }\frac{|⟨E{S}_j(s)|{\partial }_{s}H(s)|GS(s)⟩|}{\Delta E_j(s{)}^2}\ll t_a.$$

$$H_{\mathrm{QUBO}}=-\sum _{k=1}^N{h}_k{S}_k-\sum _{1\le j<k}{J}_{jk}{S}_j{S}_k,$$

$$H_{\mathrm{Ising}}=-\sum _{k=1}^N{h}_k{\sigma }_k^z-\sum _{1\le j<k}{J}_{jk}{\sigma }_j^z{\sigma }_k^z.$$

$$H_{\mathrm{trans}}=-\sum _{k=1}^N{\sigma }_k^x,$$

$$H_{\mathrm{SQUID}}=$$

$$-E_{C_{\mathrm{J}}}{\partial }_{{\phi }_{\mathrm{J}}}^2+E_{L_{\mathrm{J}}}\frac{({\phi }_{\mathrm{J}}-{\phi }_{\mathrm{J}}^x{)}^2}{2}$$

$$-E_J\cos (\phi )\cos \left(\frac{{\phi }_{\mathrm{J}}}{2}\right),$$

$$H_{\mathrm{int}}$$

$$+\frac{M}{L_{\mathrm{eff}}}{E}_L\left({\phi }_2-{\phi }_2^x\right)\left({\phi }_0-{\phi }_0^x\right)$$

$$+\frac{M^2}{L{L}_{\mathrm{eff}}}{E}_L\left({\phi }_1-{\phi }_1^x\right)\left({\phi }_2-{\phi }_2^x\right),$$

$$H_{\mathrm{total}}$$

$$+E_{L_{\mathrm{eff}}}\frac{{\phi }_0^2}{2}-E_{C_0}{\partial }_{{\phi }_0}^2+E_{L_{\mathrm{J},0}}\frac{({\phi }_{\mathrm{J},0}-{\phi }_{\mathrm{J},0}^x{)}^2}{2}-E_{C_{\mathrm{J},0}}{\partial }_{{\phi }_{\mathrm{J},0}}^2-E_{J_0}\cos ({\phi }_0)\cos \left(\frac{{\phi }_{\mathrm{J},0}}{2}\right)$$

$$+\frac{M}{L_{\mathrm{eff}}}{E}_L\left({\phi }_1-{\phi }_1^x\right){\phi }_0+\frac{M}{L_{\mathrm{eff}}}{E}_L\left({\phi }_2-{\phi }_2^x\right){\phi }_0+\frac{M^2}{L{L}_{\mathrm{eff}}}{E}_L\left({\phi }_1-{\phi }_1^x\right)\left({\phi }_2-{\phi }_2^x\right).$$

$$H_{\mathrm{int}}^{\mathrm{eff}}=C({\phi }_{\mathrm{J},0}^x){\phi }_1{\phi }_2,$$

$$\displaystyle H_{\mathrm{co}}=\frac{E_{L_{\mathrm{J},0}}^{\prime}}{2}\Bigg{(}\varphi_{\mathrm{J},0}$$

$$-E_{C_{\mathrm{J},0}}{\partial }_{{\phi }_{\mathrm{J},0}}^2,$$

$$\displaystyle\langle\varphi_{\mathrm{J},0}\rangle=\varphi_{\mathrm{J},0}^{x}-\frac{2E_{J_{0}}\sin\left(\varphi_{\mathrm{J},0}^{x}/2\right)}{4E_{L_{\mathrm{J},0}}+E_{J_{0}}\cos\big{(}\varphi_{\mathrm{J},0}^{x}/2\big{)}}.$$

$$H^{\mathrm{eff}}$$

$$T$$

$${\beta }_{\mathrm{eff}}=\frac{E_{J_0}}{E_{L_{\mathrm{eff}}}}\cos \left(\frac{⟨{\phi }_{\mathrm{J},0}⟩}{2}\right),$$

$$\cos \left({\phi }_0-\frac{M}{L(1+{\beta }_{\mathrm{eff}})}\left({\phi }_1-{\phi }_1^x+{\phi }_2-{\phi }_2^x\right)\right)$$

$$\displaystyle H^{\mathrm{eff}}=\sum\limits_{i=1}^{2}\Bigg{[}$$

$$\displaystyle+\sum\limits_{i=1}^{2}\Bigg{[}$$

$$\displaystyle-\Big{(}E_{C_{\mathrm{0}}}$$

$$H(s)=-A(s)({\sigma }_1^x+{\sigma }_2^x)-B(s)(h_1{\sigma }_1^z+h_2{\sigma }_2^z+J{\sigma }_1^z{\sigma }_2^z).$$

$$|\uparrow ⟩$$

$$|\downarrow ⟩$$

$${\tilde{E}}_L=E_L\left(1+\frac{M^2}{L{L}_{\mathrm{eff}}}\frac{{\beta }_{\mathrm{eff}}}{1+{\beta }_{\mathrm{eff}}}\right)$$

$${\phi }_i^x=h_i\gamma \frac{2e{I}_p(s)M^2}{L_{\mathrm{eff}}},$$

$$H^{\mathrm{eff},q}\approx$$

$$-\frac{E_L}{{\tilde{E}}_L}\frac{I_p^2(s)M^4}{L{L}_{\mathrm{eff}}^2}\frac{{\beta }_{\mathrm{eff}}}{1+{\beta }_{\mathrm{eff}}}\gamma \left(h_1{\sigma }_2^z+h_2{\sigma }_1^z\right)$$

$$+\frac{E_L^2}{{\tilde{E}}_L^2}\frac{I_p^2(s)M^2}{L_{\mathrm{eff}}}\frac{{\beta }_{\mathrm{eff}}}{1+{\beta }_{\mathrm{eff}}}{\sigma }_1^z{\sigma }_2^z.$$

$$-\gamma =-\underset{{\phi }_{\mathrm{J},0}^x}{\max }\frac{{\beta }_{\mathrm{eff}}}{1+{\beta }_{\mathrm{eff}}}\frac{E_L^2}{{\tilde{E}}_L^2}\le -J\gamma \le \underset{{\phi }_{\mathrm{J},0}^x}{\max }\frac{{\beta }_{\mathrm{eff}}}{1+{\beta }_{\mathrm{eff}}}\frac{E_L^2}{{\tilde{E}}_L^2}=\gamma .$$

$$\frac{{\beta }_{\mathrm{eff}}({\phi }_{\mathrm{J},0}^x)}{1+{\beta }_{\mathrm{eff}}({\phi }_{\mathrm{J},0}^x)}\frac{E_L^2}{{\tilde{E}}_L^2({\phi }_{\mathrm{J},0}^x)}=-J\gamma$$

$$H^{\mathrm{eff},q}\approx$$

Full text

quantum computation; superconducting qubits; quantum annealing; SQUID; product-formula algorithm; optimization problems

I Introduction

Theoretically, an ideal quantum computer is described in terms of two-level systems, the qubits Nielsen and Chuang (2010). However, almost all currently popular technologies such as ion traps Cirac and Zoller (1995); Monroe et al. (1995); Häffner et al. (2003); Hanneke et al. (2009); Schindler et al. (2013); Ballance et al. (2016), quantum dots Loss and DiVincenzo (1998); Levy (2002), or superconducting circuits Wendin (2017); Krantz et al. (2019) employ physical devices which are only approximately described by two-level systems Wu and Lidar (2002). Among these, trapped ions may be described in the two-level approximation under conditions discussed in Ref. Leibfried et al. (2003). Single-electron quantum dots can be described as two-level systems if the orbital wave function is neglected. However, two-electron Levy (2002); Petta et al. (2005) and three-electron systems DiVincenzo et al. (2000); Gaudreau et al. (2011) confined in quantum dots are again only approximately described by two-level systems, and leakage out of the computational space may need to be taken into account DiVincenzo et al. (2000); Medford et al. (2013); Cerfontaine et al. (2016).

For superconducting circuits, it depends on the particular circuit design how well the two relevant energy levels, which define the qubit states, are separated from the higher energy levels. For instance, the phase qubit Martinis et al. (2002); Steffen et al. (2006) and the transmon Koch et al. (2007), an extension of the charge qubit Shnirman et al. (1997); Bouchiat et al. (1998); Nakamura et al. (1999), have rather small anharmonicities. Thus, leakage to higher energy levels is a major issue when performing gate operations Motzoi et al. (2009); Lucero et al. (2010); Gambetta et al. (2011); Wallman et al. (2016); Willsch et al. (2017); Wood and Gambetta (2018) and is typically alleviated by the use of pulse-shaping techniques Motzoi et al. (2009); Chow et al. (2010); Lucero et al. (2010); Gambetta et al. (2011); Chen et al. (2016); McKay et al. (2017). Flux qubits, on the other hand, have a strong anharmonicity and are less prone to excitations to higher energy levels Poletto et al. (2009); Yoshihara et al. (2014); Billangeon et al. (2015); Chen et al. (2016); Krantz et al. (2019) as long as the qubit is not driven too strongly Ferrón and Domínguez (2010).

In this study we focus on the superconducting quantum interference device (SQUID)-based flux qubit Chiarello (2000); Makhlin et al. (2001); Poletto et al. (2009), which is used in the D-Wave quantum annealer Harris et al. (2009a, 2010a). Due to the large superconducting loop needed for the SQUID-based qubit, it is sensitive to flux noise which limits the coherence time Wendin (2017).

Other flux qubits, commonly used for gate-based quantum computing, are the three-junction qubit Mooij et al. (1999); Orlando et al. (1999); van der Wal et al. (2000); Grajcar et al. (2004), the C-shunt flux qubit You et al. (2007); Steffen et al. (2010); Yan et al. (2016), and the fluxonium qubit Manucharyan et al. (2009); Pop et al. (2014); Nguyen et al. (2019). For the three-junction qubit and the fluxonium qubit, the large inductance is realized by using two or more Josephson junctions. In this way, the loop size and thus the sensitivity to flux noise can be reduced Wendin and Shumeiko (2007). The C-shunt flux qubit is a capacitively shunted variant of the flux qubit with improved coherence when operated away from the degeneracy point You et al. (2007). There is an ongoing discussion about the role of decoherence during quantum annealing Childs et al. (2001); Sarandy and Lidar (2005); Ashhab et al. (2006); Amin et al. (2009a); Boixo et al. (2013); Dickson et al. (2013).

In this paper we address three questions. First, we study the flux dynamics of the SQUIDs used in the D-Wave quantum annealer, addressing the issue of how well these dynamics are captured by a qubit model.

Second, we investigate whether the presence of the SQUID used as a tunable coupler in the D-Wave device affects the performance of the quantum annealing process. To answer the first two questions, we study the dynamics of the flux degrees of freedom of two SQUIDs functioning as qubits and their tunable coupler, a third SQUID, by solving the time-dependent Schrödinger equation (TDSE) for the model Hamiltonian in terms of flux degrees of freedom. So far, studies including higher energy levels have been limited to four-level qudits Johnson et al. (2011); Amin et al. (2013).

The approach that we adopt in this paper is to start from an idealized model of the SQUIDs which does not take into account fabrication variations of circuit elements or stray fluxes induced by the control lines. That is, in the idealized model, the two SQUIDs functioning as qubits are equal and the complete system can be regarded as a perfect device.

The third question we consider is to what extent the data produced by a D-Wave quantum annealer can be described by quantum annealing of the qubit model including environment effects. In order to do so, we study the dynamics of a two-qubit system interacting with an environment of two-level systems, representing, e.g., a heat bath Zhao et al. (2016); De Raedt et al. (2017) or a collection of defects described by non-interacting two-level systems Shnirman et al. (2005); Müller et al. (2009); Cole et al. (2010), by solving the corresponding TDSE.

The structure of the paper is as follows. In Sec. II we give a short introduction to quantum annealing and its relation to optimization problems. Section III contains a description of the SQUID-based model which we simulate. The mapping of the model onto the qubit model is given in Sec. IV. The flux dynamics of the SQUID model are simulated by solving the TDSE using the method described in Sec. V. The results of the simulation are presented in Sec. VI. In Sec. VII we describe the two different models for the bath of two-level systems coupled to the qubit system and discuss the simulation and its results in comparison to data obtained from a D-Wave 2000Q quantum annealer. We conclude with a summary in Sec. VIII.

II Theoretical background

In general, the Hamiltonian describing a quantum annealing process can be written as

[TABLE]

where $H_{\mathrm{init}}$ is the initial Hamiltonian whose ground state defines the state in which the system is prepared initially, $H_{\mathrm{final}}$ denotes the Hamiltonian at the end of the annealing process and whose ground state is the one to be determined, $s=t/t_{a}\in[0,1]$ is the rescaled (dimensionless) time, and $t_{a}$ is the total annealing time. The functions $A(s)$ and $B(s)$ determine the energy scale (in our case GHz) and the annealing scheme. They satisfy $|A(0)|\gtrsim 1$ and $B(0)\approx 0$, and $A(1)\approx 0$ and $|B(1)|\gtrsim 1$ with respect to the corresponding energy scale.

From the adiabatic theorem Born and Fock (1928), it follows that the system stays in the ground state of the instantaneous Hamiltonian $H(s)$ during the annealing process if $t_{a}\rightarrow\infty$ such that for $s=1$ the system is in the ground state of $H_{\mathrm{final}}$. Let $\Delta E_{j}(s)$ denote the difference between the energy of the ground state $\ket{GS(s)}$ and the $j$-th excited state $\ket{ES_{j}(s)}$ at the rescaled time $s$. A finite $t_{a}$ can be sufficient for the system to stay in the ground state if Amin (2009); Albash and Lidar (2018)

[TABLE]

Quantum annealing can be used to solve optimization problems that can be mapped onto the Hamiltonian $H_{\mathrm{final}}$. The class of so-called quadratic unconstrained binary optimization (QUBO) problems can be mapped onto the Ising-spin model of the form

[TABLE]

where $N$ is the number of binary variables $S_{k}\in\{-1,1\}$, and $h_{k}$ and $J_{jk}$ are dimensionless real numbers defining the particular QUBO. The set of variables $\{S_{k}\}$ that minimizes Eq. (3) gives the solution of the QUBO problem.

Quantum annealing can, at least in principle, find (one of) the ground state(s) of the Ising-spin Hamiltonian Eq. (3) Kadowaki and Nishimori (1998) or, equivalently, solve the corresponding QUBO problem. For this purpose, the two-valued variables $S_{k}$ are replaced by the Pauli $Z$ matrices $\sigma_{k}^{z}$ with eigenvalues $\pm 1$ and eigenstates $\ket{\uparrow}$ and $\ket{\downarrow}$ such that Eq. (3) transforms into

[TABLE]

The product states of the $\sigma^{z}$ eigenstates define the so-called computational basis and are eigenstates of Eq. (4). The ground state of the Hamiltonian Eq. (4) can then be found by quantum annealing with $H_{\mathrm{final}}$ replaced by $H_{\mathrm{Ising}}$ in Eq. (1). For quantum annealing, the simplest choice for $H_{\mathrm{init}}$ is the Hamiltonian of spins in a transverse field Kadowaki and Nishimori (1998)

[TABLE]

where $\sigma^{x}$ is the Pauli $X$ matrix. The ground state of this Hamiltonian is given by the product state $\ket{+\dots+}$, with $\ket{+}=(\ket{\uparrow}+\ket{\downarrow})/\sqrt{2}$.

Equation (4) is used to formulate optimization problems for the quantum annealer manufactured by D-Wave Systems Inc. Harris et al. (2009a). By design, the parameters of the final Hamiltonian are restricted to $h_{k}\in[-2,2]$ and $J_{jk}\in[-1,1]$, and the qubit connectivity is given by the Chimera graph such that, in the notation of Eq. (4), some $J_{jk}$ have to be set to zero Harris et al. (2010b). In the following sections, we discuss the SQUID model and describe the mapping of the SQUID model onto the qubit model in terms of Eqs. (4) and (5).

III SQUID Model

In this section we introduce the model Hamiltonian of the three-SQUID system that is used to simulate the flux dynamics during quantum annealing. Two of the three SQUIDs serve as qubits, each qubit subspace being defined by projection onto the two lowest energy states of the individual SQUIDs. The third SQUID acts as a tunable coupler between the two other SQUIDs.

Figure 1 shows the circuit of a SQUID with a compound Josephson junction (CJJ) loop. It is used as a building block for the flux qubits and the effective coupling between them in the D-Wave quantum annealer. Including the CJJ loop effectively leads to a tunable Josephson junction Chiarello (2000). The two qubit states correspond to the left-circulating and right-circulating persistent current in the superconducting (main) loop and the tunable Josephson junction is used to control the annealing process. For the coupler element, the tunable Josephson junction results in the tunable coupling strength Harris et al. (2009b).

The Hamiltonian of a SQUID with a CJJ loop is given by Han et al. (1989); Harris et al. (2010a); Boixo et al. (2016)

[TABLE]

where $\varphi=2\pi\Phi/\Phi_{0}=2e\Phi$ (we use $\hbar=1$) is the dimensionless magnetic flux in the main loop and $\varphi_{\mathrm{J}}$ the dimensionless magnetic flux in the CJJ loop with $\Phi_{0}$ denoting the magnetic flux quantum and $e$ the electron charge. In addition, $E_{C}$ and $E_{C_{\mathrm{J}}}$ are capacitive energies, $E_{L}$ and $E_{L_{\mathrm{J}}}$ are inductive energies, and $E_{J}$ is the Josephson energy. For an uncoupled SQUID, the inductive energy $E_{L}$ is given by $E_{L}=1/(4e^{2}L)$, where $L=L_{\mathrm{main}}+L_{\mathrm{J}}/4$ Harris et al. (2010a) is the total inductance. The potential of the flux variable $\varphi$, $V(\varphi)=E_{L}(\varphi-\varphi^{x})^{2}/2-E_{J}\cos(\varphi_{\mathrm{J}}/2)\cos(\varphi)$, can be either a single potential well or a double-well potential depending on the value of $\varphi_{\mathrm{J}}$. Thus, $\varphi_{\mathrm{J}}$ can be used to change the shape of the potential and also the barrier height between the double wells Han et al. (1989); Chiarello et al. (2005). This property is used to control the annealing process via the external flux $\varphi_{\mathrm{J}}^{x}$ Lanting et al. (2014) and to set the coupling strength of the coupling SQUID Harris et al. (2009b). The external flux $\varphi^{x}$ can be used to tilt the potential, thereby lowering one of the wells and raising the other one Chiarello et al. (2005). In terms of the qubit model (see Sec. IV), the external flux $\varphi^{x}$ can be used to set the parameter $h_{k}$ of the Ising Hamiltonian given in Eq. (4) Harris et al. (2010b).

III.1 Total Hamiltonian

So far, we have discussed the Hamiltonian of a single SQUID. In this section we introduce the Hamiltonian of the complete system consisting of the three SQUIDs. A tunable coupling constant $J_{jk}$ (see Eq. (4)) is realized by inserting a SQUID as a coupler element between the other two SQUIDs Harris et al. (2007, 2009b). For the SQUID that functions as the coupler element, we denote the flux in the CJJ loop by $\varphi_{\mathrm{J},0}$ and the one in the main loop by $\varphi_{0}$. Accordingly, energies that correspond to the coupler main loop are labeled by an index “[math]”, and those that correspond to the coupler CJJ loop by an index “J,[math]”. The external control flux $\varphi_{\mathrm{J},0}^{x}$ can be used to tune the coupling strength between the SQUIDs. We label the fluxes of the SQUIDs corresponding to the qubits with indices “$1$” and “$2$”, respectively. Since in the idealized model, the SQUIDs are equal and subject to the same annealing functions $A(s)$ and $B(s)$, their energies and their external fluxes $\varphi_{\mathrm{J}}^{x}$ are equal. Therefore, we drop the indices “1” and “2” in these cases. Although the external fluxes $\varphi_{1}^{x}$, $\varphi_{2}^{x}$, and $\varphi_{\mathrm{J}}^{x}$ depend on time, we do not write this explicitly for reasons of readability. A sketch of the complete system is shown in Fig. 2.

By coupling SQUIDs, the inductive energies change such that they are given by $E_{L}(1+M^{2}/(LL_{\mathrm{eff}}))$ for the SQUIDs representing the qubits and by $E_{L_{\mathrm{eff}}}=E_{L_{0}}L_{0}/L_{\mathrm{eff}}$ for the coupler, where $L_{\mathrm{eff}}=L_{0}-2M^{2}/L$ is the effective inductance of the coupler element and $M$ is the mutual inductance between the coupler and the other SQUIDs’ main loops. In addition to the modified Hamiltonians of the three SQUIDs, the interaction terms

[TABLE]

have to be included van den Brink et al. (2005). The tunable coupler can be operated without an external flux in the coupler main loop if the junction asymmetry is negligible, i.e., the difference between the critical currents (of the two junctions of the coupler SQUID) is much smaller than the sum of these critical currents Harris et al. (2009b). Since we do not consider junction asymmetries in the idealized model, we set $\varphi_{0}^{x}=0$. Collecting all terms, the total Hamiltonian is given by

[TABLE]

This is the final Hamiltonian for which we solve the TDSE without further simplification.

III.2 Effective coupling

The idea of inserting the coupler element is that it leads to a tunable effective coupling between the other two SQUIDs van den Brink et al. (2005); Harris et al. (2007, 2009b) such that the interaction Hamiltonian takes the form

[TABLE]

where $C(\varphi_{\mathrm{J},0}^{x})$ is the effective coupling strength tunable by the external flux $\varphi_{\mathrm{J},0}^{x}$ of the coupler CJJ loop.

To derive an approximate effective Hamiltonian that exhibits such a coupling term, we first replace the flux of the coupler CJJ loop $\varphi_{\mathrm{J},0}$ by its approximate expectation value. To obtain this expectation value, we expand the Hamiltonian of the SQUID given in Eq. (6) to second order in $\varphi_{\mathrm{J},0}-\varphi_{\mathrm{J},0}^{x}$ and set $\varphi_{0}=0$. The resulting Hamiltonian describes a shifted harmonic oscillator

[TABLE]

where $E_{L_{\mathrm{J},0}}^{\prime}=E_{L_{\mathrm{J},0}}+E_{J_{0}}\cos(\varphi_{\mathrm{J},0}^{x}/2)/4$. The expectation value of $\varphi_{\mathrm{J},0}$ in its ground state can thus be identified as

[TABLE]

With $\varphi_{\mathrm{J},0}$ replaced by $\langle\varphi_{\mathrm{J},0}\rangle$, we can find a matrix $T$ such that the transformation

[TABLE]

of the total Hamiltonian yields an effective Hamiltonian which contains an interaction term of the form Eq. (9). Choosing

[TABLE]

where the external fluxes $\varphi_{1}^{x}$ and $\varphi_{2}^{x}$ depend on time,

[TABLE]

and expanding

[TABLE]

to second order in (the products of) $\varphi_{0}$, $\varphi_{1}-\varphi_{1}^{x}$, and $\varphi_{2}-\varphi_{2}^{x}$, we obtain the effective Hamiltonian

[TABLE]

In the basis defined by the transformation Eq. (12), we obtain the term $C(\varphi_{\mathrm{J},0}^{x})\varphi_{1}\varphi_{2}$ where the dependence on $\varphi_{\mathrm{J},0}^{x}$ is given via $\langle\varphi_{\mathrm{J},0}\rangle$ in $\beta_{\mathrm{eff}}$ (see Eqs. (11) and (14)). The only coupling term between the coupler element and the other two SQUIDs that remains is the last term in Eq. (16) which is expected to be much smaller than the previous coupling terms since $E_{C}\ll E_{L}$.

Note that none of the approximations made to derive Eq. (16) affect the simulation results, as these are obtained by solving the TDSE for the Hamiltonian Eq. (8). However, as discussed in the next section, the approximation Eq. (16) is necessary to relate the external flux $\varphi_{\mathrm{J},0}^{x}$ to the coupling constant $J_{12}$ (denoted by $J$ for two qubits), which appears in the Ising Hamiltonian Eq. (4).

IV Mapping to the qubit model

In this section we investigate the mapping of the flux model Eq. (16) onto the qubit model Eq. (1) with $H_{\mathrm{final}}$ and $H_{\mathrm{init}}$ given by Eqs. (4) and (5), respectively. The two-qubit Hamiltonian reads

[TABLE]

As we will see below, reducing Eq. (16) to the generic form of Eq. (17) enforces a specific choice of the external fluxes $\varphi_{i}^{x}$ (see Eq. (21)) and gives a relation between $J$ and $\varphi_{\mathrm{J},0}^{x}$ (see Eq. (24)).

Since we have assumed the two SQUIDs to be identical, the mapping to the qubit model is the same for both, and therefore we omit the SQUID indices in this section. The two lowest-energy states $\ket{g}$ and $\ket{e}$ of each SQUID for $\varphi^{x}=0$ define the computational subspace Harris et al. (2010a). We obtain them by diagonalizing the first part in square brackets of Eq. (16) in $\varphi$- and $\varphi_{\mathrm{J}}$-space (see Sec. V for the definition of the discretized basis).

Note that the first summand given in Eq. (16) contains an effective change of the inductive energy depending on the value chosen for $\varphi_{\mathrm{J},0}^{x}$ (because $\beta_{\mathrm{eff}}$ depends on it). Therefore, the definition of the computational subspace changes with the coupling strength. This leads to slightly different annealing schemes, i.e., a dependence of $A(s)$ and $B(s)$ on $\varphi_{\mathrm{J},0}^{x}$, as observed experimentally Harris et al. (2009b, 2010a).

The computational basis states $\ket{\uparrow}$ and $\ket{\downarrow}$ are defined as the eigenstates of the operator $\varphi$ (and thus of the second part in square brackets in Eq. (16)) inside the computational subspace $\mathrm{span}\{\ket{g},\ket{e}\}$. We obtain

[TABLE]

where $|a|=|b|=1/\sqrt{2}$ Harris et al. (2010a) and $u(\varphi,\,\varphi_{\mathrm{J}})$ and $d(\varphi,\,\varphi_{\mathrm{J}})$ are the resulting amplitudes of the states $\ket{\uparrow}$ and $\ket{\downarrow}$ in $\varphi$- and $\varphi_{\mathrm{J}}$-space. Note that $\ket{g}$ and $\ket{e}$ depend on the time-dependent external flux $\varphi_{\mathrm{J}}^{x}$, implying that the definition of the computational states changes with time. The projection of the operator $\widetilde{E}_{L}\varphi$ with

[TABLE]

has eigenstates $\ket{\uparrow}$ and $\ket{\downarrow}$ with eigenvalues $\pm I_{p}(s)/2e$, respectively. Thus, in this subspace, $\widetilde{E}_{L}\varphi$ is represented by $I_{p}(s)\sigma^{z}/2e$ and the first contributions in square brackets in Eq. (16) are mapped to $-\Delta(s)\sigma^{x}/2$, where $\Delta(s)=E_{1}(s)-E_{0}(s)$ is the energy gap between the ground state $\ket{g}$ and the first excited state $\ket{e}$.

To derive the coupling terms, we write the SQUID indices $i$ again. For the terms in $\sigma_{i}^{z}$ and $\sigma_{1}^{z}\sigma_{2}^{z}$ to scale with the same annealing function $B(s)$ Harris et al. (2010b), $\varphi_{i}^{x}$ has to be set to

[TABLE]

where $\gamma=\mathrm{max}_{\varphi_{\mathrm{J},0}^{x}}\beta_{\mathrm{eff}}E_{L}^{2}/(1+\beta_{\mathrm{eff}})\widetilde{E}_{L}^{2}$. Disregarding the contribution of the last term in Eq. (16), we find that the Hamiltonian for $\varphi_{0}$ effectively decouples from the qubit Hamiltonian and thus the effective qubit Hamiltonian can be written as

[TABLE]

For all $J\in[-1,1]$, we have

[TABLE]

Thus, and because $E_{L}^{2}\beta_{\mathrm{eff}}(\varphi_{\mathrm{J},0}^{x})/(\widetilde{E}_{L}^{2}(\varphi_{\mathrm{J},0}^{x})(1+\beta_{\mathrm{eff}}(\varphi_{\mathrm{J},0}^{x})))$ is monotonic, it is possible to find $\varphi_{\mathrm{J},0}^{x}$ such that

[TABLE]

for all $J\in[-1,1]$, and Eq. (22) becomes

[TABLE]

which has the structure of an Ising model in a transverse field. Comparing Eq. (25) to Eqs. (4) and (5), we can identify $A(s)=\Delta(s)/2$ and $B(s)=\gamma I_{p}^{2}(s)M^{2}/L_{\mathrm{eff}}$ and see that

[TABLE]

where the last term only adds a small contribution since $M^{2}\ll LL_{\mathrm{eff}}$.

V Simulation

This section starts with a brief description of the numerical technique used to perform the simulation of the three-SQUID model. Then we discuss the choice of the model parameters that appear in Eq. (8) and explain the method by which we numerically extract the annealing scheme and the qubit-qubit coupling $J$, which appears in Eq. (17), from the dynamics of the fluxes.

For the simulation of the time evolution of the system defined by Eq. (8), we use the Suzuki-Trotter product-formula algorithm Suzuki (1985); De Raedt (1987) to numerically solve the TDSE

[TABLE]

The time-dependent Hamiltonian is discretized such that the state vector $|\psi(t)\rangle$ can be updated by a time step $\tau$ to $|\psi(t+\tau)\rangle$ using the time-evolution operator $U(t,t+\tau)=\exp(-i\tau H(t+\tau/2))$. To implement the algorithm, we fix a basis for the description of $\ket{\psi(t)}$ and a decomposition of the Hamiltonian $H(t)=\sum_{k}A_{k}(t)$ such that

[TABLE]

is a good approximation for sufficiently small $\tau$ and the update of the state vector can be performed with two-component updates only. For a detailed description of how to choose the $A_{k}$, see Ref. De Raedt (1987). In our simulation, we use the second-order approach given by

[TABLE]

Note that there is no need to diagonalize the Hamiltonian or to store the full matrices representing the Hamiltonian or the time-evolution operator.

For the description of the state $\ket{\psi}$, the fluxes $\varphi_{i}$ through the main loops are discretized, i.e., the wave function is defined at $\lambda_{i}$ discrete points $\varphi_{i\mathrm{min}}+l_{i}\Delta\varphi_{i}$, $l_{i}=0,\dots,\lambda_{i}-1$. By studying the convergence of the numerical results as a function of $\lambda_{i}$ and $\Delta\varphi_{i}$, we find that $\lambda_{1}=\lambda_{2}=47$ and $-2.0\leq\varphi_{1},\varphi_{2}\leq 2.0$, and $\lambda_{0}=31$ and $-1.0\leq\varphi_{0}\leq 1.0$ provide a good compromise between accuracy and computational work to solve the TDSE. Since $E_{L_{\mathrm{J}}}\gg E_{J}$ and $E_{L_{\mathrm{J,0}}}\gg E_{J_{0}}$, the Hamiltonian for $\varphi_{\mathrm{J},i}$ resembles an oscillator with small anharmonicity. In the harmonic-oscillator basis, the evolution of $\varphi_{\mathrm{J},i}$ can be well described with the three lowest states and thus the fluxes $\varphi_{\mathrm{J},i}$ through the CJJ loops can be discretized in the harmonic-oscillator basis and labeled by $\ket{m_{i}}$, $m_{i}=0,1,2$. In summary, the state $\ket{\psi}$ is represented by

[TABLE]

where the amplitudes $\phi_{l_{0},m_{0},l_{1},m_{1},l_{2},m_{2}}$ are stored as an array of $\lambda_{1}\lambda_{2}\lambda_{0}\times 3^{3}\approx 1.85\times 10^{6}$ complex double-precision numbers. To store this array, approximately $30\,\mathrm{MB}$ of memory is needed. Parallelization of the state updates is implemented using OpenMP. Testing with decreasing time steps $\tau$ and studying the convergence, we find that for $\tau=1.5\times 10^{-5}\,\mathrm{ns}$ the results are sufficiently accurate. Due to this small time step, one quantum annealing run of $5\,\mathrm{ns}$ takes about 16 hours on a 24-core node of the supercomputer JURECA Jülich Supercomputing Centre (2018).

V.1 Parameters

The parameters of the Hamiltonian Eq. (8) and the values of the time-dependent flux $\varphi_{\mathrm{J}}^{x}(s)$, which determines the annealing scheme, were provided to us by D-Wave Systems Inc. and are typical values of the D-Wave 2000Q processor D-Wave Systems Inc. . The device parameters used in the simulation are slightly modified and listed in Table 1 and $\varphi_{\mathrm{J}}^{x}(s)$ is plotted in Fig. 3. The external fluxes $\varphi_{i}^{x}$ and $\varphi_{\mathrm{J},0}^{x}$ are computed from Eq. (21) and by solving Eq. (24) numerically for $\varphi_{\mathrm{J},0}^{x}(J)$.

Using the provided parameters, we compute the annealing scheme of a single SQUID by exact diagonalization of Eq. (6) with $\varphi^{x}=0$ and computing $\Delta(s)$ and $I_{p}(s)$ as described in Sec. IV. However, the resulting annealing scheme (data not shown) does not match the data of the annealing scheme provided to us by D-Wave Systems Inc. D-Wave Systems Inc. (see Fig. 4 (dashed lines)). Better agreement between the two annealing schemes was found by using $E_{C}=4.68\,\mathrm{GHz}$ (which was computed from the provided capacitance directly) instead of $E_{C}=5.85\,\mathrm{GHz}$ (value provided by D-Wave Systems Inc.), $E_{L_{\mathrm{J}}}=54538\,\mathrm{GHz}$ instead of $E_{L_{\mathrm{J}}}=73388\,\mathrm{GHz}$, and $M=15.97\,\mathrm{pH}$ instead of $M=13.7\,\mathrm{pH}$, see Fig. 4 (solid lines). The disagreement in $B(s)$ for small $s$ could not be removed by slight variation of the model parameters. Changing $\varphi_{i}^{x}$ would reduce the disagreement for the single-qubit terms, but at the same time Eq. (21) would be violated, effectively yielding different functions for the single-qubit and two-qubit $\sigma^{z}$ terms. Thus, we decided to keep $\varphi_{i}^{x}$ as given by Eq. (21).

V.2 Estimation of the coupling strength and the annealing scheme

In order to map the full state Eq. (30) to the computational space, we have to trace out the degrees of freedom of the coupler element and project the resulting reduced density matrix onto the computational subspace. To do so, we discretize Eqs. (18) and (19) to obtain

[TABLE]

for a single qubit, and accordingly the product states for the two-qubit states, where $u_{lm}$ and $d_{lm}$ are the discretizations of $u(\varphi,\,\varphi_{\mathrm{J}})$ and $d(\varphi,\,\varphi_{\mathrm{J}})$, respectively. Since this projection is not a unitary transformation, the trace of the projected density matrix $\rho^{\mathrm{comp}}$ will be less than one if there is leakage to higher levels, i.e., if excitations to states outside the computational subspace occur. The deviation of the trace from one is a measure for the amount of leakage to higher levels.

To obtain the effective coupling strength and the annealing scheme we proceed as follows. We start with the ideal qubit Hamiltonian (for simplicity, with $h_{1}=h_{2}=0$)

[TABLE]

where $\Delta(s)$ and $C(s)$ are to be determined by comparison with the data obtained by simulating the model Eq. (8). For fixed $s\in[0,1]$, the evolution determined by $H_{2}(s)$ of the initial state $\ket{++}=\left(\ket{\uparrow}+\ket{\downarrow}\right)\otimes\left(\ket{\uparrow}+\ket{\downarrow}\right)/2$ and the expectation values $\langle\sigma_{1}^{\alpha}\sigma_{2}^{\beta}\rangle$ for $\sigma_{i}^{\alpha},\sigma_{i}^{\beta}\in\{\mathbb{I}_{i},\sigma_{i}^{x},\sigma_{i}^{y},\sigma_{i}^{z}\}$ in the evolved state can be calculated analytically. On the other hand, for any time $t$, we can compute these expectation values directly from the simulation of the time evolution of the initial state $\ket{{++}}$ expressed in flux degrees of freedom using Eqs. (31) and (32) at a fixed value for $s$. The time evolution is governed by the full Hamiltonian Eq. (8) based on the flux degrees of freedom with fixed $s$. In this case, the expectation values are computed by $\mathrm{Tr}(\rho^{\mathrm{comp}}\sigma_{1}^{\alpha}\sigma_{2}^{\beta})$. We can then estimate $\Delta(s)$ and $C(s)$, and thus the effective coupling strength and the annealing scheme, by fitting the analytical expressions to the simulation data.

VI Results

In this section, we present the results obtained from the simulation of the flux model described by the Hamiltonian Eq. (8). First, we show that the mapping between $J$ and $\varphi_{\mathrm{J,0}}^{x}$ given by Eq. (24) leads to the desired effective coupling strength. Subsequently, we discuss the effective annealing scheme obtained by using the procedure described in Sec. V.2. We check the results of the simulation based on the flux model and the results of the qubit model against each other by comparing the probabilities during and at the end of the annealing process. Finally, we briefly discuss the data obtained from the D-Wave quantum annealer in comparison to the simulation results.

VI.1 Effective coupling and annealing scheme

In order to assess the mapping between $J$ and $\varphi_{\mathrm{J,0}}^{x}$ using Eq. (24), we first study the effective mutual inductance $M_{\mathrm{eff}}$ as a function of $J$. We utilize the method described in Sec. V.2 for various values of $J$ and $s=1$ (such that $\Delta(s)\approx 0$) to determine the coupling strength. In this case, the analytical result for the expectation value $\langle\sigma_{1}^{y}\sigma_{2}^{z}\rangle=\sin(2C(1)t)$ can be used for fitting. The obtained value for $C(1)$ for each $J$ is then mapped onto the effective mutual inductance $M_{\mathrm{eff}}(J)=C(1)/I_{p}^{2}(1)=-J\gamma M^{2}/L_{\mathrm{eff}}$ and plotted against $J$. The result for the effective inductance $M_{\mathrm{eff}}(J)$ between the qubits is presented in Fig. 5 and shows good agreement between the theoretical linear curve from the approximation and the simulation result. For $J$ in the range $[-1,1]$, we can reach all values for $M_{\mathrm{eff}}$ in $[-M_{\mathrm{eff,max}},M_{\mathrm{eff,max}}]$ to good precision and have thus obtained a transformation $\varphi_{\mathrm{J,0}}^{x}\leftrightarrow J$ such that the mapping $J\leftrightarrow M_{\mathrm{eff}}$ is linear. Therefore, we can expect that the mapping onto the qubit model and the resultant mapping $J\leftrightarrow\varphi_{\mathrm{J,0}}^{x}$ work reasonably well.

To assess the effective annealing scheme, we use the method described in Sec. V.2 for different values $s^{\prime}\in[0,1]$, using the analytical expression of the expectation value

[TABLE]

for the fitting of $\Delta(s^{\prime})$ and $C(s^{\prime})$. Figure 6 shows the effective annealing scheme (data points) obtained in this way.

We find that the data points in Fig. 6 deviate from the annealing scheme for an uncoupled qubit (solid lines, obtained by using Eq. (6)), but they are in better agreement with the annealing scheme obtained by using $\widetilde{E}_{L}$ (see Eq. (20)) instead of $E_{L}$ (dashed lines). Note that for the computation of the annealing scheme, the single-SQUID Hamiltonian Eq. (6) is mapped onto the effective Hamiltonian $H^{\mathrm{eff},q}\approx-A(s)\sigma_{1}^{x}-B(s)h_{1}\sigma_{1}^{z}$. In the qubit model (Eq. (17)), this gives the same function for $B(s)$ as the term proportional to $\sigma_{1}^{z}\sigma_{2}^{z}$ in the case of two coupled qubits. Because of the choice for $\varphi_{1}^{x}$ (see Eq. (21)), this should also be the case for the SQUID model if the mapping to the qubit model works well enough.

We find, in agreement with our analytical calculation, that the effective coupling between the SQUIDs induces a shift in the inductive energy, leading to shifts in the annealing scheme. Including this shift, the effective annealing scheme can be well described by the single-SQUID annealing scheme. The influence of the coupling on the inductive energy was also observed in experiments Harris et al. (2009b).

VI.2 Comparison to the qubit model

The next step is to compare the overall performance and the final probabilities between the real-time simulation with the Hamiltonian given in Eq. (8) and the qubit Hamiltonian given in Eq. (17).

As mentioned in Sec. V, the amount of leakage to higher excited states can be computed by projecting the density matrix onto the computational subspace. The projected density matrix can also be used to obtain the probabilities of the computational basis states. As an illustration, in Fig. 7 we show the results for $J=-1$, $h_{1}=0.96$, and $h_{2}=0.94$. For this choice of parameters, the ground state of Eq. (4) is $\ket{{\uparrow\downarrow}}$. The total annealing time was set to $t_{a}=5\,\mathrm{ns}$ for the simulations of both the flux model and the qubit model.

This annealing time is much less than typically used on the D-Wave processors (order of microseconds), but for comparison of the results of the flux simulation with the qubit description only, this difference is unimportant. Note that in this section our aim is to scrutinize the validity of the qubit model as a description of the flux dynamics governed by the Hamiltonian Eq. (8), not to compare simulation results with experiments performed on the D-Wave quantum annealer (see Secs. VI.3 and VII below).

As seen from Fig. 7, there are small deviations from the probabilities obtained from the qubit representation. Some leakage which has its maximum at about $s=0.6$, where the change in the probabilities of the computational states is strongest, can also be observed. In general, the evolutions of both the full model and the qubit model show the same features.

In the following we refer to the probability of finding the system at the end of the annealing process in the ground state of Hamiltonian Eq. (4) as success probability. For the example case shown in Fig. 7, the success probability for the flux model and the qubit model differ only slightly.

In this example, the success probability is higher for the qubit model. However, Fig. 8 shows that there are also cases in which the success probability is lower for the qubit model. Note that the annealing process does not start with equal probability for all states because we start the annealing in the ground state of the system instead of in the state $\ket{++}=\ket{+}_{1}\otimes\ket{+}_{2}$, since for $B(s=0)>0$, the ground state of the qubit model Eq. (4) is not exactly the state $\ket{++}$ but a superposition of all basis states. A simulation of the qubit model comparing the annealing processes with the two different initial states shows deviations during the annealing process, but there is no significant difference in the success probability (data not shown).

In summary, we observed an influence of the coupling on the annealing scheme and some amount of leakage to higher levels. The important question, however, is whether these effects have consequences on the final success probability. Figure 8 shows the success probability for many different problems (defined in Appendix A) as a function of the minimal energy gap between the ground state and the first excited state during the annealing process, computed from the qubit model. As can be seen in Fig. 8, for most of the investigated cases, the effects on the success probability of using a subspace of a larger system as the qubit instead of an ideal qubit representation are rather small. The data points generated by the simulation based on Eq. (8) (black open circles) show, apart from a few exceptions, only small deviations from the data points of the ideal qubit model (red closed circles). Note that the success probability can be enhanced as well as reduced compared to the ideal qubit model.

Another interesting observation can be made in Fig. 8. For large minimal energy gaps $\delta E=\min_{s}E_{1}(s)-E_{0}(s)$, the system shows Landau-Zener behavior Landau (1932); Zener (1932). For small minimal energy gaps $\delta E$, the success probabilities form two clusters, one cluster of probability approximately equal to $0.5$ and a second one of probability approximately equal to $0.3$. Considering the energy spectra of these instances, we can separate them into three classes. The first class would be the one where the energies of all three excited states come close to the ground-state energy and the success probability clusters at roughly $0.3$. The second class of problem instances has a spectrum similar to the one shown in Fig. 9(a). Only the energy of the first excited state comes close to the ground-state energy. The energies of the second and third excited states are much higher. For this class, the success probability clusters around $0.5$. The third class of problems has a spectrum similar to the one shown in Fig. 9(b). The energies of the first and second excited states come close to the ground-state energy. For this class, the success probability depends on the particular problem instance. The reason is that for degenerate ground states, the probabilities to find these states are not necessarily equal Matsuda et al. (2009) and for problem instances which are close to these degenerate cases with unequal probabilities, this imbalance may have an influence when non-adiabatic transitions occur. All instances which show larger deviations between the success probabilities obtained from the qubit model and the flux model simulations (see Fig. 8) belong to the third class. A possible explanation for these deviations might be that, in some cases, due to the presence of the additional states in the flux model, these unequal probabilities in the computational subspace are different during the evolution than in the qubit model. We leave a more detailed study of this phenomenon for future work.

As noted in Ref. Albash et al. (2015); D-Wave Systems Inc. (2019) and confirmed by our analytical calculation in Sec. IV, there is some crosstalk between the qubits (last term in Eq. (26)), leading to small offsets in the parameters $h_{i}$. Furthermore, a dependence of the annealing scheme on the parameter $\varphi_{\mathrm{J},0}^{x}$ was found, also leading to small discrepancies between the ideal qubit representation and the full system. Additionally, for the mapping of $J$ to $\varphi_{\mathrm{J},0}^{x}$, we had to draw on an approximate analytical calculation which may be another source for the small differences between the results obtained from the two models. Nevertheless, the results fit very well. Interestingly, the coupler element, which can be viewed as part of the environment and might be the source of additional noise, does not cause significant deviations in the results compared to the results of the qubit description.

VI.3 Comparison to D-Wave 2000Q data

Because we find good agreement between the system described by the Hamiltonian Eq. (8) and the qubit model, we compare the success probability for both these systems with the percentage of successful runs on the D-Wave 2000Q quantum annealer. In Appendix A (see Table 3), we present the data obtained by at least ten repetitions of annealing simultaneously, 992 (976) copies of the two-qubit problems distributed over the D-Wave DW_2000Q_2 (DW_2000Q_2_1) chip, for an annealing time $t_{a}=20\,\mu\mathrm{s}$. Postprocessing and autoscaling have been turned off for all experiments on the D-Wave 2000Q.

We find that although the annealing time on the D-Wave is much larger than for our simulations ($20\mu$s instead of $5$ns), a large fraction of the D-Wave data seems to agree (approximately) with the corresponding success probabilities obtained from the simulation of the SQUID model Eq. (8) and its two-level approximation. This agreement is probably accidental. The annealing time of $5\,$ns was chosen to keep the real time to solve the TDSE of the SQUID model within acceptable limits as well as having some variation in the success probability at the end of the annealing process without having to use too-small values for, or differences between, the parameters $J$ and $h_{i}$. In spite of the large difference in annealing times, the good agreement suggests that in the D-Wave device there are physical processes at work that affect the annealing, processes which are not incorporated in the SQUID model Eq. (8) or the corresponding qubit model Eq. (26).

Concrete evidence for the relevance of such processes is shown in Table 2, where we present D-Wave data for three different cases whose energy spectra are shown in Figs. 9(a)–(c). Because the spectra of these cases differ significantly, we assume that they are a representative subset of the cases studied previously. Data characterizing the problem instances such as the minimal energy gap $\delta E$ and the problem gap $\Delta_{p}$ of the final Hamiltonian are listed in Table 2 as well as the frequency of runs finding the ground states ($\ket{\downarrow\uparrow},\ket{\uparrow\downarrow},$ and $\ket{\uparrow\downarrow}$ for the three cases, respectively) on the D-Wave machine for four different annealing times. The results reported in Table 2 were obtained by putting 992 copies of the two-qubit problems on the Chimera graph and performing 1000 annealing runs.

Table 2 shows some additional interesting facts. First, recall that for the shortest annealing time considered, i.e., $1\,\mu\mathrm{s}$, simulation of the quantum annealing process in the qubit description yields the ground state with probability one for the three cases listed. Clearly, as Table 2 shows, this is not the case for the D-Wave data, not even if we increase the annealing time to $1\,\mathrm{ms}$, as is most evident for case (c). We emphasize that these differences in the observed frequencies for finding the ground state are not due to poor statistics nor can they be attributed to the differences in the minimal gaps $\delta E$. Correlating these observations with the energy level spectra shown in Fig. 9 suggests that the observed reduction of the frequency for finding the ground state may be related to the distribution of energy levels for $s\rightarrow 1$. However, the mechanism that causes the observed change of frequencies when going from case (a) to (c) cannot be found within the description of the quantum dynamics of the two-qubit system. Explaining these experimental observations requires considering additional physical processes.

The first process that comes to mind is the interaction of the qubits with their environment at non-zero temperature. In the next section, we address this issue by solving the TDSE of the two-qubit model Eq. (26) coupled to a bath of two-level systems, complementing previous work that investigated the effects of finite temperature on quantum annealing Amin et al. (2008); Johansson et al. (2009); Amin et al. (2009b, a); Amin and Brito (2009); Dickson et al. (2013); Amin (2015).

VII System coupled to a bath

The annealing process of the isolated two-qubit system may be understood in terms of the adiabatic theorem. However, in the real world, the system modeling the two qubits is in contact with an environment at finite temperature. In this section we scrutinize the extent to which the coupling to a heat bath, i.e., the presence of thermal fluctuations, affects the annealing process. This, we hope, may shed light on the annealing behavior observed on the D-Wave machine in the cases studied.

We do not assume the heat bath to be Markovian but instead we solve the TDSE of the system comprising the two-qubit system and the bath. In order to be able to perform such simulations, it is necessary to keep these models simple. Therefore, it would be remarkable to find good quantitative agreement between the results of the simulations and those obtained with the D-Wave machine. Thus, the goal here is limited to find out if such models can reproduce, qualitatively, the trends observed in the D-Wave data.

We model the heat bath as a collection of two-level systems Phillips (1972); Anderson et al. (1972) which might represent e.g., defects in the material. Such models have been used to discuss noise and dephasing in superconducting resonators and circuits Shnirman et al. (2005); Müller et al. (2009); Cole et al. (2010); Burnett et al. (2014); Faoro and Ioffe (2015); Lisenfeld et al. (2015); de Graaf et al. (2018); Müller et al. (2019). We assume that this heat bath is at thermal equilibrium, with a temperature given by the operating temperature of the D-Wave machine. We denote the inverse of this operating temperature by $\beta^{*}=0.588\,\mathrm{ns}$ (in units of $\hbar=k_{B}=1$), corresponding to a temperature of $T\approx 13\,\mathrm{mK}$.

The Hamiltonian of the system (S) + bath (B) reads

[TABLE]

where $\lambda$ controls the overall strength of the system-bath interaction. The time evolution during the quantum annealing process of the closed quantum system defined by the Hamiltonian Eq. (35) is obtained by solving the TDSE given in Eq. (27) with the initial state

[TABLE]

The method to prepare the thermal state $|\Phi(\beta)\rangle$ and other technical details are discussed in Appendix B.

The system Hamiltonian is given by

[TABLE]

and changes with time as a function of the annealing variable $s=t/t_{a}$. We consider two extreme cases for $H_{\mathrm{B}}$ and $H_{\mathrm{SB}}$.

VII.1 Model I

In the first case, the bath is modeled as a ring of two-level systems represented by the Pauli matrices $\bm{\mu}_{n}=(\mu^{x}_{n},\mu^{y}_{n},\mu^{z}_{n})$, described by the Hamiltonian

[TABLE]

Here and in the following $N_{\mathrm{B}}$ denotes the number of bath particles. The couplings $r^{x}_{n}$, $r^{y}_{n}$, and $r^{z}_{n}$ are taken to be uniform random numbers in the range $[-1,+1]$ and $K$ determines the spectral range of $H_{\mathrm{B}}$. For random couplings it is unlikely that the model Eq. (38) is integrable (in the Bethe-ansatz sense) or has any other special features such as conserved magnetization. The bath Hamiltonian Eq. (38) with random couplings has the property that the distribution of nearest-neighbor energy levels is Wigner-Dyson-like Zhao et al. (2016). Extensive simulation work on spin baths with very different degrees of connectivity Jin et al. (2013) suggests that as long as there is randomness in the system-bath coupling and randomness in the intrabath coupling, the simple model Eq. (38) may be considered as a generic spin bath Zhao et al. (2016). The Hamiltonian describing the interaction of the two-qubit system with the bath is taken to be

[TABLE]

where $n$ and $m$ are chosen randomly from the set $\{1,\ldots,N_{\mathrm{B}}\}$ such that $n\not=m$. The $r^{\alpha}_{n,1}$ and $r^{\alpha}_{m,2}$ are real-value random numbers in the range $[-1,+1]$.

VII.2 Model II

In this case, the bath is modeled as a collection of non-interacting two-level systems given by the Hamiltonian

[TABLE]

where the parameter $\Omega$, together with the random numbers $r_{n}^{x},\,r_{n}^{y},\,r_{n}^{z}\in[-1,1]$, determines the level splitting of each two-level system. The interaction between the qubits and the two-level systems of the bath is chosen such that each qubit interacts with a different half of the bath. This type of interaction is modeled by the Hamiltonian

[TABLE]

VII.3 Parameters

Obviously, to compare the simulation results with D-Wave results it is necessary to solve the TDSE for the same time interval as used on the D-Wave machine. In practice, this requirement puts a severe constraint on the size of the problems for which one can solve the TDSE within a reasonable amount of real time. Simulating the annealing process over $1\,\mu\mathrm{s}$ (the shortest annealing time possible on the D-Wave machine) for a system comprising $N_{B}=16$ on a BullSequana X1000 supercomputer (JUWELS Jülich Supercomputing Centre (2019)) takes about 4 hours using 40 compute cores. Performing a simulation for $N_{B}=28$ two-level systems in the bath takes about one week (400000 time steps of $0.0025\,\mathrm{ns}$ using 5120 compute cores). Earlier work which studied the TDSE dynamics of two spins coupled to a spin bath De Raedt et al. (2017) shows that the results for $N_{B}=16$ and $N_{B}=28$ primarily differ in the statistical fluctuations on the data (see also Appendix B). Therefore, we have chosen to perform most simulations with $N_{B}=16$, repeating runs with different random numbers to collect statistics, and use a few runs with $N_{B}=28$ as an additional check on the data.

We use the annealing schedule of the DW_2000Q_2 chip (see Fig. 10) which is the machine that we used for our experiments.

In the case of model I, the initial state of the bath is prepared by projection with the operator $\exp(-\beta H_{B})=\exp(-\beta K(H_{B}/K))$, as explained in Appendix B. From Eqs. (43) and (44) it is clear that baths with the same $\beta K$ (and the same $r_{n}^{x}$, $r_{n}^{y}$, and $r_{n}^{z}$; see Eq. (38)) will have the same thermal equilibrium properties. Therefore, we may use $K$ as an adjustable parameter to “calibrate” the temperature of the bath with respect to the operating temperature of the D-Wave machine on which we performed our experiments. For instance, for a fixed choice of $r_{n}^{x}$, $r_{n}^{y}$ and $r_{n}^{z}$, baths with $(K=5\,\mathrm{GHz},\,\beta=0.2\,\mathrm{ns}\,(T\approx 38.2\,\mathrm{mK}))$ and $(K=5/3\,\mathrm{GHz},\,\beta=0.6\,\mathrm{ns}\,(T\approx 12.7\,\mathrm{mK}))$ have the same thermal properties. On the other hand, $K$ sets the time scale of the dynamics of the two-level systems of the bath. Simulations with $K$ in the range $[1\,\mathrm{GHz},5\,\mathrm{GHz}]$ (data not shown) reveal that the primary quantity of interest, the success probability of the two-qubit system at $s=1$, does not change significantly with $K$ (in the mentioned range and for the same value of $\lambda$). This leaves only the system-bath interaction $\lambda$ as a parameter to fit the simulation data to the D-Wave data.

In the case of model II, $\Omega$ plays the role of $K$ in model I, i.e., $\beta\Omega$ determines the thermal equilibrium properties of the bath. Note that for modest values of $N_{B}$, model II is too simple to act as a genuine heat bath, but as a model for a few defects interacting with the SQUIDs, it can be a realistic choice Shnirman et al. (2005); Müller et al. (2009); Cole et al. (2010). Therefore, in this case, we set $\beta\approx\beta^{*}$ and use $\Omega$ and $\lambda$ as fitting parameters.

VII.4 Comparison to D-Wave 2000Q data

Figures 11(a)–(c) depict the results of the simulation with the heat bath (closed circles) averaged over ten different random initializations of the heat bath with $N_{B}=16$. Results for $N_{B}=28$ (asterisks) show that for each value of $\beta$, the averages of 10 samples of $N_{\mathrm{B}}=16$ data are in good agreement with the data obtained from one $N_{\mathrm{B}}=28$ sample. The solid line indicates the probability $p_{0}$ of finding the ground state of the isolated qubit system in thermal equilibrium for $H_{\mathrm{S}}$ at $s=1$, i.e., $p_{0}=\exp(-\beta E_{0})/Z$, where $E_{0}$ is the ground-state energy and $Z=\mathrm{Tr}(\exp(-\beta H_{S}))$ is the partition function. Qualitatively, the simulation data obtained using model I (closed circles) nicely match the equilibrium line. The deviations from the equilibrium line may be due to the freeze-out where thermal transitions stop Amin (2015) and/or too-short annealing times and/or the magnetic Foehn effect Saito and Miyashita (2001). However, the simulation data do not match the data generated on the D-Wave machine (crosses). Assuming that the qubit system on the D-Wave machine is in thermal equilibrium, we would infer from Figs. 11(a)–(c) that $\beta\approx 0.2\,\mathrm{ns}$ (corresponding to $T\approx 38.2\,$mK), which is about a factor of 3 smaller than the inverse operational temperature of about $\beta^{*}=0.588\,\mathrm{ns}$ (corresponding to $T\approx 13\,$mK).

For model II, we have searched the parameter space $1/8\,\mathrm{GHz}\leq\Omega\leq 2\,\mathrm{GHz}$ and $1/2\,\mathrm{GHz}\leq\lambda\leq 2\,\mathrm{GHz}$ for sets of ($\Omega,\lambda$) which would fit the D-Wave results best. These data are shown in Figs. 11(a)–(c) as open circles ($\Omega=0.125\,\mathrm{GHz},\lambda=1\,\mathrm{GHz}$) and open squares ($\Omega=0.125\,\mathrm{GHz},\lambda=1.5\,\mathrm{GHz}$). For cases (b) and (c) (see Table 2), the former fit remarkably well to the D-Wave data. However, we have not found a set ($\Omega,\lambda$) which fits all D-Wave data very well.

VIII Summary

We simulated the full system of three SQUIDs resembling two qubits and one tunable coupler element as used in the quantum annealer built by D-Wave Systems Inc. and found that the higher energy levels as well as the presence of the coupler element have observable effects on the annealing process which however do not have a strong influence on the final success probability compared to the ideal qubit model. In contrast to the investigation of the influence of the higher levels and resonators present in current systems for gate-based quantum computing Willsch et al. (2017), we found that, apart from a few exceptions, the effects are small for the cases of quantum annealing examined.

Furthermore, we investigated in detail the derivation of the qubit representation to obtain expressions for $\varphi_{i}^{x}$ and $J$ that led to satisfying results (see Figs. 5–8). The simulation results are in good agreement with effects encountered in this derivation such as the change in the annealing scheme depending on the choice of $\varphi_{\mathrm{J},0}^{x}$, which is also supported by findings in experiments Harris et al. (2009b).

Simulation data for the SQUID model Eq. (8) and the corresponding two-level approximation Eq. (26) for an annealing time $t_{a}=5\,\mathrm{ns}$ show remarkably good agreement with the D-Wave data obtained with an annealing time $t_{a}=20\,\mu\mathrm{s}$, also in those cases for which the success probability is substantially less than one. Although this agreement might be accidental, it suggests that the dynamics of the D-Wave quantum annealer are more complicated than what can be described by a closed-system model of the SQUIDs including higher levels and the tunable coupler.

Therefore, as a first step, we have studied quantum annealing in the presence of a heat bath. We solve the TDSE of the two-qubit system (Eq. (26)) plus bath (Eq. (38)) for three cases with qualitatively different energy spectra of the two-qubit system (see Fig. 9). Comparing D-Wave and simulation results for an annealing time of $1\,\mu\mathrm{s}$, we found that the simulation data for the success probabilities of the two-qubit systems quite nicely agree with the corresponding thermal equilibrium values but also that these probabilities are significantly larger than those obtained with the D-Wave annealer. We have not found a common set of parameters $(\beta,K,\lambda)$ of the two-qubit-bath model that reproduces the D-Wave results for the three different cases considered.

Modeling the environment as a collection of non-interacting two-level defects (see Eq. (40)) was found to yield a much more appropriate description of the D-Wave data. Although we could readily find values of the bath parameters $\Omega$ and $\lambda$ for which the solution of the TDSE yields results that are in excellent agreement with D-Wave data for two of the three different cases considered, we could not find a similar level of agreement with these data for all three cases if we impose the elementary requirement that the bath parameters $\Omega$ and $\lambda$ do not depend on the $J$ and $h$ that define the problem Hamiltonian.

Unlike in the case of a time-independent problem, the exchange of energy between the two qubits in the time-dependent (annealing) field and the bath of two-level systems seems to be an intricate process which, according to our simulation data, depends on the model parameters in a complicated manner. We leave a detailed study of this challenging problem for future research.

Acknowledgements.

We would like to thank Mauricio Reis of D-Wave Systems Inc. for providing us with technical information and Seiji Miyashita for helpful discussions. Access and compute time on the D-Wave machine located at the headquarters of D-Wave Systems Inc. in Burnaby (Canada) were provided by D-Wave Systems Inc. The authors gratefully acknowledge the Gauss Centre for Supercomputing e.V. for funding this project by providing computing time on the GCS Supercomputer JUWELS at Jülich Supercomputing Centre (JSC) and the computing time granted through JARA on the supercomputer JURECA at Forschungszentrum Jülich. D.W. was supported by the Initiative and Networking Fund of the Helmholtz Association through the Strategic Future Field of Research project “Scalable solid state quantum computing” (ZT-0013).

Appendix A Data

Table 3 contains a list with the parameters $h_{i}$ and $J$, the minimal energy gap $\delta E$, and the success probabilities for the problems used to generate Fig. 8.

Appendix B Numerical Solution of the TDSE

The numerical solution of the TDSE for a pure state of $N_{\mathrm{B}}+2$ two-level systems requires computational resources (memory and CPU time) proportional to $2^{N_{\mathrm{B}}+2}$. For a brute force calculation of thermal expectation values $\mathbf{Tr}(e^{-\beta\mathcal{H}}{\cal A}(t))/\mathbf{Tr}(e^{-\beta\mathcal{H}})$ this number changes to $2^{N_{\mathrm{B}}+2}\times 2^{N_{\mathrm{B}}+2}$. Fortunately, this increase in cost can be avoided by making use of random-state technology, reducing the cost to that of solving the TDSE for one pure state Hams and De Raedt (2000). If $|\Phi\rangle$ is a pure state, picked uniformly from the $D=2^{{N_{\mathrm{B}}+2}}$-dimensional unit hypersphere, one can show that for Hermitian matrices $X$ Hams and De Raedt (2000),

[TABLE]

and that the statistical errors resulting from approximating $\mathbf{Tr}(X)$ by $D\langle\Phi|X|\Phi\rangle$ are small if $D$ is large Hams and De Raedt (2000). For large baths, this property of the random pure state renders the problem amenable to numerical simulation.

In the case at hand, we proceed as follows. First, we generate a thermal random state of the bath system, meaning that we compute the pure state

[TABLE]

where $\beta$ denotes the inverse temperature. For any bath observable ${\cal A}(t)$ we have Hams and De Raedt (2000)

[TABLE]

The initial state of the whole system is then a product state of the ground state of the two qubits at $s=0$ and the thermal pure state $|\Phi(\beta)\rangle$, i.e.,

[TABLE]

Since the Hamiltonian Eq. (35) depends explicitly on time, we can only solve Eq. (27) numerically by time stepping. For this purpose, we use a Suzuki-Trotter product-formula-based algorithm De Raedt and Michielsen (2006). This algorithm employs a decomposition in terms of unitary matrices and is unconditionally stable by construction. All our simulations for the two qubits coupled to a heat bath were carried out with the massively parallel quantum-spin dynamics simulator (in-house software), which is based on the same computational kernel as the massively parallel quantum computer simulator De Raedt et al. (2019).