Experimental and Theoretical Study of Thermodynamic Effects in a Quantum Annealer

TL;DR

This paper combines experiments and simulations to analyze thermodynamic effects in a quantum annealer, revealing insights into noise effects, state transfer mechanisms, and potential improvements in sampling uniformity.

Contribution

It introduces a novel approach that interpolates between quantum dynamics and thermodynamics, providing new insights into noise effects and state transfer in quantum annealers.

Findings

Enhanced thermodynamic signatures with anneal pause parameter

Identification of a special eigenstate bridging separated ground states

Potential to improve sampling uniformity in quantum annealing

Abstract

Quantum devices are affected by intrinsic and environmental noises. An in-depth characterization of noise effects is essential for exploiting noisy quantum computing. To this end, we studied the energy dissipative behavior of a quantum annealer via experiments and numerical simulations. Our investigation adopts a recently proposed technique that interpolates between pure quantum dynamics and pure thermodynamics. Experiments were conducted on a quantum annealer with an anneal pause function, which inserts a thermal relaxation period into the annealing schedule by pausing the transverse field, which is a source of quantum fluctuation. After investigating the special Hamiltonian that characterizes the quantum thermodynamics of the system, we then observed enhancement of thermodynamic signature depending on the anneal pause parameter. The time development of the state vector, observed in the open quantum simulation, provides rich information for investigating phenomena beyond energy-gap analysis. We identified a special eigenstate bridges ground states far-separated in Hilbert space and the transfer probabilities from one ground state to another. This finding can improve the sampling uniformity by reducing the sampling bias in finding the classical ground states in the quantum annealer. Our study does not only characterize the open quantum phenomenon of the specific Hamiltonian but also demonstrates the usefulness of the method in investigating noisy quantum devices.

Equations

$$\frac{d}{dt}|\psi (t)⟩=-i\mathcal{H}(t)|\psi (t)⟩,$$

$$\mathcal{H}(t)=B(t){\mathcal{H}}_c+A(t){\mathcal{H}}_q.$$

$${\mathcal{H}}_c=-\sum _{(ij)}{J}_{ij}{\sigma }_i^z{\sigma }_j^z-\sum _i{h}_i{\sigma }_i^z,$$

$${\mathcal{H}}_q=-\sum _i{\sigma }_i^x.$$

$$\frac{d{P}_i(t)}{dt}=\sum _j{\mathcal{L}}_{ij}{P}_j(t),$$

$${\mathcal{L}}_{ij}=\begin{cases}\ \Big{[}1+e^{\frac{E_{i}-E_{j}}{T}}\Big{]}^{-1}&(\mbox{single-spin flip})\\ \ -\sum_{k\neq i}{\cal L}_{ki}&(i=j)\\ \ 0&(\mbox{otherwise})\end{cases},$$

$$|\psi (t)⟩=\sum _i{c}_i(t)|i⟩.$$

$$c_i(0)$$

$$P_i(0)$$

$$|\tilde{\psi }(t+dt)⟩$$

$$=\sum _i\sqrt{r_i(t+dt)}\frac{c_i(t+dt)}{|c_i(t+dt)|}|i⟩,$$

$${\tilde{P}}_i(t+dt)$$

$$|{\tilde{c}}_i(t+dt){|}^2={\tilde{P}}_i(t+dt)$$

$$\ket{\psi(t+dt)}=\\ \sum_{i}\sqrt{(1-\alpha)\bigl{|}c_{i}(t)-i\sum_{j}\mathcal{H}(t)_{ij}c_{j}(t)dt\bigr{|}^{2}+\alpha\Bigl{(}|c_{i}(t)|^{2}+\sum_{j}\mathcal{L}_{ij}|c_{j}(t)|^{2}dt\Bigr{)}}\\ \times\frac{c_{i}(t)-i\sum_{j}\mathcal{H}(t)_{ij}c_{j}(t)dt}{|c_{i}(t)-i\sum_{j}\mathcal{H}(t)_{ij}c_{j}(t)dt|}\ket{i}.$$

$$\ket{\psi(t+dt)}=\sum_{i}c_{i}(t)\ket{i}+\mathcal{F}\big{(}\{c_{i}(t)\},\mathcal{H}(t),\mathcal{L},\alpha\big{)}dt+O(dt^{2}),$$

$$\ket{\psi(t+dt)}-\ket{\psi(t)}=\mathcal{F}\big{(}\{c_{i}(t)\},\mathcal{H}(t),\mathcal{L},\alpha\big{)}dt+O(dt^{2}).$$

$$\frac{d}{dt}\ket{\psi(t)}=\mathcal{F}\big{(}\{c_{i}(t)\},\mathcal{H}(t),\mathcal{L},\alpha\big{)}.$$

$$C$$

$$S$$

$$P_c$$

$$P_s$$

$$\mathcal{H}={\mathcal{H}}_c+{\mathcal{H}}_q=\left(\begin{array}{cc}-h& -\Gamma \\ -\Gamma & h\end{array}\right)$$

$$\mathcal{L}=\frac{1}{e^{\beta h}+e^{-\beta h}}\left(\begin{array}{cc}-e^{-\beta h}& e^{\beta h}\\ e^{-\beta h}& -e^{\beta h}\end{array}\right),$$

$$\mathcal{L}=\left(\begin{array}{cc}0& 1\\ 0& -1\end{array}\right).$$

$$|\psi (t)⟩=\left(\begin{array}{c}u(t)\\ d(t)\end{array}\right)\text{mbox}and\mathbf{P}(t)=\left(\begin{array}{c}|u(t){|}^2\\ |d(t){|}^2\end{array}\right)$$

$$u(t+dt)=$$

$$\displaystyle\times\frac{u(t)+\big{[}ihu(t)+i\Gamma d(t)\big{]}dt}{\big{|}u(t)+\big{[}ihu(t)+i\Gamma d(t)\big{]}dt\big{|}}$$

$$d(t+dt)=$$

$$\displaystyle\times\frac{d(t)+\big{[}i\Gamma u(t)-ihd(t)\big{]}dt}{\big{|}d(t)+\big{[}i\Gamma u(t)-ihd(t)\big{]}dt\big{|}}$$

$$u(t+dt)-u(t)$$

$$d(t+dt)-d(t)$$

$$\frac{d}{dt}u(t)$$

$$\frac{d}{dt}d(t)$$

$$\displaystyle\frac{d}{dt}\big{[}u(t)u^{*}(t)\big{]}=\frac{d}{dt}\rho_{11}$$

$$\displaystyle\frac{d}{dt}\big{[}u(t)d^{*}(t)\big{]}=\frac{d}{dt}\rho_{12}$$

$$(\ast \text{mbox}Definitionofthesystem\ast )$$

Full text

I Introduction

Quantum mechanics is sufficiently powerful to explore all possible combinations of variables in quadratic unconstrained binary optimization (QUBO) problems. To exploit this power, researchers have investigated quantum annealing (QA) in the Ising model Kadowaki and Nishimori (1998); Kadowaki (1998); Brooke et al. (1999); Farhi et al. (2001); Santoro (2002); Santoro and Tosatti (2006); Das and Chakrabarti (2008); Morita and Nishimori (2008); Tanaka et al. (2017). In the framework of the Ising model, QUBO problems can be represented with spin-1/2 particles as the binary variables, whose arbitrary relations are encoded as $p$-body ($p\geq 2$) interactions and their local fieldsBillionnet and Jaumard (1989); Biamonte (2008); Lucas (2014). Using a hardware implementation of QA developed by D-Wave Systems Inc., in which the spin-1/2 variables are superconducting flux qubits, we can evaluate algorithms for QUBO problems in real-world settingsJohnson et al. (2011).

Reducing the thermal fluctuations and other noise sources would lengthen the coherence time and improve the reliability of quantum information processing on qubits, but is challenging in practice. At the same time, the performance to solve QUBO problems by QA can be improved by thermal relaxation after anticrossing, which recovers the ground state probability from excited statesAmin et al. (2008); Venuti et al. (2017); Dickson et al. (2013). Analyzing open quantum systems from this perspective would provide new insights, not only for improving the coherence time and evaluating device limitations but also for understanding the synergy between quantum and thermal effects, thereby improving device performance.

Open quantum systems have been studied by perturbation approaches such as the Redfield and Lindblad equationsRedfield (1965); Kossakowski (1972), which handle the microscopic interactions between the system and the bath. Alternatively, we have analyzed open quantum systems by interpolating between quantum dynamics and thermodynamicsKadowaki (2018). In the interpolated dynamics (ID), the coupling strength between the system and its environment is parameterized by an interpolation ratio, negating the need for implementing the system–bath interactions. This method can analyze the system under any parameter condition. At one extreme of the interpolation, ID recovers the closed quantum dynamics; at the other extreme, ID provides the classical thermodynamics. Between these two extremes, the dynamics are mixed. When the interpolation ratio is small (i.e., the coupling is weak), the system is dominated by quantum dynamics and is perturbed by thermodynamics. The dynamics of a two-level system derived by this method generalize the optical Bloch equationsBloch (1946); Arecchi and Bonifacio (1965). At the weak-coupling and asymptotic limits, the Bloch equations and our model equations become equivalent.

Differential equations of the dynamics are derived from a continuous-time limit of repeated mixing process of the two solutions from quantum dynamics and thermodynamics. This methodology is similar to stochastic Shrödinger equations (SSE) by taking a continuous-time limit of repeated quantum measurements of the environmentPellegrini (2008, 2010). In SSE, a system interacts with “a copy” of the environment for a short period, and quantum measurement of the environment is conducted and repeats the same process with “other copies” of the environment. The measurement of the environment affects the system indirectly by a projection as they have interacted before the measurement. Note that the solution is non-deterministic due to quantum jumps (stochastic projections introduced by quantum measurements). Solutions of the equations are called quantum trajectories. Continuous-time limit of the equations results in differential equations, and depending on the measurement type, either jump-type or diffusive-type differential equations are obtained.

While the system is continuously interfered by the indirect continuous measurement in SSE, another continuous interference is adopted in ID. Here, we briefly review ID, the method for open quantum systemsKadowaki (2018). We construct ID from Schrödinger equation and classical master equation, which is determined by diagonal elements of the Hamiltonian in the Schrod̈inger equation. Schrödinger equation in the natural unit ($\hbar=1$) is

[TABLE]

where the time-dependent Hamiltonian is defined by

[TABLE]

In QA, we adopt ${\mathcal{H}}_{c}$ is a problem Hamiltonian representing QUBO problems using diagonal operator $\sigma^{z}$,

[TABLE]

and ${\mathcal{H}}_{q}$ is represented by off-diagonal operator $\sigma^{x}$,

[TABLE]

Time dependence of the Hamiltonian is determined through $A(t)$ and $B(t)$. These scheduling parameters $A(t)$ and $B(t)$ change monotonically from 0 to 1 and from 1 to 0 respectively. Therefore, the initial Hamiltonian is $\mathcal{H}_{q}$, which is the local transverse field, and the ground state is a superposition of all possible states. The final Hamiltonian is $\mathcal{H}_{c}$ at the end of the dynamics, and the ground state realizes the solution of the QUBO problem. The actual curves of the parameters in this study are depicted in Fig. 1. $A$ and $B$ are functions of not $t$ directly but indirectly through $s(t)$, where $0\leq s\leq 1$. We can design the shape of $s(t)$.

Master equation describes time-development of a probability distribution,

[TABLE]

where $P_{j}(t)$ is the probability of the $j$-th state and ${\mathcal{L}}_{ij}$ is the transition rate matrix. This matrix is characterized by the temperature $T$ and energy levels $\{E_{i}\}$ of ${\mathcal{H}}_{c}$,

[TABLE]

where we use the natural unit ($k_{B}=1$).

Solutions for Schrödinger equation can be expressed as a linear combination of basis states with complex coefficients,

[TABLE]

We assume that both dynamics start from a compatible states, a trivial ground states of the given initial Hamiltonian, i.e. a superposition of all possible $z$-basis states for Schrödinger equation and uniform distribution for master equation,

[TABLE]

where $N$ is the number of all possible states. This assumption satisfies the relation $|c_{i}(0)|^{2}=P_{i}(0)$ for any $i$.

Similar to the short-time development by the total Hamiltonian in SSE, the two systems develop respectively for a short period. Although SSE introduces the system-environment interaction in the total Hamiltonian, ID introduces not directly but indirectly as a result of the latter part of this methodology. To interact the two systems indirectly, we construct a new state vector based on the solution of Schrödinger equation with the solution of master equation as a source of thermal interference. The new state vector with a parameter $\alpha$, which controls the strength of the thermal interference, and an associated probability distribution are defined by

[TABLE]

where $r_{i}(t)=(1-\alpha)|c_{i}(t)|^{2}+\alpha P_{i}(t)$. With these definitions, the same relation assumed in the initial state also satisfies

[TABLE]

for any $i$. Through $r_{i}(t)$, the new state vector $\ket{\tilde{\psi}(t)}$ carries all information needed for ID. As we interpolate solutions between Schrödinger and master equations parametrized by $\alpha$, the new state vector becomes a solution of Schrödinger equation in the case of $\alpha=0$, and master equation in $\alpha=1$.

We adopt $\ket{\tilde{\psi}(t+dt)}$ and $\{\tilde{P}_{i}(t+d)\}$ as a state vector and a probability distribution at $t+dt$ insted of $\ket{\psi(t+dt)}$ and $\{P_{i}(t+d)\}$ for further time-development by Schrödinger equation and master equation. This artificial “jump” form $\ket{\psi}$ to $\ket{\tilde{\psi}}$ (and from $\{P_{i}\}$ to $\{\tilde{P}_{i}\}$) is correspond to the stochastic projection by measurement in SSE. While SSE is non-deterministic, ID is deterministic. Difference equations (discrete jump process) can be obtained by substituting $c_{i}(t+dt)=c_{i}(t)-i\sum_{j}\mathcal{H}_{ij}c_{j}(t)dt$ from Eq. (1) and $P_{i}(t+dt)=P_{i}(t)+\sum_{j}\mathcal{L}_{ij}P_{j}(t)dt$ from Eq. (5) for Eq. (10), and applying the assumption, i.e, replacing $\ket{\tilde{\psi}(t+dt)}\to\ket{\psi(t+dt)}$ and $P_{i}(t)\to|c_{i}(t)|^{2}$. Then, we have

[TABLE]

A series expansion for $dt$ of this equation is

[TABLE]

where $\mathcal{F}(\cdots)$ is the first order coefficient. This can be written as a difference equation,

[TABLE]

Continuous-time limit gives differential equations of ID,

[TABLE]

SSE and ID can be expressed in a wave function form as well as a density matrix form, because the system keeps pure statePellegrini (2008, 2010). We reported the density matrix representation of ID for a two-level system Kadowaki (2018), and a wave function representation is provided in Appendix A. The differential equations (30) and (31) consist of terms from Shrödinger equation (the first and the second terms) and master equation (the third term) as well as additional non-linear terms in the right-hand side (RHS). The thermal bath is included in the dynamics not explicitly but indirectly through master equation. As a result, system–bath interaction is introduced and the interpolation parameter $\alpha$ controls the strength of energy dissipation. Continuous interference of a wave function by replacing with an interpolated state in our proposed dynamics and quantum jumps by measurements in stochastic Schrödinger equation have similar roles in terms of introducing thermal (environmental) effects.

In the present paper, we conduct anneal pause experiments on the eight-spin quantum-signature modelBoixo et al. (2013); Albash et al. (2015) depicted in Fig. 2 using the D-Wave quantum annealer. Here we consider the ferromagnetic case, and set all the interactions as $J_{ij}=1$. Due to the frustration caused by competition between the ferromagnetic interaction and the local field, the ground states degenerate into a cluster of 16 ground states with four up-spins in the core, and an isolated ground state with all eight core and outer spins being down, as follows:

[TABLE]

where four spins in the left(right) side represents core(outer) spins and $\updownarrow$ stands for both of up ($\uparrow$) and down ($\downarrow$) spin configurations. Probabilities of the clustered and the isolated states are defined by

[TABLE]

Note that $P_{c}$ is the average of 16 states while $P_{s}$ is calculated from a single state. The original studyBoixo et al. (2013) found that the isolated state probability is suppressed in quantum dynamics and enhanced in thermodynamics. A ratio $P_{s}/P_{c}$ is 1 at the thermodynamic equilibrium. In other words, all ground states can be obtained with equal probability. If the ratio is lower (higher) than 1, the isolated (clustered) state is suppressed and less probability. Therefore the ratio reflects the imbalance between two types of ground states. In the case that the ratio is close to 1, all possible ground state configurations can be obtained efficiently, and the number of annealing trials can be minimized. To this end, understanding the phenomena in the model helps to mitigate the suppression of the isolated state probability, and the ultimate goal is to figure out a general methodology to achieve $P_{s}/P_{c}\sim 1$. This equal probability sampling is called “fair sampling”Matsuda et al. (2009); Chancellor (2017); Mandrà et al. (2017); Könz et al. .

The anneal pause depicted in Fig. 3 is a special annealing schedule that pauses the change of the transverse field for a specified time, in other words, the waiting period. During the anneal pause period, the quantum mechanics keep the system essentially in the same state as the Hamiltonian does not change. Therefore, thermodynamics dominate despite the low temperature and sufficient isolation from the environment. By changing the pausing level of the transverse field, we can change the Hamiltonian during the period, thereby controlling the thermodynamics. Numerical simulations are expected to probe the detailed dynamics of the above-described open quantum system, beyond merely analyzing the energy-gap structure. In particular, we investigate how the probability transfer from the clustered ground states to the isolated ground state mitigates the bias in the classical ground-state probabilities.

The remainder of this study is organized as follows: Section II describes the setup and results of the D-Wave experiments of the quantum-signature model. Section III provides the analytical expression of the differential equations of the model by expanding the interpolation method, performs the numerical simulations, and conducts data analysis of the dynamics. Finally, Section IV summarizes and discusses our results.

II D-Wave experiments

The anneal time $\tau_{\mathrm{anneal}}$ plus anneal pause period $\tau_{\mathrm{pause}}$ settings in the anneal pause experiments were varied as 1 + 1, 2 + 2, 5 + 5, and 10 + 10 $\mu$s. Under each condition, 49 anneal pause parameters ($s_{\mathrm{pause}}=$ 0.02, 0.04, $\dots$, 0.98), normalized by the anneal time, were tested. For example, $s_{\mathrm{pause}}=$ 0.46 during $\tau_{\mathrm{anneal}}$ = $\tau_{\mathrm{pause}}$ = 5 $\mu$s means that the annealing is paused at $t=$ 0.46 $\times$ 5 = 2.3 $\mu$s and lasts for 5 $\mu$s. The remaining 2.7 $\mu$s is scheduled for the latter part of the annealing. As a reference, we also performed continuous annealing without pause, designated $s_{\mathrm{pause}}=0.0$ in the results. The quantum-signature Hamiltonian can be mapped onto the unit Chimera graph of the D-Wave 2000Q device. In this way, more than 250 copies are mapped and annealed simultaneously. Each job submitted to the device comprised 1000 individual runs of the given Hamiltonian. The runs were averaged to give the final results. To remove the experimental biases, we repeated the measurements in two ways: First, we mapped different patterns in the physical qubits to reduce the pattern bias. For each unit Chimera graph, we selected ten mapping patterns from 144 (= 4! $\times$ 4! / 4) patterns implementing the quantum-signature Hamiltonian. Second, we ran ten experiments, randomizing the job order by shuffling in each run, to cancel out the general biases depending on the order of experiments, such as the unexpected drift of the operation temperature, and other time-dependent biases. By the above process, we generated 4\times(49+1)$$\times 250$$\times 1000$$\times 10$$\times 10=$$5\times 10^{9} data points. The quantum processing time and communication time of the whole process was approximately one hour.

Figure 4 plots the $P_{s}/P_{c}$ ratio versus anneal pause $s_{\mathrm{pause}}$ for various anneal time and pause periods $\tau_{\mathrm{anneal}}+\tau_{\mathrm{pause}}$. During longer annealing time, more thermal effects accumulate through the system–environment coupling; hence, increasing the annealing time increases the ratio. In addition, the ratio increased in the 0.3–0.7 range of $s_{\mathrm{pause}}$ under all four anneal-time conditions. The thermal fluctuations drive the system more efficiently in this region than outside the region. Although the isolated ground state is suppressed by quantum dynamics, it recovers during the anneal pause. This mitigation of the imbalance between $P_{s}$ and $P_{c}$ is favorable for fair sampling. The mechanism of the $s_{\mathrm{pause}}$-dependent increase of the ratio will be investigated in the next section.

III Numerical simulations

In our previous paper, we developed a method that represents open quantum systems by interpolating between quantum mechanics and thermodynamics. Energy in the system dissipates not through an external thermal bath, but via the incorporated dynamics governed by the master equation. We reported the dynamics of a single-spin system and larger systems by analytical and numerical expressions, respectively. In the latter approach, we did not take a continuous-time limit to obtain final merged differential equations. In the present paper, we expand our method to the analytical expression of larger systems. The nonlinear differential equations of the quantum-signature Hamiltonian are symbolically obtained using Mathematica (version 11.3).

The RHS of Eq. (17) is obtained by the following steps: (1) Construct the Schrödinger and master equations based on the given Hamiltonian in the $z$-basis. (2) Using Eq. (10), calculate the total derivative of the interpolated wave function from the total derivatives of $c_{i}(t)$ and $P_{i}(t)$. (3) Extract the first order coefficients of $dt$ by series expansion. The coefficients of $dt$ are the RHSs of the differential equations of the system. Appendix B provides the Mathematica code which generates C code of RHS of the dynamics, $\mathcal{F}(\cdots)$ in Eq. (17). The final RHS expression constitutes $2^{n}$ complex variables: $n+1$ variables in each equation including own state and $n$ neighbor states connected by a single-spin flip. The system size $n$ is eight in this study. So that, we have 256 equations and each equation has nine terms in RHS.

As the energy scale of the Hamiltonian is normalized ($|J_{ij}|\sim 1$ and $|h_{i}|\sim 1$), the simulation parameters to be specified are the coupling constant (interpolation parameter) $\alpha$ and the temperature $T$. To match the operation temperature of the D-Wave annealer (approximately 15 mKKing et al. ), we chose $T=0.3$ GHz ($\sim$ 14.4 mK). We conducted grid search of the coupling constant $\alpha$ to fit approximately the experimental peak ratio 0.0235 at $\tau_{\mathrm{anneal}}$ = $\tau_{\mathrm{pause}}$ = 1 $\mu$s to the experimental results. This process was started from the grid search region [0.001, 0.04] followed by [0.003, 0.005]. Finally, we chose $\alpha$ = 0.0045 and the peak ratio is 0.0238, while the ratio is 0.0196 and 0.0290 for $\alpha$ = 0.004 and 0.005. Although $s_{\mathrm{pause}}$ values for the peak of curves in the experiment and the numerical simulation does not match in our analysis, we prioritize to fit the peak ratio. The effects of varying $T$ and $\alpha$ around the chosen values are shown in Fig. 5.

The 256 equations were solved by the ordinary differential equation (ODE) solver scipy.integrate.ode in SciPy. The Adams method in the VODE solver was employed with nsteps = 200000 and atol = rtol = 1.136871e-13. The ODE time step was selected as 0.01. Smaller values of the time step (0.005, 0.002 and 0.001) did not change the results. In each time step, the total length of the state vector was normalized to 1 to prevent accumulation of systematic drift. The anneal scheduling parameters $A(s)$ and $B(s)$ of the D-Wave 2000Q applied to the numerical simulations were provided by D-Wave Systems Inc.

The simulation result (Fig. 6) reproduces the two fundamental characteristics in the experimental plot of Fig. 4; i.e., the anneal time and transverse field dependences of $P_{s}/P_{c}$. The probability range was larger in the numerical simulation than in the experiments, and the detailed curve shapes slightly differed between the two sets of results. To understand the phenomena that $P_{s}/P_{c}$ increases in a certain region (0.34–0.6 for $\tau_{\mathrm{anneal}}$ + $\tau_{\mathrm{pause}}=$ 1 + 1 $\mu$s), we investigated the energy-gap structure. The energy gaps are calculated from the instantaneous eigenstates and eigenvalues. If the system evolves adiabatically, it remains in the eigenstates with the lowest eigenvalue. The instantaneous energy-gaps for the lowest 20 excited eigenstates as functions of $s$ are plotted in Fig. 7. The instantaneous first excited eigenstate (blue dashed curve) around $s$ = 0.34 converges with three additional (second, third, and fourth) excited eigenstates. The thermal energy level of the system (0.3 GHz, dotted horizontal line in red) exceeds the energy of the first excited eigenstate near $s$ = 0.60. Above this value of $s$, change in the energy-gap structure is not dominant because the thermal fluctuation level is higher than the energy gap. In other words, the system freezes out, and the probabilities do not change during the anneal pause periodAmin (2015).

As the anneal control parameter $\tau_{\mathrm{anneal}}$ + $\tau_{\mathrm{pause}}$ increased, the simulated system exhibited more thermodynamic behavior, i.e. larger $P_{s}/P_{c}$ was observed as in Fig. 6. Static analysis, such as energy-gap structure analysis, is of limited help in understanding such behavior. Although additional static analyses can provide further information, they cannot inform beyond the dynamical analysis, as discussed later.

During the anneal pause, the quantum dynamics with a time-independent Hamiltonian do not essentially drive the system. Therefore, the probability should be transferred from the clustered ground states to the isolated state by the thermodynamic effects dominantly, which increases $P_{s}/P_{c}$. We investigate the dynamics in instantaneous eigenstates (Fig 8) to narrow down excited states involved in the dynamics, then focus on $z$-basis states (Fig. 9) to understand probability transfer among the isolated and clustered ground states.

Figure 8(a) shows how the simulated wave-function changes during the anneal pause at $s_{\mathrm{pause}}$ = 0.46 in terms of overlaps against the instantaneous eigenstates. In Fig. 8(a), changes in probability defined as difference of squared overlaps $|\braket{a}{\psi(t_{2})}|^{2}-|\braket{a}{\psi(t_{1})}|^{2}$ ($\{\ket{a}\}$: instantaneous eigenstates, $t_{1}$ and $t_{2}$: time of pause start and end) are plotted. The probability of the instantaneous ground state (the eigenstate with the lowest eigenvalue) decreased whereas those of the instantaneous fifth and 17th excited eigenstates increased. The same phenomena, probability change in these eigenstates, are commonly observed in the region 0.34–0.6, where $P_{s}/P_{c}$ ratio increased (see Fig. 8(b)). This observation implies these eigenstates have an essential role in the anneal pause period. In addition, at the end of the annealing, the fifth excited eigenstate was identified as the isolated state, as noted in an earlier study Boixo et al. (2013) whereas the 17th excited eigenstate was not well characterized.

To characterize the 17th excited eigenstate, the overlaps between the instantaneous 17th excited eigenstates at $s$ = 0.46 and the $z$-basis states in Fig. 9(a). This plot helps to identify essential spin configurations in $\sigma^{z}$ representation including the isolated and clustered ground states of $\mathcal{H}_{c}$. The $z$-basis states highlighted by different shapes and colors have relatively large overlaps, and they are potentially involved in the probability transfer among the isolated and clustered ground states. While Fig. 9(a) is a static analysis, Fig. 9(b) is a dynamic analysis and depicts the probability changes in the $z$-basis states during the anneal pause in the numerical simulation. Highlighted points in both figures from static and dynamic analyses are consistently non-zero and relatively large compared to the non-highlighted points (black dots). This consistency supports that these figures contain essential information to characterize the 17th excited eigenstate.

The $z$-basis states highlighted in Fig. 9 are located in the path from the clustered to isolated ground states in Fig. 10. The intermediate states can be grouped by the magnetization of their core spins. Groups with core-spin magnetizations of 2, 1, 0, -1, and -2 are labeled as CL, E1, E2, E3, and ISO, respectively.

The temporal probability behaviors of these states in the open ($\alpha$ = 0.0045) and closed ($\alpha$ = 0) quantum systems are displayed in panels (a) and (b) of Fig. 11, respectively. During the annealing period from $t$ = 0.46 to 1.46 $\mu$s, the wave function varied in the open system but remained static (except for the rotating phases) in the closed system. Referring to the single-spin flip chain shown in Fig. 10, the E1 group had a relatively large probability at the start of the anneal pause because its members were directly connected to the clustered ground states. The lower probabilities of groups E2 and E3 reflect their distances from the clustered ground state. During the first half of the anneal pause period, the slope of the probability versus time plot was larger in E3 than in the E2 and E1 groups, supporting a probability transfer toward the isolated ground state.

IV Summary and Discussion

Quantum annealers, such as the D-Wave machine, and other quantum devices are operated at the finite temperature and are coupled (albeit weakly) with the environment. By understanding such quantum devices as open quantum systems, we can expect to improve the control of systems.

The anneal pause function implemented in the D-Wave 2000Q can probe the thermodynamic effects in open quantum systems. We investigated how the anneal pause modulates the quantum-signature model. The $P_{s}/P_{c}$ ratio in the model reflects the balance between quantum mechanics and thermodynamics. From the experimental results in Fig. 4, we found dependences on both the total time of the dynamics and the strength of the transverse field during the anneal pause. The total time dependence can be explained by accumulative thermodynamic intervention from the environment such that long (slow) dynamics enhance the thermodynamic signature. The dependence on the strength of the transverse field during the anneal pause is nontrivial and is not well explained by the energy gap structure. To address these difficulties, we must magnify the phenomenon at the wave function scale in a quantum simulation analysis.

Performing a fine-grained analysis of the system, we employed a recently developed method for open quantum systemsKadowaki (2018). This method introduces energy dissipation by interpolating between the quantum dynamics and the thermodynamics. Although the method provides a numerical interpolation protocol, an accurate analysis requires an analytical form of the nonlinear differential equations that merge the whole process of the protocol. In the present paper, we established a method for deriving the nonlinear differential equations of multi-spin systems, supported by a symbolic computation program.

By properly selecting the parameters, temperature $T$ and coupling $\alpha$, in the dynamics of the quantum-signature model, we reproduced the two fundamental features observed in the experiment. The temperature was selected to approximate the operational temperature of the real quantum annealer whereas the coupling strength was fitted to the experimental results. Another potential parameter, which was not considered in this study, is the timescale difference between the quantum and thermodynamic phenomena. This parameter was omitted to reduce a large number of experiments and simulations when fitting multiple parameters simultaneously, but must be considered in a quantitative comparison.

From the simulation results, we could extract detailed information on the system dynamics. The analyses revealed the essential path driven by the 17th excited eigenstate, which transfers probability from the clustered ground states to the isolated ground state (the thermalization mechanism in this model). Owing to the many large overlaps between the excited eigenstate and the 11 states in the $z$-basis, a chain of single spin-flips establishes between the clustered and isolated ground states. This phenomenon informs us that thermalization in a quantum system can help to retrieve the classical ground states than pure quantum dynamics. Mitigating unfair sampling can improve statistical machine learning Korenkevych et al. .

For this purpose, we do not have an analytical method to select the best $s_{\mathrm{pause}}$ to control $P_{s}/P_{c}\sim 1$ at this point. Although we need to conduct a grid search on $s_{\mathrm{pause}}$, we found that the $P_{s}/P_{c}$ ratio can be improved more than one digit by combining longer annealing period and appropriate anneal pause setting (see the experimental results in Fig. 4). Further investigation can provide a better method to select $s_{\mathrm{pause}}$. In addition to the previous findings of thermal relaxation after anticrossingAmin et al. (2008); Venuti et al. (2017); Dickson et al. (2013), we found another role of thermal relaxation that transfers probabilities among ground states, then mitigates imbalance of probabilities. This is an example of how quantum and thermal effects work synergistically in quantum annealing.

For feasible open-quantum simulations, the method solves the nonlinear differential equations of $2^{n}$ complex variables originating from $n$-qubits without introducing the thermal bath variables. Relaxation of the system by the thermal bath is incorporated directly to the differential equations of the system through the interpolation with master equation parametrized by $\alpha$. This construct reduces the computational cost by removing the Hilbert space expansion for the thermal bath. Conversely, nonlinear differential equations are more computationally costly to construct and solve than linear ones. In addition, which systems can be resolved by this method is only partially answered in the present and the previous studies. We demonstrated the qualitative usefulness of this method by reproducing the trivial and nontrivial phenomena observed in open quantum systems. In the current formulation, the system–bath coupling is assumed as $\sigma^{z}$ interaction as our master equation is formulated by $\sigma^{z}$. Beyond this assumption, expansion of the method will be required.

In future work, the potential of the fair-sampling feature should be confirmed in experiments and theoretical analyses using open-quantum simulations. A former study of different models found that the classical ground state tends to accumulate more probabilities when the number of free spins is increasedMatsuda et al. (2009). The same conclusion was reached in our current analysis. Thus, we expect that fair sampling will emerge in the modelsMatsuda et al. (2009) and in more general cases. The biased sampling problem and possible performance improvements using thermal fluctuations have also been discussedMandrà et al. (2017). The authors suggested hybrid architectures that encourage thermal fluctuations, as also supported by our findings.

Our method is not only limited to simulate quantum annealers but also can be applied to enhance the understanding of quantum devices, such as noisy intermediate-scale quantum devices. Such further applications would reveal the usefulness and limitations of the method.

Acknowledgment

We are grateful to S. Tanaka, M. J. Miyama, H. Irie and S. Okada and D-Wave support team for their helpful discussions and suggestions. We are also grateful to the reviewer for pointing out the relation between our proposed method and stochastic Schrödinger equation. Part of the work is financially supported by MEXT KAKENHI Grant No. 15H03699 and 16H04382, and by ImPACT and JST START.

Appendix A Derivation for a two-level system

We demonstrate to derive differential equations of ID for a two-level system. The system consists of a spin interacting with the virtical field $\mathcal{H}_{c}=-h\sigma^{z}$ and the transverse field $\mathcal{H}_{q}=-\Gamma\sigma^{x}$. Thus, a matrix representation of the Hamiltonian is,

[TABLE]

Using the diagonal part of the Hamiltonian, we have the following transition rate matrix,

[TABLE]

where $\beta$ is the inverse temperature. For simplicity, we take a zero temperature limit, $\beta\to\infty$. Then, the matrix is

[TABLE]

A wave function of Schrödinger eqation and a probability distribution of master equation are two components, and we assume that they can be expressed using the same functions, $u(t)$ and $d(t)$,

[TABLE]

The following equations are obtained by substituting them for Eq. (14),

[TABLE]

By using series expansion for $dt$ we have difference equations,

[TABLE]

Finally differential equations of the two level system are derived as,

[TABLE]

One can confirm that these equations are equivalent to the density matrix representation previously reportedKadowaki (2018),

[TABLE]

Appendix B Generation of C code for simulation

This Mathematica code generates C code to calculate RHS of the interpolated dynamics for the two-level system in Appendix A. The eight-spin quantum-signature model can be calculated by defining $H$, $L$, $\psi$, and $P$ in the same manner.

[TABLE]