Improving quantum annealing of the ferromagnetic $ p $-spin model through pausing

TL;DR

Pausing during quantum annealing of the ferromagnetic p-spin model significantly enhances success probability by leveraging thermal relaxation, with optimal timing leading to up to 45% improvement and potential two-order magnitude gains.

Contribution

This study demonstrates that controlled pausing in quantum annealing can substantially improve performance on the ferromagnetic p-spin model, extending previous experimental findings.

Findings

Optimal pausing point is 60% beyond the avoided crossing time.

Pausing yields up to 45% increase in success probability.

Performance can improve by up to two orders of magnitude.

Abstract

The probability of success of quantum annealing can be improved significantly by pausing the annealer during its dynamics, exploiting thermal relaxation in a controlled fashion. In this paper, we investigate the effect of pausing the quantum annealing of the fully-connected ferromagnetic $ p $-spin model. We numerically show that (i) the optimal pausing point is 60% longer than the avoided crossing time for the analyzed instance, and (ii) at the optimal pausing point, we register a 45% improvement in the probability of success with respect to a quantum annealing with no pauses of the same duration. These results are in line with those observed experimentally for less connected models with the available quantum annealers. The observed improvement for the $ p $-spin model can be up to two orders of magnitude with respect to an isolated quantum dynamics of the same duration.

Equations

$$H_{\text{p}}=-\frac{n}{2}{\left(\frac{1}{n}\sum _{i=1}^n{\sigma }_i^z\right)}^p.$$

$$H_0=-\frac{1}{2}\sum _{i=1}^n{\sigma }_x.$$

$$H_Q(s)=A(s)H_0+B(s)H_{\text{p}},$$

$$H_{\text{S}}=-\frac{n}{2}|0⟩⟨0|.$$

$$V_{\text{QB}}=g\sum _{i=1}^n{\sigma }_i^z\otimes \sum _k\left(a_k+a_k^{†}\right),$$

$$\gamma (\omega )=2\text{uppi}\eta \frac{\omega e^{-\omega /{\omega }_{\text{c}}}}{1-e^{-\beta \omega }},$$

$$g^2\sum _k\delta (\omega -{\omega }_k)=\eta \omega e^{-\omega /{\omega }_{\text{c}}}.$$

$$\frac{d\rho(t)}{dt}=i\bigl{[}\rho(t),H_{Q}(t)+H_{\text{LS}}(t)\bigr{]}+\mathcal{D}\bigl{[}\rho(t)\bigr{]},$$

$$L_{ab}(t)=⟨{\epsilon }_a(t)⟩{\sum }_k{\sigma }_z^k{\epsilon }_b(t)|{\epsilon }_a(t)⟩⟨{\epsilon }_b(t)|,$$

$$L_0(t)=\sum _a⟨{\epsilon }_a(t)⟩{\sum }_k{\sigma }_z^k{\epsilon }_a(t)|{\epsilon }_a(t)⟩⟨{\epsilon }_a(t)|.$$

$$\mathcal{D}\bigl{[}\rho(t)\bigr{]}=\sum_{\alpha}\Bigl{[}C_{\alpha}(t)\rho(t)C_{\alpha}^{\dagger}(t)-\frac{1}{2}\left\{C^{\dagger}_{\alpha}(t)C_{\alpha}(t),\rho(t)\right\}\Bigr{]}.$$

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

$$+\sum _{\alpha }{C}_{\alpha }(t)\rho (t)C_{\alpha }^{†}(t),$$

$$H_{\text{eff}}(t)=H_Q(t)+H_{\text{LS}}(t)-\frac{i}{2}\sum _{\alpha }{C}_{\alpha }^{†}(t)C_{\alpha }(t).$$

$$|\tilde{\psi }(t)⟩=U(t,0)|\psi (0)⟩\approx \prod _{k\ge 0}{e}^{-i{H}_{\text{eff}}(k\delta t+\delta t/2)\delta t}|\psi (0)⟩,$$

$${\sigma }_{\text{MC}}^2(t)=\frac{1}{M(M-1)}\sum _{m=1}^M{\left[O_m(t)-\bar{O}(t)\right]}^2,$$

$$p_0$$

$$\equiv 1-⟨\tilde{\psi }(t+\text{Updelta}t)⟩\tilde{\psi }(t+\text{Updelta}t),$$

$$\rho (s_T)=\left(\begin{array}{cc}{\rho }_{11}(s_T)& 0\\ 0& {\rho }_{22}(s_T)\end{array}\right),$$

$$\frac{1}{\tau }\frac{d{\rho }_{11}(s)}{ds}={\Gamma }_{2\to 1}{\rho }_{22}(s)-{\Gamma }_{1\to 2}{\rho }_{11}(s),$$

$$\rho_{11}(s)=\rho_{11}(s_{T})-C\bigl{[}1-e^{-(s-s_{T})/s_{1}}\bigr{]},$$

$$\Phi (l_{\text{p}})={\Phi }_{\text{sat}}\left[1-\alpha e^{-(l_{\text{p}}-l_0)/T_{\text{r}}}\right],$$

$$H_{\text{LS}}=\sum _{a\ne b}S({\omega }_{ba})L_{ab}^{†}{L}_{ab}+S(0)\sum _{ab}{L}_{aa}^{†}{L}_{bb},$$

$$\displaystyle\begin{aligned} \mathcal{D}\bigl{[}\rho(t)\bigr{]}&=\sum_{a\neq b}\gamma(\omega_{ba})\Bigl{(}L_{ab}\rho(t)L^{\dagger}_{ab}-\frac{1}{2}\big{\{}L^{\dagger}_{ab}L_{ab},\rho(t)\big{\}}\Bigr{)}\\ &\quad{}+\gamma(0)\sum_{ab}\Bigl{(}L_{aa}\rho(t)L^{\dagger}_{bb}-\frac{1}{2}\left\{L^{\dagger}_{aa}L_{bb},\rho(t)\right\}\Bigr{)},\end{aligned}$$

$$S(\omega )=\mathcal{P}{\int }_{-\infty }^{+\infty }\frac{d{\omega }^{\prime }}{2\pi }\frac{\gamma ({\omega }^{\prime })}{\omega -{\omega }^{\prime }},$$

$$\beta {\omega }_{\text{c}}\gg 1;$$

$$\frac{1}{{\omega }_{\text{c}}\log (\beta {\omega }_{\text{c}})}<\min \left\{2{\tau }_B,\frac{{\tau }_{B}h}{\tau }\left(\frac{1}{{\Delta }^2}+\frac{{\tau }_B^2}{\tau }\right)\right\}.$$

$$\biggl{(}\sum_{i}\sigma^{z}_{i}\biggr{)}\ket{w}=(n-2w)\ket{w}.$$

$$H_{\text{p}}=-\frac{n}{2}\sum _{w=0}^n{(1-2w/n)}^p|w⟩⟨w|.$$

Full text

i) the optimal pausing point is $60\text{\,}\mathrm{\char 37\relax}$ longer than the avoided crossing time for the analyzed instance, and ii) at the optimal pausing point, we register a $45\text{\,}\mathrm{\char 37\relax}$ improvement in the probability of success with respect to a quantum annealing with no pauses of the same duration. These results are in line with those observed experimentally for less connected models with the available quantum annealers. The observed improvement for the $p$-spin model can be up to two orders of magnitude with respect to an isolated quantum dynamics of the same duration.

adiabatic quantum computation, quantum annealing, open quantum systems, $p$-spin model

I Introduction

Adiabatic quantum computation (AQC) Farhi et al. (2000); Albash and Lidar (2018a); Santoro et al. (2002) can be realized implementing quantum annealing (QA) algorithms Kadowaki and Nishimori (1998); Dickson et al. (2013) on quantum devices operating at low temperatures Harris et al. (2010, 2011). It is considered a heuristic technique for solving NP-hard optimization tasks that classical digital computers cannot solve efficiently. Defining and detecting a quantum speed up is a long-standing and delicate issue Rønnow et al. (2014); Amin (2015), however it is known that quantum annealing, in some cases, performs better than its classical counterpart, thermal annealing Santoro et al. (2002); Brooke et al. (1999); Albash and Lidar (2018b).

The main idea underlying AQC is to map the solution of a hard computational problem to the ground state of an Ising Hamiltonian Lucas (2014). The two states of Ising spins represent the two logical states of classical bits. At the starting time, the system is prepared in the ground state of a transverse field Hamiltonian, describing quantum spins fluctuating between their two classical values. These quantum fluctuations allow to explore the phase space through quantum tunneling. The Hamiltonian is slowly deformed to target the Ising Hamiltonian at the final time, hence the final ground state encodes the solution of the computational problem of interest. The adiabatic theorem of quantum mechanics ensures that a quantum system, prepared in its ground state at time $t=0$, will evolve remaining in its instantaneous ground state if the Hamiltonian is slowly varied in time, on a time scale proportional to the inverse of the minimal energy gap separating the ground state from the rest of the spectrum. In this framework, the time-to-solution is determined solely by the minimal energy gap.

However, in real annealing devices, mostly based on superconducting electronics Harris et al. (2010, 2011), such procedure runs in the presence of a finite-temperature environment, inducing dissipation. The current understanding of this problem is getting more and more detailed Amin (2015); Job and Lidar (2018); Marshall et al. (2017), and it is now clear that the time-to-solution is determined by a delicate balance between adiabaticity and thermal processes Amin et al. (2008); Mishra et al. (2018); Albash and Lidar (2015); Sarandy and Lidar (2005); Nalbach and Thorwart (2009); Javanbakht et al. (2015); Arceci et al. (2017); Smelyanskiy et al. (2017).

In general, the environment is expected to be detrimental for quantum annealing Arceci et al. (2018); Keck et al. (2017), despite the intrinsic robustness of AQC against decoherence Childs et al. (2001). However, sometimes the thermal environment can improve quantum annealing performances Smelyanskiy et al. (2017); Passarelli et al. (2018); Cangemi et al. (2018). This idea is interestingly pursued in a recent paper (Ref. Marshall et al. (2018)), in which the authors propose to pause the annealing at wisely chosen times, in order to take advantage also of thermal processes, in the search for the ground state of the final Hamiltonian. The possibility to pause the annealing dynamics is already implemented in the D-Wave 2000Q quantum computer dwa , and allowed the authors of Ref. Marshall et al. (2018) to experimentally test their hypothesis for instances of short range Ising models, easy to embed on the Chimera graph Choi (2008, 2011). However, the recently announced Pegasus graph Dattani et al. (2019), with increased connectivity, can pave the way for embedding also $p$-body all-to-all interacting models on the next generation of commercial quantum annealers, with less ancillary bits with respect to the Chimera graph Chancellor et al. (2016, 2017).

Indeed, computationally hard problems may require $p$-body interactions (with $p>2$), as for instance in Boolean satisfiability Cook (1971) or in the Grover search Jörg et al. (2010); Grover (1996). Grover-like search can be described by the so-called fully-connected ferromagnetic $p$-spin model (with large and odd $p$), whose embedding on the Chimera graph can be cumbersome and inefficient (in particular for large values of $p$), hence it is hardly studied with real devices. Nonetheless, this model has a remarkable relevance for many reasons. Originally introduced in Refs. Derrida (1981); Gross and Mezard (1984), the $p$-spin model is extensively studied in the context of spin-glasses Knysh (2016); Bapst et al. (2013); Bapst and Semerjian (2012) and quantum optimization Seoane and Nishimori (2012); Susa et al. (2018a); Ohkuwa et al. (2018). Its mean field character allows to recover analytic results in the thermodynamic limit for large $p$. Moreover, the existence of an exact ground state, combined with the non-triviality of its phase diagram, makes it a natural candidate for benchmarking the performances of quantum annealing in an exactly solvable, realistically hard case Wauters et al. (2017).

In this paper, we contribute to this topic investigating the quantum annealing of the fully-connected ferromagnetic $p$-spin model. In particular, we focus on the possible advantages coming from pausing the annealing dynamics, as proposed in Ref. Marshall et al. (2018) for Ising-type models. This paper is organized as follows. In Section II, we present the $p$-spin model Hamiltonian and address the computation of the dynamics in the presence of a realistic dissipation model. We adopt a Born-Markov approximation and use a quantum master equation in the Lindblad form, simulated using a Monte Carlo wave function method Yip et al. (2018). In Section III, we present our results, analyzing the effect of time and duration of the pause on the quantum annealing of a specific instance of the $p$-spin model. Finally, we present our conclusions in Section IV. We discuss the validity of the Born-Markov approximation in Appendix A.

II Model

II.1 Ferromagnetic $p$-spin model

In this work, we adopt natural units and fix $\hslash=k_{\text{B}}=1$. The dimensionless Hamiltonian of the $p$-spin model for $n$ qubits reads

[TABLE]

For $p$ odd, its ground state is ferromagnetic with all spins up-aligned. For $p$ even, the system is $Z_{2}$-symmetric, and there are two degenerate ground states. The limit $n\to\infty$, $p\to\infty$ (odd, with $p\leq n$) of this model encodes a Grover-like search Jörg et al. (2010); Grover (1996).

The Hamiltonian of Eq. (1) is reached via a quantum annealing starting from the dimensionless transverse field Hamiltonian $H_{\text{0}}$:

[TABLE]

We build the time-dependent Hamiltonian

[TABLE]

where $s=t/\tau\in[0,1]$ is the dimensionless time and $\tau$ is the annealing time. $A(s)$ and $B(s)$ encode typical annealing schedules Job and Lidar (2018); dwa , and are depicted in Fig. 1. $H_{Q}(s)$ is spherically symmetric, and both the ground state of $H_{0}$ and the ground state of $H_{\text{p}}$ belong to the symmetry subspace of maximum spin, having dimension $N=n+1$.

In order to obtain a quantum advantage for this adiabatic Grover-like search, the annealing schedule should be chosen carefully. Indeed, it is known that, using a linear schedule, the time to solution for search-like problems scales as $2^{n}$, which is the same as in classical approaches Roland and Cerf (2002). D-Wave’s annealing schedule is not linear, but still qualitatively different from the optimal annealing schedule for these classes of problems, reported in Refs. Roland and Cerf (2002); Morley et al. (2019) (where the same quadratic speedup of the Grover Grover (1996) algorithm is shown).

We fix the Hamiltonian parameters $n=20$ and $p=19$ such that the $p$-spin model Hamiltonian has a low lying spectrum very close to the one of the search Hamiltonian $H_{\text{S}}$:

[TABLE]

where $\ket{0}$ is the ferromagnetic ground state with all spins up-aligned. The comparison between the $p$-spin and the search Hamiltonian is shown in Fig. 2, where the eigenvalues of $H_{\text{p}}$ and $H_{\text{S}}$ are compared, as a function of the Hamming weight. The Hamming weight of a state is the number of $-1$ in its wave function, expressed in the single qubit basis $\set{1,-1}$, and is a measure of the distance of the given state from the ferromagnetic ground state (whose Hamming weight is of course zero). In Fig. 2, we compare the Hamming weights of the $p$-spin model with $p=3$, $p=7$, and $p=19$ with those of the search Hamiltonian of Eq. (4). The blue circles are for $p=3$, the red squares are for $p=7$, the green up triangles are for $p=19$, and the orange down triangles are for the search Hamiltonian of Eq. (4). Case $p=19$ and search show very similar spectral properties. In fact, the eigenvalues of the $p$-spin model with $p=19$ are nonzero only close to the target state (and to the very excited state with opposite magnetization). By contrast, $p=3$ and $p=7$ show qualitatively different behaviors.

In the thermodynamic limit $n\to\infty$, the $p$-spin model is subject to a dynamical quantum phase transition (QPT) Heyl (2018), at zero temperature, separating the para- and the ferromagnetic phases. The QPT occurs at $s=s_{\Delta}$, i. e., the time of the avoided crossing between the two lowest-lying energy levels. The QPT is second-order for $p=2$ and first-order for $p>2$ Bapst et al. (2013). The latter reduces the efficiency of adiabatic computations as the minimal gap $\Delta$ closes exponentially as a function of $n$ Bapst and Semerjian (2012). Several techniques are known to mitigate the detrimental effects of first-order QPTs in the quantum annealing of the $p$-spin model, such as non-stoquastic AQC Seki and Nishimori (2012); Seoane and Nishimori (2012); Nishimori and Takada (2017); Susa et al. (2017), inhomogeneous driving Susa et al. (2018b, a), or reverse annealing Chancellor (2017); Perdomo-Ortiz et al. (2011); Ohkuwa et al. (2018); dwa ; King et al. (2018). We will not discuss these techniques here.

II.2 Environment and dissipation

We model the environment as a collection of independent harmonic oscillators with Hamiltonian $H_{B}=\sum_{k}\omega_{k}a_{k}^{\dagger}a_{k}$, linearly coupled to the qubit system via the interaction potential

[TABLE]

where $g$ is the qubit-bath coupling energy and $a_{k}$ ($a_{k}^{\dagger}$) annihilates (creates) a boson in mode $k$ Caldeira and Leggett (1983); Leggett et al. (1987). We purposely use this operator as it does not break the spin symmetry, in contrast with more realistic dephasing models, where $g$ becomes $g_{ik}$. The annealer-environment Hamiltonian is $H(t)=H_{Q}(t)+H_{B}+V_{\text{QB}}$.

Noise in superconducting quantum annealers is modeled as the sum of a low frequency $1/f$ contribution and a high frequency Ohmic spectrum dwa ; Amin et al. (2009); Boixo et al. (2016); Massarotti et al. (2015), having the form

[TABLE]

where $\omega_{\text{c}}$ is a high frequency cutoff, $\beta=1/T$ and $\eta$ is a dimensionless coupling defined via

[TABLE]

The spectral function in Eq. (6) satisfies the quantum detailed balance (i. e., the Kubo-Martin-Schwinger) condition Albash et al. (2012). We will neglect the $1/f$ noise and fix $\omega_{\text{c}}=$1\text{,}\mathrm{THz}, $T=$12.1\text{\,}\mathrm{mK}$=$1.57\text{\,}\mathrm{GHz}, and $\eta=$1\text{\times}{10}^{-3}$$. With these parameters, the annealer-bath coupling strength is weak compared with the energy scales of the device, as shown in Appendix A. Thus, we can trace over bath states and adopt a Born-Markov approximation for the dynamics of the reduced system $\rho(t)$ Breuer and Petruccione (2007).

We employ a time-dependent Markovian quantum master equation (QME) in the Lindblad form to ensure complete positivity Lindblad (1976); Gorini et al. (1976); Albash et al. (2012). It reads

[TABLE]

where the explicit form of the Lamb shift $H_{\text{LS}}$ and the dissipator super-operator $\mathcal{D}$ is given in Appendix A. They are both expressed in terms of Lindblad operators,

[TABLE]

where $\ket{\varepsilon_{a}(t)}$ are the instantaneous eigenvectors of $H_{Q}(t)$. Each Lindblad operator $L_{ab}$ induces a quantum jump of frequency $\omega_{ba}(t)=\varepsilon_{b}(t)-\varepsilon_{a}(t)$, from $\ket{\varepsilon_{b}(t)}$ to $\ket{\varepsilon_{a}(t)}$, for $a\neq b$. Lindblad operators with $a=b$ induce dephasing.

II.3 Monte Carlo wave function

Time-dependent Lindblad equations can be reformulated as stochastic Schrödinger equations for the qubit wave function Yip et al. (2018). This approach is similar to the quantum jump method in quantum optics Mølmer et al. (1993). The wave function is evolved non-unitarily using an effective non-Hermitian Hamiltonian, which causes a loss of probability for the reduced system. The Lindblad operators make the system jump instantaneously between pairs of energy eigenstates, at random times and with rate $\gamma(\omega)$, thus restoring the lost probability. The QME is recovered by taking averages over a number $M$ of different stochastic realizations (or trajectories), with error scaling as $M^{-1/2}$. The advantage of this method, known as Monte Carlo wave function (MCWF), is twofold. On the one hand, we work with ket states rather than density matrices, hence our software is less memory demanding. On the other hand, the trajectories are independent, thus the technique is easily parallelizable. In what follows, we fix $M=5000$ trajectories as this number provides good accuracy in computations, e. g., a relative error $\mathop{}\!\delta{\left\langle O\right\rangle}/\left\langle O\right\rangle\sim$1\text{,}\mathrm{\char 37\relax}$$ on all measured observables $O$.

In order to unravel the QME through MCWF, Eq. (8) has to be rewritten in a more suitable form. We follow the algorithm described in Ref. Yip et al. (2018). We collect all dephasing operators into a single one:

[TABLE]

Then, neglecting accidental degeneracies, we denote each non-zero frequency $\omega_{ba}$ as $\omega_{\alpha}$, with $\alpha=1,\dots,n(n-1)$. We label the corresponding Lindblad operators $L_{ab}(t)$ (with $a\neq b$) using the same index $\alpha$. Defining also $\omega_{\alpha=0}=0$, we now have $n(n-1)+1$ Lindblad operators, each uniquely accompanied by its frequency $\omega_{\alpha}$, with $\alpha=0,\dots,n(n-1)$.

Defining $C_{\alpha}(t)=\sqrt{\gamma(\omega_{\alpha})}L_{\alpha}(t)$, the dissipator in Eq. (8) can be written compactly as

[TABLE]

We rewrite Eq. (8) as

[TABLE]

where the non-Hermitian Hamiltonian $H_{\text{eff}}(t)$ reads

[TABLE]

The non-Hermitian Hamiltonian of Eq. (13) generates a non-unitary dynamics. Given a starting state $\ket{\psi(0)}$, we discretize the time interval using a time step $\delta t$ and evolve the ket state using a first-order Trotter decomposition of the evolution operator,

[TABLE]

where states with a tilde are not normalized. In fact, the environment shifts the energy levels of the spectrum of the Hamiltonian $H_{Q}(t)$ and causes a decay of the norm of the state of the reduced system. To ensure convergence, we progressively reduce $\delta t$ until all analyzed observables are unaffected by the choice within Monte Carlo errors at all times. For an observable $O(t)$, the MCWF error $\sigma_{\text{MC}}(t)$ is defined as

[TABLE]

where $O_{m}(t)$ is the value assumed by $O(t)$ in the $m$th trajectory and $\bar{O}(t)$ is the mean value over all trajectories.

Our numerical MCWF is implemented in Fortran 90 and employs waiting time distributions and the integrated algorithm described in Ref. Yip et al. (2018), computationally faster than step-by-step stochastic evolutions. We report here the pseudo code for the single quantum trajectory; the reader can find more details in the original reference.

Draw a random number $r$ uniformly distributed in $[0,1]$. 2. 2.

Starting from a normalized state $\ket{\psi(0)}$, evolve non-unitarily using $H_{\text{eff}}(t)$ until $\braket{\tilde{\psi}(t^{*})}{\tilde{\psi}(t^{*})}=r$. At $t=t^{*}$, a quantum jump occurs. In fact, the norm squared of the unnormalized state $\ket{\tilde{\psi}(t+\mathop{}\!\Updelta t)}$ at a time $t+\mathop{}\!\Updelta t$ is the cumulative probability that no jumps have occurred during the time interval $[t,t+\mathop{}\!\Updelta t]$:

[TABLE]

where $\lambda(t)$ is the jump rate and $\mathop{}\!\Updelta t$ needs not to be infinitesimal. 3. 3.

Draw another random number $\mu$ uniformly in $[0,1]$ and select the quantum jump. The probability of jump $\alpha$ is given by $P_{\alpha}=\Pi_{\alpha}/\sum_{\alpha^{\prime}}\Pi_{\alpha^{\prime}}$, where $\Pi_{\alpha}=\braket{\tilde{\psi}(t^{*})}{C_{\alpha}^{\dagger}C_{\alpha}}{\tilde{\psi}(t^{*})}$. The index of the occurring jump is the smallest non-negative integer $m\leq n(n-1)$ satisfying $\sum_{\alpha=0}^{m}P_{\alpha}\geq\mu$. 4. 4.

Update the state as $\ket{\tilde{\psi}(t^{*}+0^{+})}=C_{m}\ket{\tilde{\psi}(t^{*})}$ and renormalize. Draw another random number $r$ and use $\ket{\psi(t^{*}+0^{+})}$ as the new normalized starting state. Repeat steps $14$ until the wanted annealing time.

II.4 Pausing the quantum annealing

A feature of the D-Wave 2000Q allows to pause the quantum annealing at user defined dimensionless time $s_{\text{p}}$ for a time duration $l_{\text{p}}$. During the time interval $l_{\text{p}}$, the system evolves with time-independent Hamiltonian $H_{Q}(s_{\text{p}})$, subject to dissipation. As discussed in Ref. Marshall et al. (2018), pausing has effect on an isolated quantum system only if the pause is inserted at the minimal gap. Instead, open quantum systems may benefit from pausing later on during the dynamics. The environment provides relaxation channels towards the ground state. For the analyzed instances, the authors found an enhancement in the fidelity $\Phi$, i. e., the ground state occupation probability at $s=1$, if $s_{\text{p}}$ is close to $s_{\Delta}$, with the maximum improvement being for $s_{\text{p}}$ approximately $10\text{\,}\mathrm{\char 37\relax}15\text{\,}\mathrm{\char 37\relax}$ of the total annealing time longer than $s_{\Delta}$. The authors interpret this result as a thermal-induced fidelity gain, concentrated around $s_{\Delta}$ as the relaxation rate is maximum in this region. Consistently with this interpretation, they register no perceptible effects far enough from the minimal gap, either because the tunneling amplitude is large (e. g., at early times) or because the dynamics is frozen out (e. g., at long times) Amin (2015). Our aim is to study the same effect for the $p$-spin model, using the QME in Eq. (8), unraveled using MCWF.

We present our results in the next Section. In the literature, the dissipative dynamics of the $p$-spin model has been studied for systems up to $n=16$ qubits and $p\leq 7$ Passarelli et al. (2018, 2019); Nakamura and Tanimura (2018). Here, we study a system of $n=20$ qubits, with $p=19$. With our parameters, the minimal gap is $\Delta\approx$0.14\text{,}\mathrm{GHz}$$, one order of magnitude less than the working temperature of the device. The avoided crossing between the two lowest energy levels occurs at $s_{\Delta}\approx 0.334$. In Fig. 3, we report the first $L=10$ eigenvalues and show the minimal gap $\Delta$, obtained via numerical diagonalization of $H_{Q}(s)$.

We choose the annealing time $\tau=$100\text{,}\mathrm{ns}. This is one order of magnitude smaller than the minimal annealing time of currently available quantum annealers, and might be experimentally accessible with next-generation D-Wave devices Job and Lidar ([2018](#bib.bib13)). In the unitary limit, this choice of $\tau$ makes the dynamics non-adiabatic, with a fidelity $\Phi\approx$5.51\text{\times}{10}^{-3}. Of course, one may decide to choose a longer annealing time to obtain a higher fidelity. For instance, with $\tau^{\prime}=10\tau$ one obtains $\Phi\approx$5.39\text{\times}{10}^{-2}$$, i. e., an improvement of one order of magnitude. However, in Refs. Passarelli et al. (2018, 2019), the authors have shown that in this non-adiabatic regime the low temperature environment can improve the annealing performances. Moreover, pausing the system can further increase the fidelity within the same time scales, as shown in the next Section.

III Results

III.1 Quantum annealing with no pauses

The annealing time $\tau$ is smaller than $h/\Delta^{2}$, where $h=\max_{s,a,b}\braket{\varepsilon_{a}(s)}{\partial_{s}H_{Q}(s)}{\varepsilon_{b}(s)}\approx$1.2\text{,}\mathrm{THz}, hence the quantum annealing is not adiabatic Albash *et al.* ([2012](#bib.bib63)). Landau-Zener diabatic transitions Landau ([1932](#bib.bib70)); Zener ([1932](#bib.bib71)) excite the qubit system and abruptly reduce the fidelity at the avoided crossing between the ground state and the first excited state. In Fig. [4](#S3.F4), we show the ground state population $\rho_{11}(s)$ in the energy eigenbasis as a function of time, during a dynamics driven by the Hamiltonian in Eq. ([3](#S2.E3)). The red line with squares represents the unitary dynamics ($\eta=0$) and the blue line with circles the dissipative one ($\eta=$1\text{\times}{10}^{-3}).

In the unitary case, the ground state population drops almost to zero after the avoided crossing at $s=s_{\Delta}$. Far from $s_{\Delta}$, the dynamics is frozen due to the large level spacing and the population remains constant. The fidelity of the algorithm is $\Phi\approx$5.51\text{\times}{10}^{-3}$$.

In the dissipative case, after $s=s_{\Delta}$, the bath-induced relaxation partially compensates diabatic transitions and improves the fidelity ($\Phi\approx 0.799$) with respect to the isolated case. Instead, bath-induced excitations right before the gap account for the small decrease in the plateau that can be seen for $s\lesssim s_{\Delta}$. Far from $s=s_{\Delta}$, thermal processes are exponentially suppressed and the ground state occupation probability is constant, in analogy with the isolated dynamics Marshall et al. (2018).

In what follows, we adopt a two-level model to give a simple description of the small decrease in $\rho_{11}(s)$ before $s_{\Delta}$, highlighted in Fig. 4. Let us consider the ground state and the first excited state, and focus on the region $s\in[s_{T},s_{\Delta}]$, where $s_{T}$ is the time before $s_{\Delta}$ where the instantaneous gap is equal to the temperature $T$. At shorter times, the thermal processes have scarce influence on the dynamics. At $s=s_{T}$, the system is almost completely in its ground state. As our simplifying hypothesis, we suppose that for $s\in[s_{T},s_{\Delta}]$ the gap is constant and equal to $\Delta$. In the energy eigenbasis, the density matrix of the qubit system at $s=s_{T}$ is diagonal as $B(s\leq s_{T})\ll A(s\leq s_{T})$, and can be written as

[TABLE]

where $\rho_{11}(s_{T})\approx 1$ and $\rho_{22}(s_{T})=1-\rho_{11}(s_{T})\approx 0$. The population transfer due to thermal processes can be effectively modeled by a classical master equation of the form

[TABLE]

where $\Gamma_{1\to 2}=\gamma(-\Delta)=e^{-\beta\Delta}\gamma(\Delta)$ and $\Gamma_{2\to 1}=\gamma(\Delta)$. The solution to this ME reads

[TABLE]

where $C=\rho_{11}(s_{T})-\tau s_{1}\gamma(\Delta)$ and s_{1}=1/\tau\sqrt{\gamma(\Delta)}\bigl{(}1+e^{-\beta\Delta}\bigr{)}. The ground state population at $s\to s_{\Delta}^{-}$ is then $\rho_{11}(s_{\Delta}^{-})\approx 0.975$, in agreement with numerical simulations.

III.2 Quantum annealing with a pause

Thermal effects can be enhanced by pausing the quantum annealing. We tested various pause lengths $l_{\text{p}}$, from $100\text{\,}\mathrm{ns}900\text{\,}\mathrm{ns}$. The maximum pause length $l_{\text{p}}$ gives a total annealing time $\tau^{\prime}=10\tau$, whose corresponding fidelity, in the absence of pauses, is $\Phi\approx$5.39\text{\times}{10}^{-2} for the unitary dynamics and $\Phi\approx 0.664$ for the dissipative one, with $\eta=$1\text{\times}{10}^{-3}. The effects of pausing substantially increases this fidelity. These results are $l_{\text{p}}$-independent, although there are small quantitative differences. In Fig. 5a, we show the fidelity $\Phi$ as a function of the pause time $s_{\text{p}}\in[0,1]$, for different values of the pause duration, l_{\text{p}}=\text{100\text{,}\mathrm{ns}400\text{,}\mathrm{ns}900\text{,}\mathrm{ns}}. Fig. 5b is a zoom in the region around the optimal pausing time $s_{\text{opt}}\approx 0.55$. Four different regions can be distinguished for all $l_{\text{p}}$.

When $s_{\text{p}}<s_{\Delta}$, the fidelity is not affected by pausing. Here, the tunneling amplitude is large compared with the thermal relaxation rate, and the system evolves quantum mechanically with little influence from the environment. 2. 2.

When $s_{\text{p}}\approx s_{\Delta}$, thermal processes are more frequent and, correspondingly, the relaxation rate is maximum. In this region, $\exp(-\beta\Delta)\sim 1$, hence the excitation rate $\Gamma_{1\to 2}$ is comparable with the decay rate $\Gamma_{2\to 1}$. When $s_{\text{p}}\lesssim s_{\Delta}$, most of the population is in the adiabatic ground state. Thus, transitions from the ground state to the first excited state are more probable than reverse processes. This imbalance causes a decrease in the observed fidelity after the pause, and the worst pausing time is sharply localized right before $s_{\Delta}$. The situation is reversed for $s_{\text{p}}\gtrsim s_{\Delta}$, and here the fidelity is slightly enhanced. This effect is more pronounced for longer $l_{\text{p}}$. 3. 3.

For $s_{\text{p}}=s_{\text{opt}}$, we observe a peak in the success probability for any $l_{\text{p}}>$100\text{,}\mathrm{ns}$$. The peak height increases with increasing $l_{\text{p}}$, following a saturation law of the form

[TABLE]

with fitted parameters $\Phi_{\text{sat}}=0.976\pm 0.007$, $\alpha=0.160\pm 0.005$, and $T_{\text{r}}=$4.1\text{,}\mathrm{ns}$\pm$0.4\text{,}\mathrm{ns}, and $l_{0}$ fixed to $l_{0}=$100\text{\,}\mathrm{ns}. In particular, $T_{\text{r}}$ is related to the thermal relaxation time of the many-body system, and $\Phi_{\text{sat}}$ is an estimate of the maximum fidelity that can be achieved by pausing the dynamics at the optimal point. The fidelity shows a peak almost at the time $s_{\text{p}}=s_{\text{opt}}$ independently of the pause duration $l_{\text{p}}$. By contrast, $l_{\text{p}}$ influences the time at which the fidelity goes back to its baseline value. In particular, for $l_{\text{p}}=$900\text{,}\mathrm{ns}$$, we register a $20\text{\,}\mathrm{\char 37\relax}$ increase of the fidelity with respect to the dissipative dynamics with no pause and total annealing time $\tau$, and a $45\text{\,}\mathrm{\char 37\relax}$ increase with respect to a dissipative dynamics of total annealing time $\tau^{\prime}$. The optimal pausing point $s_{\text{opt}}$ is independent of the pause length. 4. 4.

When $s_{\text{p}}\gg s_{\Delta}$, the fidelity is not influenced much by pausing. The ground state is well-separated in energy from the other levels (see Fig. 3). The eigenstates of the Hamiltonian are almost diagonal in the $\sigma_{z}$ basis and the qubit-bath coupling operator has exponentially small off-diagonal matrix elements. Perturbation theory predicts that thermal processes are exponentially suppressed, and the dynamics is frozen.

All the results are summarized in Fig. 6, where we show the fidelity as a function both of pause time and pause duration. Fig. 6b is a zoom in the region around $s_{\Delta}$. The heat maps show even more evidently that around an optimal pausing time $s_{\text{opt}}$ the fidelity abruptly increases, almost independently of $l_{\text{p}}$. This is evident in Fig. 6, where the dark shadow shows up at $s_{\text{p}}\approx s_{\text{opt}}$ and $l_{\text{p}}>$100\text{,}\mathrm{ns}. At shorter pausing lengths $l_{\text{p}}<$100\text{\,}\mathrm{ns}, this phenomenon is no longer visible. The white vertical line around $s_{\text{p}}\approx s_{\Delta}$ reflects a sharp decrease of the fidelity and can be better visualized in Fig. 6b. For the fully-connected $p$-spin model, the largest fidelity enhancement occurs for $s_{\text{opt}}\approx 1.65\,s_{\Delta}$, or, equivalently, $s_{\text{opt}}\approx s_{\Delta}+0.22$, while the worst pausing point is $s_{\text{wor}}\approx s_{\Delta}$.

Fig. 7 shows the differences in the dynamics of $\rho_{11}(s)$ when a pause of $l_{\text{p}}=$500\text{,}\mathrm{ns} is inserted at $s_{\text{opt}}$ (blue curve with a circle) or $s_{\text{wor}}$ (red curve with a square). When $s_{\text{p}}=s_{\text{opt}}$, the ground state population grows monotonically during the pause because of thermal relaxation. In fact, thermal excitations outside the ground state are suppressed as the spectral gap is large compared to $T$, e. g., $\varepsilon_{2}(s_{\text{opt}})-\varepsilon_{1}(s_{\text{opt}})\approx$100\text{\,}\mathrm{GHz}. On the other hand, for $s_{\text{p}}=s_{\Delta}$, the system can be excited due to $\Delta\ll T$ and $\rho_{11}(s)$ is reduced with respect to the previous case. Excitations and decays are almost equally probable, and produce evident noisy oscillations of $\rho_{11}(s)$ around a stationary value $\bar{\rho}_{11}\approx 0.25$ during the pause.

We repeated the simulations of this Section using Hamiltonian (4) as target Hamiltonian. Results are presented in Appendix B.

IV Conclusions

Quantum annealing is designed for solving NP-hard optimization tasks faster than classical techniques. This class of problems can often be mapped onto Ising Hamiltonians, and solved heuristically using the available hardware. The last generation of D-Wave devices has made possible to fine tune the annealing schedule, through reverse annealing, quenches, qubit-specific annealing offsets, and pausing. These advanced strategies allow to study several interesting problems. Currently, one of the main limitations is the sparse connectivity of the Chimera graph, since in many cases a cost function involves long range and multiple-body interactions among qubits. These problems are hard to embed efficiently in modern quantum annealers, as they require a large number of ancilla qubits and minor or parity embedding Choi (2008, 2011); Chancellor et al. (2017). In the context of coherent Ising machines Yamamoto et al. (2017), there have already been several proposals of all-to-all architectures based on superconducting electronics Puri et al. (2017); Nigg et al. (2017) or quantum optical systems Inagaki et al. (2016a); McMahon et al. (2016); Inagaki et al. (2016b), that should allow to overcome these limitations. At present, numerical simulations are fundamental to capture the physical behavior of fully-connected systems, while these promising platforms are developed and refined.

In this paper, we studied the quantum annealing of the ferromagnetic $p$-spin model for $n=20$ and $p=19$. The embedding of this model requires a large number of ancillae, for instance, using the decomposition of $p$-body interactions into $2$-body interactions described in Ref. Leib et al. (2016), or adopting the perturbative scheme described in Ref. Dodds et al. (2018), with a number of auxiliary qubits equal to the number of logical qubits. Instead, due to the spin symmetry of this model, we are able to simulate its dissipative dynamics using MCWF. Using a realistic model of dissipation, we showed that at anti-adiabatic times quantum annealing can be improved in the presence of a low temperature environment. More importantly, the beneficial thermal effects can be enhanced by pausing the annealing around an optimal pausing point. In fact, we showed that a quantum annealing of total time $\tau^{\prime}=$1000\text{,}\mathrm{ns} is circa $45\text{\,}\mathrm{\char 37\relax}$ more effective in reaching the $p$-spin ground state if a pause of length $l_{\text{p}}=$900\text{\,}\mathrm{ns} is inserted at $s=s_{\text{opt}}$ during an annealing of duration $\tau=$100\text{,}\mathrm{ns}, compared with a regular quantum annealing of time $\tau^{\prime}$ and no pause ($\Phi\approx 0.957$ versus $\Phi\approx 0.664$). In the isolated case, the fidelity is two orders of magnitude smaller ($\Phi\approx$5.39\text{\times}{10}^{-2}). We would need an annealing time $\tau\approx$50\text{,}\mathrm{\SIUnitSymbolMicro s}$$ to achieve a fidelity similar to that at the optimal pausing point. These results are in qualitative agreement with those of Ref. Marshall et al. (2018), and indicate that exploiting the environment in a controlled fashion allows to substantially reduce the time-to-solution compared with standard quantum annealing.

V Acknowledgments

We acknowledge the CINECA award under the ISCRA initiative, for the availability of high performance computing resources and support.

Appendix A Validity of the Lindblad approximation

The time-dependent Lindblad quantum master equation can be derived under very general assumptions Albash et al. (2012). The $\text{qubit}+\text{bath}$ Hamiltonian reads $H(t)=H_{Q}(t)+H_{B}+V_{\text{QB}}$, where all terms are defined in Section II. The bosonic bath is at equilibrium at temperature $T$ and is described by the correlation function $\mathcal{B}(t)=\left\langle X(t)X(0)\right\rangle$, where $X=\sum_{k}(a_{k}+a_{k}^{\dagger})$ and $X(t)$ is in the interaction picture with respect to $H_{B}$. The modulus of $\mathcal{B}(t)$ decays as $\exp(-t/\tau_{B})$ at small times and as $(t/\tau_{M})^{-2}$ at longer times, with $\tau_{B}=\beta/2\uppi+O(\omega_{\text{c}}^{-1})$ and $\tau_{M}=\sqrt{2\beta/\omega_{\text{c}}}$ Albash et al. (2012).

In the weak coupling regime, the Born, Markov and rotating wave approximations lead to a quantum master equation of the form reported in Eq. (8). This is equivalent to disregarding correlations and memory effects, and to ensuring complete positivity. The Lamb shift Hamiltonian and the dissipator super-operator in Eq. (8) are expressed as follows:

[TABLE]

where

[TABLE]

$\mathcal{P}$ denotes the principal value, and we omitted the time dependence from Lindblad operators and frequencies for brevity. We report here the constraints that all parameters must satisfy for the approximation to be valid. More details and extended discussion can be found in Ref. Albash et al. (2012).

• •

The adiabatic condition requires $\tau\gg h/\Delta^{2}$, where $h=\max_{s,a,b}\braket{\varepsilon_{a}(s)}{\partial_{s}H_{Q}(s)}{\varepsilon_{b}(s)}$. However, extensive use of adiabatic QMEs out of the adiabatic regime show results in reasonable agreement with other approaches.

• •

The perturbative corrections to the unitary dynamics due to the dissipator super-operator $\mathcal{D}$ must be small compared with $\Delta$. This implies $g^{2}\tau_{B}\ll\Delta$.

• •

The Markov approximation requires $\lvert\mathcal{B}(t)\rvert$ to decay more rapidly than the typical relaxation time scale $1/g$. Thus, $g\tau_{B}\ll 1$.

• •

The changes of the instantaneous eigenbasis of $H_{Q}(t)$ must be small during a time scale of the order of $\tau_{B}$. This requires $\tau\gg h\tau_{B}^{2}$.

Additionally, for Ohmic dissipation, $\omega_{\text{c}}$ must satisfy the following constraints:

[TABLE]

Except for the adiabatic condition, which we purposely violated as discussed in Section II, all other conditions are satisfied by our parameters.

Appendix B Comparison with the search Hamiltonian

In Section II.1, we discussed the similarity between the $p$-spin system with $n=20$ and $p=19$, and the typical search Hamiltonian of Eq. (4) (see Fig. 2). Indeed, these two systems can be mapped onto each other in the limit of large and odd $p$. Here we will quickly review this mapping.

We work in the subspace with maximum total spin, i. e., $S=n/2$. Within this subspace we have $n+1$ wave functions identified by $\ket{w}\equiv\ket{n-2w}$, with $w=0,1,\dots,n$. They are eigenstates of the $z$-component of the total spin, $\sum_{i}\sigma_{i}^{z}$, satisfying

[TABLE]

In this basis, Hamiltonian (1) is represented as

[TABLE]

Its ground state is the state $\ket{0}$. In the limit $n\to\infty$ and $p\to\infty$ ($p\leq n$), Hamiltonians (1) and (4) have all the same eigenvalues for $w=0,1,\dots,n-1$, and differ only by their last eigenvalue. This state has very large energy and does not affect the outcome of the computation, with our parameters. In this sense, the Hamiltonian (1) is analogous to the Hamiltonian (4).

In order to highlight the similarities between the $p$-spin and the search Hamiltonian, we repeat part of the analysis presented in Section III, using the search Hamiltonian. In Fig. 8, we report the first $L=10$ eigenvalues of the Hamiltonian $H_{Q}(s)$, as a function of the dimensionless annealing time $s$. Comparing the eigenvalues of Fig. 8 with the ones of the $p$-spin model (Fig. 3), we note a very similar behavior except for two minor differences. The first one is that, differently from the $p$-spin case, in the search model all the excited energy levels at $s=1$ have zero energy, and there are no anti-crossing involving pairs of excited eigenstates. The second difference is that the time of the minimal gap is slightly shifted ($s_{\Delta}^{\text{search}}\approx 0.335$), and the minimal gap is slightly reduced ($\Delta\approx$0.12\text{,}\mathrm{GHz}$$). However, these differences do not affect the dynamics significantly.

Using an annealing time $\tau=$100\text{,}\mathrm{ns}$$, we investigated the effects of pausing also using the search Hamiltonian as target, during a dissipative dynamics with the same parameters as those of Section III. At $s=1$, we measured the fidelity as a function of the pausing point $s_{\text{p}}$ for several values of the pause duration $l_{\text{p}}$ (e. g., for l_{\text{p}}=\text{100\text{,}\mathrm{ns}400\text{,}\mathrm{ns}900\text{,}\mathrm{ns}}). Our results are reported in Fig. 9. Panel 9b is a focus around the optimal pausing point, which in this case is $s_{\text{opt}}^{\text{search}}\approx 0.70$.

Similarly to the $p$-spin case, we distinguish four different regions.

When $s_{\text{p}}<s_{\Delta}^{\text{search}}$, the level spacing is large and the dynamics is unaffected by pauses. The baseline fidelity is slightly larger than the $p$-spin one ($\Phi\approx 0.849$ versus $\Phi\approx 0.799$). 2. 2.

When $s_{\text{p}}\approx s_{\Delta}^{\text{search}}$, the mean level spacing is one order of magnitude smaller than the temperature. Thermal processes are frequent and deplete ($s_{\text{p}}\lesssim s_{\Delta}^{\text{search}}$) or repopulate ($s_{\text{p}}\gtrsim s_{\Delta}^{\text{search}}$) the adiabatic ground state. 3. 3.

When $s_{\text{p}}\approx s_{\text{opt}}^{\text{search}}$, the fidelity increases. At difference with the $p$-spin case, the behavior around the optimal pausing point is noisier and a secondary maximum appear for $l_{\text{p}}=$900\text{,}\mathrm{ns}$$. 4. 4.

For $s_{\text{p}}\approx 1$, the fidelity goes back to its baseline.

The qualitative behavior of these data is very similar to that of the $p$-spin model (see Fig. 5).

We also computed all the diagonal elements of the qubit density matrix at the end of the annealing, for both the search Hamiltonian of Eq. (4) and the $p$-spin Hamiltonian of Eq. (1), for the same starting state and all equal parameters. In Fig. 10, we show the populations of the first $L=6$ states, as a function of the Hamming distance from the ferromagnetic ground state, using a semilogarithmic scale. The left panel is for the $p$-spin Hamiltonian and the right panel is for the search Hamiltonian. In both panels, and for each Hamming weight, the blue bars on the left refer to a dissipative quantum annealing with $\tau=$100\text{,}\mathrm{ns} and no pauses, whereas the red bars refer to a quantum annealing of the same duration, and a pause of $l_{\text{p}}=$900\text{\,}\mathrm{ns} inserted at the optimal pausing point ($s_{\text{opt}}\approx 0.55$ and $s_{\text{opt}}^{\text{search}}\approx 0.70$, respectively). The black dashed line represents the (all equal) populations of the starting state, i. e., $P=1/N$, with $N=n+1$. The result shows striking similarities between the two cases. However, some minor differences can be found. In particular, in the case of the $p$-spin Hamiltonian the population of the first excited state, with Hamming weight $w=1$, is slightly larger with respect to the corresponding one using the search Hamiltonian. However, in both cases the final population of these states is at least one order of magnitude smaller than the ground state population. All other states are even less populated. Moreover, if a pause is inserted at the optimal pausing point, the $p$-spin results become even closer to the ones of the search task, as the population of the state with $w=1$ drops below $P=1/N$, and all other populations are about three orders of magnitude smaller than the ground state occupation probability.