% decade, $n_d= {\mathcal N}_T \ln{(10)}/\ln{(\gamma_M/\gamma_m)}$, as follows
% ${\mathcal A} = \overline{b^2} P_0/4 = \overline{b^2}  n_d/(2 \ln(10))$.
Recall the form of $S(ω)$ in Eq.~\eqref{eq:totalpowerspectrum}.
By finding the correct parameters like $b$, $γ_M$, and $γ_m$ to best fit with the experimental data, one can find out the spectrum $S(ω)$. This is what we do in this section. We do not know what the power spectrum of the IBM quantum computer is.

Another direction is, given the power spectrum, simulate the qubit free evolutions and evolutions under dynamical decoupling. In this case, we know exactly what $S(ω)$ is in a given range of $ω$. With this information, one can deduce the parameters like $b$, $γ_M$, and $γ_m$, and simulate the free evolutions and evolutions under dynamical decoupling with such parameters, etc.
%The spectrum exhibits a crossover to $\omega^{-2}$ behavior at $\omega \sim \gamma_M$.
Fig.~\ref{fig:ibmqeo} illustrates both directions.
\begin{figure}[h!]
\centering
\includegraphics[width=17.5cm]{1f/big_picture_1f_Copy_3.pdf}
\caption{The processes of determining the spectrum~\cite{bylander2011noise} and performing simulations based on the spectrum.}
\label{fig:ibmqeo}
\end{figure}

% \subsection{Choice of circuit Hamiltonian and persistent operator}
% Although most use the qubit $2$-level model to perform the open system $1/f$ simulation, one can use the full Hamiltonian and persistent operators. Table.~\ref{tab:hf} lists all possible choices of Hamiltonian and fluctuations. Note that the logical $\overline{\sigma_z}$ is the same as  the previous section, i.e $\overline{\sigma_z} = \ket{\overline{0}}\bra{\overline{0}} - \ket{\overline{1}}\bra{\overline{1}}$. Note that the last two rows list the qubit $\text{ntruc}$-level models (with level of truncation $=\text{ntruc}$), which allow faster simulations.

% \begin{table}[h]
% \centering
% \begin{tabular}{ |c|c|c| }
%  \hline
%  & Hamiltonian & Fluctuation \\
%  \hline
%  Full & $H_{\text{circ}}$ &$I_{p}$ (Fluctuations on all energy levels)  \includegraphics[width=18mm]{1f/ipschematic.png}\\
%  \hline
%   Full & $H_{\text{circ}}$ &$\overline{\sigma_z}$   \\
%   \hline
%   $\text{ntruc} \times \text{ntruc}$ truncation & $[H_{\text{circ}}]_{\text{ntruc}\times \text{ntruc}}$ &$[I_{p}]_{\text{ntruc}\times \text{ntruc}}$   \\
%      \hline
%   $\text{ntruc} \times \text{ntruc}$ truncation & $[H_{\text{circ}}]_{\text{ntruc}\times \text{ntruc}}$ &$\sigma_z$\\
%   \hline
% \end{tabular}
% \caption{}
% \label{tab:hf}
% \end{table}

\subsection{Fitting results with IBM machine}
\label{sec:fitIBM}
We try to fit the IBM experimental data with our $1/f$ noise simulation. The goal is to reduce the sum of L2 norm of the difference of the simulation data and the experimental data, i.e. $∑_i s_i-e_i_2$, and extract the best fluctuator parameters. ($b, γ_min, γ_max$, etc)

We perform the simulations in a rotated frame, where the frame of reference is the original system Hamiltonian.
% D.D.
% \begin{figure}[h!]
%     \centering
%     \includegraphics[width=7.5cm]{1f/screendd.png}
%     \caption{XYXY sequence}
%     \label{fig:ibmqeo}
% \end{figure}

\subsubsection{Free induction decay: Fluctuation $A= I_p$}
We try to fit the simulations with the IBM experimental data of free induction decay, where the gates along evolution are the identity gates.

The fluctuation chosen in our simulations is $A= I_p$. We fix the standard  deviation of fluctuator strength to be $Δb /⟨b ⟩= 0.2$. The initial value of each fluctuator is randomly sampled according to the thermal equilibrium distribution $δp_eq= 0$ (i.e., each fluctuator is equally likely to be initialized in the state $1$ and state $-1$). The number of fluctuators is $20$.
We first output the best-fit parameters ($b, γ_min, γ_max$, etc) to the free induction decay curve with duration of $t_f = 5405$ns, with the initial state being $|+⟩$. The best fit parameters of fluctuators are $b = 0.12(9.2) ×10^-3, γ_min = 0.05 GHz, γ_max = 5 GHz$ (Fig.~\ref{fig:free_fit1}).
With these best-fit parameters, we can perform simulations with a new initial state. The simulation data of initial state $|1⟩$ also compares well with the free induction decay experimental data (Fig.~\ref{fig:free_fit2}). However, the simulation data of initial state $|0⟩$ does not compare well experiments (Fig.~\ref{fig:free_fit3}). If the initial state of the transmon qubit is the ground state, the experimental data suggests that it stays in the ground state, while the simulation data from $1/f$ noise suggests it goes to the maximally mixed state.

\begin{figure}[h!]
\centering
\subfigure[]{\includegraphics[width = 0.493\columnwidth]{1f/transmonibmplot_plus_0_Ip_2.pdf}\label{fig:free_fit1}}
\subfigure[]{\includegraphics[width = 0.493\columnwidth]{1f/transmonibmplot_one_0_Ip.pdf}\label{fig:free_fit2}}
\subfigure[]{\includegraphics[width = 0.493\columnwidth]{1f/transmonibmplot_zero_0_Ip.pdf}\label{fig:free_fit3}}
\caption{$t_f = 5405$ns, $\overline{b} = 0.12(9.2) \times 10^{-3}, \gamma_{\min} = 0.05~\text{GHz}, \gamma_{\max} = 5~\text{GHz}$, fluctuations in $I_p$, $20$ fluctuators, $200$ trajectories. Free induction decay. Rotated frame.}
\end{figure}

% \begin{figure}[h!]
%     \centering
%     \includegraphics[width=15.5cm]{1f/transmonibmplot_one_0_Ip.pdf}
%     \caption{$t_f = 5405$ns, $\overline{b} = 0.12(9.2) \times 10^{-3}, \gamma_{\min} = 0.05~\text{GHz}, \gamma_{\max} = 5~\text{GHz}$, fluctuations in $I_p$, $20$ fluctuators, $200$ trajectories. Free induction decay. Rotated frame.}
%     \label{fig:free_fit2}
% \end{figure}

% \begin{figure}[h!]
%     \centering
%     \includegraphics[width=15.5cm]{1f/transmonibmplot_zero_0_Ip.pdf}
%     \caption{$t_f = 5405$ns, $\overline{b} = 0.12(9.2) \times 10^{-3}, \gamma_{\min} = 0.05~\text{GHz}, \gamma_{\max} = 5~\text{GHz}$, fluctuations in $I_p$, $20$ fluctuators, $200$ trajectories. Free induction decay. Rotated frame.}
%     \label{fig:free_fit3}
% \end{figure}

% \clearpage
% \subsubsection{Free induction decay: Hamiltonian $=0$ , Fluctuation = $\overline{\sigma_z}$}
% The goal is to output parameters ($\overline{b}, \gamma_{\min}, \gamma_{\max}$, etc).
% \begin{figure}[h!]
%     \centering
%     \includegraphics[width=15.5cm]{1f/transmonibmplot_plus_0_sigmaz_2.pdf}
%     \caption{$t_f = 5405$ns, $\overline{b} = 0.043(9.2) \times 10^{-3}, \gamma_{\min} = 0.01~\text{GHz}, \gamma_{\max} = 1~\text{GHz}$, fluctuations in $\overline{\sigma_z}$, $20$ fluctuators, $200$ trajectories. Free induction decay.}
%     \label{fig:ibmqeo}
% \end{figure}

% \begin{figure}[h!]
%     \centering
%     \includegraphics[width=10.5cm]{1f/transmonibmplot_one_0_sigmaz.pdf}
%     \caption{$t_f = 5405$ns, $\overline{b} = 0.1(9.2) \times 10^{-3}, \gamma_{\min} = 0.05~\text{GHz}, \gamma_{\max} = 5~\text{GHz}$, fluctuations in $\overline{\sigma_z}$, $20$ fluctuators, $200$ trajectories. Free induction decay.}
%     \label{fig:ibmqeo}
% \end{figure}

% \begin{figure}[h!]
%     \centering
%     \includegraphics[width=10.5cm]{1f/transmonibmplot_zero_0_sigmaz.pdf}
%     \caption{$t_f = 5405$ns, $\overline{b} = 0.1(9.2) \times 10^{-3}, \gamma_{\min} = 0.05~\text{GHz}, \gamma_{\max} = 5~\text{GHz}$, fluctuations in $\overline{\sigma_z}$, $20$ fluctuators, $200$ trajectories. Free induction decay.}
%     \label{fig:ibmqeo}
% \end{figure}

%    \centering
%    \includegraphics[width=15.5cm]{transmonibmplot_1_0_Ip_2.pdf}
%    \caption{$t_f > 5405$ns, $\overline{b} = 0.001104, \gamma_{\min} = 0.05~\text{GHz}, \gamma_{\max} = 5~\text{Hz}$, fluctuations in $I_p$, $20$ fluctuators, $200$ trajectories.}
%    \label{fig:ibmqeo}

% \clearpage
% \subsubsection{Circuit Hamiltonian: on , Fluctuation = $I_p$}
% \begin{figure}[h!]
%     \centering
%     \includegraphics[width=15.5cm]{1f/transmonibmplot_plus_H_sigmaz_2_good_zoom.pdf}
%     \caption{$t_f =2560$ns, $\overline{b} = 0.143(9.2) \times 10^{-3}, \gamma_{\min} = 0.01~\text{GHz}, \gamma_{\max} = 1~\text{GHz}$, fluctuations in $I_p$, $20$ fluctuators, $200$ trajectories. Hamiltonian: on}
%     \label{fig:ibmqeo}
% \end{figure}

% \begin{figure}[h!]
%     \centering
%     \includegraphics[width=15.5cm]{1f/transmonibmplot_one_H_Ip.pdf}
%     \caption{$t_f = 2560$ns, $\overline{b} = 9.2 \times 10^{-2}, \gamma_{\min} = 100~\text{GHz}, \gamma_{\max} = 1000~\text{GHz}$, fluctuations in $I_p$, $20$ fluctuators, $200$ trajectories. Hamiltonian: on}
%     \label{fig:ibmqeo}
% \end{figure}

% \clearpage
% \subsubsection{Circuit Hamiltonian: on , Fluctuation = $\overline{\sigma_z}$}
% \begin{figure}[h!]
%     \centering
%     \includegraphics[width=15.5cm]{1f/transmonibmplot_plus_H_sigmaz_2.pdf}
%     \caption{$t_f =2560$ns, $\overline{b} = 0.143(9.2) \times 10^{-3}, \gamma_{\min} = 0.01~\text{GHz}, \gamma_{\max} = 1~\text{GHz}$, fluctuations in $\overline{\sigma_z}$, $20$ fluctuators, $200$ trajectories. Hamiltonian: on}
%     \label{fig:ibmqeo}
% \end{figure}

%    \centering
%    \includegraphics[width=15.5cm]{transmonibmplot_0_sigmaz.pdf}
%    \caption{$t_f > 5405$ns, $\overline{b} = 0.001104, \gamma_{\min} = 0.05~\text{GHz}, \gamma_{\max} = 5~\text{Hz}$, fluctuations in $\overline{\sigma_z}$, $20$ fluctuators, $200$ trajectories. Free induction decay.}
%    \label{fig:ibmqeo}

% \begin{figure}[h!]
% \centerline{\includegraphics[width=17.5cm]{1f/qctrlscreenshot2.png}}
% \caption{\tiny Q-ctrl Boulder-opal: webinar-introduction-to-robust-control\\
% URL: https://docs.q-ctrl.com/boulder-opal/application-notes/webinar-introduction-to-robust-control}
% \end{figure}

% The total Stochastic Hamiltonian $H(t)$ can be formulated as fluctuation of the  longitudinal noise or persistent current operator or other multiaxis noises :
% \begin{equation}
% \label{eq:stocH}
%     H(t) =  H_{S}(\phi_x) + \frac{1}{2} \,  \sum_{i}b_i \chi_i(t) A \,,
% \end{equation}
% where $\chi_i(t)$ switches with $\gamma_i/2$, with $\gamma_i$ sampled from $1/\gamma$ distribution in the range of $\{\gamma_{\text{min}}, \gamma_{\text{max}}\}$. $H_S$ is the time-independent, transmon circuit Hamiltonian; $I_p = I_p(\phi_x)$ is the persistent current operator.

% \begin{center}
% \begin{tabular}{ |c|c|c| c|}
%  \hline
%  & Longitudinal noise & Persistent current operator & Multiaxis noise \\
%  \hline
%  fluctuation (A) & $\overline{\sigma_z}$ &$I_{p}$ & \textbf{?}  \\
%  \hline
% \end{tabular}
% \end{center}
% $\overline{\sigma_z} = \ket{\overline{0}}\bra{\overline{0}} - \ket{\overline{1}}\bra{\overline{1}}$.

% \begin{figure}[h!]
% \centerline{\includegraphics[width=2.5cm]{1f/sigmaz.png}}
% \caption{$\overline{\sigma_z}$}
% \label{fig:1f}
% \end{figure}

% For fluctuator 1/f simulation, We do not perform any truncation and simulate the populations of all the $21$ levels , defined by the transmon circuit Hamiltonian.

% \subsection{Fitting experimental data and reverse engineering of the spectrum}

% There are two directions. One direction is reconstructing the $1/\omega^{\alpha}$ spectrum by fitting the experimental data of $T_1$ and $T_2$ times of free evolution and dynamical decoupling. This is what~\cite{bylander2011noise} did: Using the standard $1/f$ fluctuators simulation, find out the parameters that fit with the experimental data from free evolutions and dynamical decoupling, and reconstruct the $1/\omega^{\alpha}$ spectrum (and find out the value of $\alpha$). For more details, recall that a single two-state fluctuator is described by a stochastic variable $\chi (t)= \pm 1$ with rates $\gamma_{1\to 2}=\gamma_{2\to 1}=\gamma/2$. The power spectrum of the noise generated by the $i$-th fluctuator is:
% \begin{equation} \label{ns2}
% S_i(\omega)=b_i^2 \cL_{\gamma_i} (\omega)\, ,
% \end{equation}
% where
% %Consider a classical stochastic variable $\chi (t)$, which at any time takes the values $\chi (t) = \pm 1$.  It is thus suitable for describing
% %a two-state fluctuator that can find itself in one of two (meta)stable states, 1 and 2, and once a while it
% %makes a switch between them.
% %rates $\gamma_{1\to 2}$ and $\gamma_{2\to 1}$, which can be different.  For simplicity, we will limit ourselves to
% %symmetric RT process: $\gamma_{1\to 2}=\gamma_{2\to 1}=\gamma/2$.
% $\cL_{\gamma_i} (\omega)$ is a Lorentzian function of frequency,
% \begin{equation} \label{Lorentzian}
% \mathcal{L}_{\gamma_i} (\omega) \equiv \frac{1}{\pi} \frac{\gamma_i}{\omega^2 +\gamma_i^2}.
% \end{equation}
% The couplings $b_i$ are distributed with a  small dispersion around an average value $\overline b $. To obtain the $1/\omega$ low-frequency behavior of the power spectrum one has to assume that the distribution function $\cL_{\gamma}$ behaves as $1/\gamma$.
% Under these conditions the total power spectrum reads
% \begin{equation} \label{ns_narrow_b}
% S(\omega) = \overline{b^2} \int_{\gamma_m}^{\gamma_M} d\gamma \, \frac{P_0}{2 \gamma} \,  \cL_{\gamma} (\omega)
% \approx    \frac{{\mathcal A}}{\omega}\, ,
% \end{equation}

% where the amplitude $\mathcal A$ can be expressed in terms of the number of \fs per noise
% decade, $n_d= {\mathcal N}_T \ln{(10)}/\ln{(\gamma_M/\gamma_m)}$, as follows
% ${\mathcal A} = \overline{b^2} P_0/4 = \overline{b^2}  n_d/(2 \ln(10))$. Therefore, by finding the correct parameters like $\overline{b}$, $\gamma_M$, and $\gamma_m$, one can find out the spectrum $S(\omega)$.
% This is a direction one can perform on IBM machines, by fitting the simulations with the experimental decay curves.

% Another direction is, given the PSD spectrum, simulate the qubit free evolutions and evolutions under dynamical decoupling. In this case, we know exactly what $S(\omega)$ is in a given range of $\omega$. With this information, one can deduce the parameters like $\overline{b}$, $\gamma_M$, and $\gamma_m$, and simulate the free evolutions and evolutions under dynamical decoupling, etc. This is the direction one can perform on QEO flux qubit.
% %The spectrum exhibits a crossover to $\omega^{-2}$ behavior at $\omega \sim \gamma_M$.
% Fig.~\ref{fig:ibmqeo} listed out both directions.
% \begin{figure}[h!]
%     \centering
%     \includegraphics[width=18.5cm]{1f/big_picture_1f.pdf}
%     \caption{The processes of determining the spectrum~\cite{bylander2011noise} and performing simulations based on the spectrum.}
%     \label{fig:ibmqeo}
% \end{figure}

% \subsection{Choice of circuit Hamiltonian and persistent operator}
% Although most use the qubit $2$-level model to perform the open system $1/f$ simulation, one can use the full Hamiltonian and persistent operators. Table.~\ref{tab:hf} lists all possible choices of Hamiltonian and fluctuations. Note that the logical $\overline{\sigma_z}$ is the same as  the previous section, i.e $\overline{\sigma_z} = \ket{\overline{0}}\bra{\overline{0}} - \ket{\overline{1}}\bra{\overline{1}}$. Note that the last two rows list the qubit $\text{ntruc}$-level models (with level of truncation $=\text{ntruc}$), which allow faster simulations.

% \begin{table}[h]
% \centering
% \begin{tabular}{ |c|c|c| }
%  \hline
%  & Hamiltonian & Fluctuation \\
%  \hline
%  Full & $H_{\text{circ}}$ &$I_{p}$ (Fluctuations on all energy levels)  \includegraphics[width=18mm]{1f/ipschematic.png}\\
%  \hline
%   Full & $H_{\text{circ}}$ &$\overline{\sigma_z}$   \\
%   \hline
%   $\text{ntruc} \times \text{ntruc}$ truncation & $[H_{\text{circ}}]_{\text{ntruc}\times \text{ntruc}}$ &$[I_{p}]_{\text{ntruc}\times \text{ntruc}}$   \\
%      \hline
%   $\text{ntruc} \times \text{ntruc}$ truncation & $[H_{\text{circ}}]_{\text{ntruc}\times \text{ntruc}}$ &$\overline{\sigma_z}$\\
%   \hline
% \end{tabular}
% \caption{}
% \label{tab:hf}
% \end{table}

\subsubsection{Dynamical decouplings: Fluctuation $A= [I_{p}]_{2\times 2}$}
We try to fit the simulations with the IBM experimental data of dynamical decoupling.

We fix the standard deviation of fluctuator strength to be $Δb /⟨b ⟩= 0.2$. The initial value of each fluctuator is randomly sampled according to thermal equilibrium distribution $δp_eq= 0$. The number of fluctuators is $20$.
Based on the fitting with free evolution, we perform simulations with the $XYXY$ dynamical decoupling sequence.
We model the dynamical decoupling as ideal pulses in terms of Pauli two-level $σ_x$ and $σ_y$ matrices. Therefore, our system Hamiltonian and fluctuators are also truncated to the first two energy levels. $A = [I_p]_2×2$.

We first find out the best parameters in ($b, γ_min, γ_max$, etc) to the free induction decay curve with duration of $t_f = 15000$ns, with initial state $|+⟩$. The best fit parameters of fluctuators are $b = 0.12(9.2) ×10^-3 GHz, γ_min = 0.05 GHz, γ_max = 5 GHz$ (Fig.~\ref{fig:dd_fit1}).
With these best-fit parameters, we can perform simulations of dynamical decoupling (XYXY sequence between each data point) under $1/f$ noise, also with initial state $|+⟩$. However, as we can see in Fig.~\ref{fig:dd_fit2}, the distance between the experimental data and simulation data diverges at long times. Fast noises, which cannot be removed by dynamical decoupling, are present in Transmon qubit, and thus the experimental dynamical decoupling curve decays faster than pure $1/f$ noise simulation.

Nevertheless, we increase the average fluctuator strength by a bit, and can obtain a better fit (Fig.~\ref{fig:dd_fit3}).

\begin{figure}[h!]
\centering
\subfigure[]{\includegraphics[width = 0.493\columnwidth]{1f/newdataddornot3.pdf}\label{fig:dd_fit1}}
\subfigure[]{\includegraphics[width = 0.493\columnwidth]{1f/newdataddornot2_new.pdf}\label{fig:dd_fit2}}
\subfigure[]{\includegraphics[width = 0.493\columnwidth]{1f/newdataddornot2_better_new.pdf}\label{fig:dd_fit3}}
\caption{Fluctuations in $[I_{p}]_{2\times 2}$, $20$ fluctuators, $300$ trajectories. Rotated frame. $\gamma_{\min} = 0.01~\text{GHz}, \gamma_{\max} = 0.1~\text{GHz}$. (a), (b): $\overline{b} = 9.2929 \times 10^{-4}~\text{GHz}$, (c) Modified averaged fluctuator strength: $\overline{b} = 1.05 \times 10^{-3}~\text{GHz}$.}
\end{figure}

\section{Conclusion of this chapter}
\label{secconclude: 1f}
In this chapter, based on the previously studied Bloch vector-magnetic field approach, we have formulated a stochastic fluctuator Hamiltonian simulation method that includes a series of fluctuator terms with different flipping frequencies and Gaussian couplings that altogether constitute the $1/f$ power spectrum. We, for the first time, show how to incorporate them into the inherently time-dependent annealing Hamiltonians, and construct an overall time-dependent stochastic Hamiltonian to perform parallelizable simulations. We thoroughly analyzed how noise with $1/f$ power spectrum affects the quantum annealing, and we also studied  temperature effects, loss of coherence ($T_1$ and $T_2$) in QA under different fluctuator parameters. We developed a flexible approach, allowing the noise fluctuators to take any form of multi-level operators and any noise-axis direction depending on the superconducting qubits and circuit decoherence process. We showed how to append these operators to two different superconducting circuit Hamiltonians. We fit our simulation results with experimental data, and showed how dynamical decoupling can extend the coherence of qubits under $1/f$ noise.

\section{Acknowledgments of this chapter}
Sec.~\ref{sec:fitIBM} would be impossible without Vinay Tripathi's experimental data and discussions.
%Other topics such as crosstalks modelling in IBM quantum machines are neglected to avoid overlapping with other authors' contributions in a future published paper.
The codes in this chapter are available at:
\url{https://github.com/USCqserver/1fnoise}

% \begin{figure}[h!]
%     \centering
%     \includegraphics[width=15.5cm]{1f/newdataddornot3.pdf}
%     \caption{ $\overline{b} = 9.2929 \times 10^{-4}~\text{GHz}, \gamma_{\min} = 0.01~\text{GHz}, \gamma_{\max} = 0.1~\text{GHz}$, fluctuations in $[I_{p}]_{2\times 2}$, $20$ fluctuators, $300$ trajectories. Rotated frame.}
%     \label{fig:dd_fit1}
% \end{figure}

% \begin{figure}[h!]
%     \centering
%     \includegraphics[width=15.5cm]{1f/newdataddornot2_new.pdf}
%     \caption{ $\overline{b} = 9.2929 \times 10^{-4}~\text{GHz}, \gamma_{\min} = 0.01~\text{GHz}, \gamma_{\max} = 0.1~\text{GHz}$, fluctuations in $[I_{p}]_{2\times 2}$, $20$ fluctuators, $300$ trajectories. Rotated frame, with dd.}
%     \label{fig:dd_fit2}
% \end{figure}

% \begin{figure}[h!]
%     \centering
%     \includegraphics[width=15.5cm]{1f/newdataddornot2_better_new.pdf}
%     \caption{ $\overline{b} = 1.05 \times 10^{-3}~\text{GHz}, \gamma_{\min} = 0.01~\text{GHz}, \gamma_{\max} = 0.1~\text{GHz}$, fluctuations in $[I_{p}]_{2\times 2}$, $20$ fluctuators, $300$ trajectories. Rotated frame, with dd.}
%     \label{fig:dd_fit3}
% \end{figure}
\chapter{Quantum feedback control}
\label{chap: feedback}
\paragraph{}
As mentioned in the conclusion~\ref{secconclude: qt} of Chapter.~\ref{chap: qt}, an extension of the trajectories method is the study of weak measurement and feedback in annealing system. We devise a quantum feedback error correction method to reverse the effect (error) of thermal excitations during quantum annealing.

\section{Introduction}
In QA, the solution of the computational problem is encoded as the ground state of the final Hamiltonian. In an ideal situation, a quantum state is initiated as the ground state of an initial Hamiltonian, and is evolved to the ground state of the final Hamiltonian. It is desirable is to maintain the state during the anneal as the ground state of the time-dependent Hamiltonian. However, due to diabatic or thermal transitions, this is not always successful, i.e., the quantum state at the end of the process is not the ground state of the final Hamiltonian.

Quantum error correction (QEC) aims to correct such quantum errors. There were previous extensive studies on quantum error correction in QA~\cite{young2013error, lidar2008towards, mishra2016performance}, mostly centered on using error correcting codes to perform majority voting at the end of the anneal. This chapter is about a particular active quantum error correction approach that can used \emph{during} an anneal: feedback-based error correction for QA.

Quantum feedback-based error correction aims to manipulate the system's quantum state or trajectory to evolve the system towards some desired outcome. Depending on the delay between detection and feedback, the effect of the feedback can be Markovian or non-Markovian. We study the following topics about feedback-based error correction on QA.

\subsection{Feedback and control}
In active quantum error correction, one may want to know when and what types of corrections to apply to the system qubit(s) during the quantum operation/process. One way to achieve this goal is to perform measurement on the quantum system at a current time-step and use the measurement result as feedback to determine/improve the error correction for future time-steps.

%For the particular case of quantum annealing, the system qubits evolve according to a time-dependent system Hamiltonian $H_{S}$. Mathematically, $H_{S} = A(t) H_{i} + B(t)H_{f}$, where $H_{i}$ is the initial system Hamiltonian and $H_{f}$ is the final system Hamiltonian, and $A(t) = 0$ at the end of anneal $t = t_{f}$, and $B(t) = 0$ at $t = 0$ (For example, see Fig.~\ref{fig: DWavenasaschedule}). Decoherence arises due to the coupling to a quantized bath (environment). The coupling is described by another system-bath Hamiltonian $H_{SB}$.
%Certain decoherence models can be assumed in order to derive a continuous equation describing evolution of system qubits. As the qubit(s) is evolving, measurement and feedback correction can be performed on the qubit(s). One would need to modify the continuous equation to describe the combined effect.

% The definition of feedback is quite broad. The feedback is one such that one observe some information and recycle it into the operation such that the operation becomes better. One example of this is machine learning. The idea of quantum feedback error correction is as follows
\subsection{Weak measurement}
One way to obtain information about quantum system is by weak measurement. Weak measurement is a type of quantum measurement that results in an observer obtaining very little information about the system on average, but also disturbs the state very little~\cite{brun2002simple}. Weak measurement is also often referred to as continuous measurement, where (weak) measurement is performed continuously in successive tiny time-steps. Weak measurement arises from the system qubit successively interacting with environmental (probe) qubits which are subsequently measured (See Fig.~\ref{fig:weakmeasurement}). Depending on the basis in which the environmental qubits is measured, one can obtain different types (jumpy/diffusive) of measurement trajectories for the system qubit being measured.

\begin{figure}
\centering
\includegraphics[scale=0.5]{feed/continuousmeasurementcenter.png}
\caption{By successively interacting the system with environmental (probe) qubits which are subsequently measured, we can obtain the weak measurement trajectories.}
\label{fig:weakmeasurement}
\end{figure}

When little information is gathered during the measurement, the quantum state is almost entirely preserved. Nevertheless, taking a sufficiently large sample of data still permits the same average information to be extracted as a textbook projective measurement that fully collapses the state. By continuously monitoring partial information about a quantum system, one can control experimental parameters in real time for a variety of purposes, such as error correction, coherence stabilization, state purification.

\section{Thermal excitation as error: detection}
In QA, thermal excitation is considered an error since we want to preserve the ground state of the system Hamiltonian. To do that in feedback-based error correction, the system evolution is continuously monitored and feedback is applied conditioned on each weak measurement result. A photo-detection current $I(t)$ of thermal excitation is obtained, $0<t<t_f$.

With the assumption of weak-coupling limit, one can derive a master equation describing the evolution of system qubit as:
\begin{align}
\frac{d}{dt} \rho(t) =  -i \left[H_{\textrm{S}}(t) + H_{\textrm{LS}}(t), \rho(t) \right] &+ D[A_{-\omega}(t)]\rho(t) + D[A_{0}(t)]\rho(t) + D[A_{\omega}(t)]\rho(t) \,.
\end{align}

\begin{equation}
\label{eq:exciteop}
A_{-\omega}(t) \propto |\eps_1(t) \rangle  \langle \eps_0(t)|
\end{equation}
is the Lindblad operator of thermal excitation in the instantaneous energy eigenbasis.
$D[ · ]$ refers to the decoherence channel (excitation/dephasing/damping). This master equation can be unraveled into a stochastic Schr\"{o}dinger equation describing the weak measurement.

The system evolution under annealing is continuously monitored and feedback is applied conditioned on each continuous measurement result. The continuous measurement results are in the form of photodetection current $I(t)$, $0≤t ≤t_f$. The photodetection current $I(t)$ corresponds to the detection of excitation of the system qubit. $I(t) = 1$ when excitation is detected, $=0$ otherwise  (See Fig.~\ref{fig:photocurrent}).

Conditioned on this current, we can obtain an evolution (unravelling in the quantum jump picture) of the \textit{conditioned} state $ρ_I(t)$.

\begin{figure}[h!]
% \includegraphics[width=0.25\textwidth]{Itillustration.png}
\centering
\includegraphics[width=0.35\textwidth]{feed/plotI.pdf}
\caption{Photodetection current $I(t)$.}
\label{fig:photocurrent}
\end{figure}

\section{Delay between detection and feedback}
After the detection of a thermal excitation ($I(t) = 1$) from ground to excited state, we then want to apply feedback to correct such errors. However, there usually is a delay between detection and feedback. The most general expression for the effect of the feedback on the conditioned state $ρ_I(t)$ is~\cite{Wiseman-Milburn-book}:
\begin{equation}
[\dot{\rho}_{\text{I}}(t)]_{\text{fb}} = \mathcal{F}[t, \mathbf{I}_{[0, t)}]\rho_{\text{I}}(t)   \,,
\end{equation}
where $ 𝐈_[0, t)$ is the complete photocurrent record from time $0$ to $t$, and $ℱ[t, 𝐈_[0, t)]$ is a superoperator.
Assuming that the feedback is of linear functional form, the feedback evolution can be expressed as:
\begin{equation}
[\dot{\rho}_{\text{I}}(t)]_{\text{fb}} = \int_{0}^{\infty} h(s) I(t-s) \mathcal{L}_f \rho_{\text{I}}(t) ds\,,
\end{equation}
where $ℒ_f$ is the Liouville superoperator representing the feedback control.
If the delay is fixed as $τ$, an equation incorporating the effect of the feedback is:
\begin{equation}
\label{Eq:implicitfeedbackevolution}
[\dot{\rho}_{I}(t)]_{\text{fb}} = I(t-\tau)\mathcal{L}_f\rho_{I}(t) \,.
\end{equation}
Note that Eq.~\eqref{Eq:implicitfeedbackevolution} is non-Markovian for finite $τ$. An analytical solution of Eq.~\eqref{Eq:implicitfeedbackevolution} is hard to obtain and numerical simulation is needed.

\section{Feedback superoperator $\mathcal{L}_f$: Lindbladian}
The feedback superoperator $ℒ_f$ can take the form of a Lindbladian or a Hamiltonian. In the form of a Lindbladian:
\begin{equation}
\mathcal{L}_f[\,\cdot\,] = F(t)\,\cdot\,F^{\dagger}(t) - \frac{1}{2}\left\{F^{\dagger}(t)F(t), \,\cdot\, \right\} \,.
\end{equation}
Physically, for the purpose of QA, this Lindbladian represents a cooling/damping channel in which the qubit is relaxed from an excited state to the ground state.

$F(t)$ can take many forms. For example, $F(t) = A^†_-ω(t-τ) = | _0(t-τ)   ⟩⟨_1(t-τ)  |$, the complex conjugate of the thermal excitation operator [Eq.~\eqref{eq:exciteop}] a time $τ$ ago.
This exactly cancels the error caused by thermal excitation if $τ= 0$. With this form of $F(t)$, we can devise a quantum feedback algorithm -- see Algorithm.~\ref{alg:feedbackalgo}. We perform such trajectories feedback simulation in Section~\ref{sec:feedback case studies}.

\begin{algorithm}[H]
\caption{Quantum feedback algorithm with delay}\label{alg:feedbackalgo}
\begin{algorithmic}
\Procedure{Feedback with delay}{$t_f, t_j = 0, \tau$}\Comment{Initialize $t_j = 0$.}
%\State $t_j\gets a\bmod b$
\While{$t_j + \tau \leq t_f$}
\State Evolve according to quantum trajectory method.
%\State $r\gets a\bmod b$
\If {Excitation ${A}_{-\omega}(t)$ occurs at $t = t_j$ to $\ket{\eps(t_j)}$ (Indicated by $I(t) = 1$)}
\State Apply feedback $F(t_j+\tau) = A^{\dagger}_{-\omega}(t_j)$ at $t_j + \tau$ to $\ket{\eps(t_j +\tau)}$.
\EndIf
\EndWhile
%\State \textbf{return} $b$\Comment{The gcd is b}
\EndProcedure
\end{algorithmic}
\end{algorithm}

\section{Case studies}
\label{sec:feedback case studies}
We consider the case of a single qubit with annealing Hamiltonian, i.e.,
$H_S(t) = A(t) σ_x + B(t) σ_z$. The $A(t) = A(1-t/t_f)$ and $B(t) = B(t/t_f)$ schedules are the linear schedules, with $A = 1.57$ GHz and $B = 2.01$ GHz (in units $ħ= 1$).

We use $F(t) = A^†_-ω(t-τ) = | _0(t-τ)   ⟩⟨_1(t-τ)  |$, with different delays $0≤τ≤t_f$.
In Fig.~\ref{fig:feedbackplotlinearsimple}, we simulate the performance of the feedback error protocol for different delay $τ$ and plot the ground state populations along the anneal.

\begin{figure}[h!]
\centering
% \includegraphics[width=0.65\textwidth]{feedbackplotlinearsimple_copy_copy_test_highlight_crop.png}
\includegraphics[width=0.7\textwidth]{feed/feedbackplotlinearsimple.pdf}
\caption{Evolution of instantaneous ground state population for a total time of $t_{f} = 100 \mu s$ and temperature $2.6$~GHz, with feedback applied with some different delays $\tau$ ($0,5,10,30 \mu s$). Results are obtained from averaging 5$k$ quantum trajectories of the quantum feedback algorithm. Error bars are $2\sigma$ over $5k$ trajectories.}
\label{fig:feedbackplotlinearsimple}
\end{figure}

For the purpose of QA we want to maximize the final ground state population, i.e. ground state population at $t_f$. This serves as a metric to measure the effectiveness of the quantum feedback error correction. Note that the final ground state population need not decrease as delay increases. There is an optimal $τ$ we can search for. In Fig.~\ref{fig:feedbackdeltalinear_new}, we can see that there is a nonzero delay $τ$ where final ground state population reaches a local maximum. Thus there is an \textit{experimentally allowed} delay time $τ$ where the final ground state population is maximized.

\begin{figure}[h!]
\centering
\includegraphics[width=0.7\textwidth]{feed/feedbackdeltalinear_new.pdf}
\caption{Linear Schedule. Final ground state population. Single qubit. Total annealing time $t_f = 100\mu s$. Dashed purple line is the final ground state population without feedback. Average over $5k$ trajectories with Lindbladian feedback. There is a non-small (thus experimentally allowed) $\tau \approx 60us$ where the Final GS population is maximized.}
\label{fig:feedbackdeltalinear_new}
\end{figure}

\section{Feedback superoperator $\mathcal{L}_f$: Hamiltonian and basis of feedback}
$ℒ_f$ can also take the form of Hamiltonian~\cite{Wiseman-Milburn-book}.
\begin{equation}
\mathcal{L}_f[\,\cdot\,] = -i\left[H_{\text{fb}}(t), \,\cdot\,\right] \,,
\end{equation}
where $H_fb(t)$ is called the feedback Hamiltonian. Some examples of $H_fb(t)$ are
\begin{itemize}
\item $\sigma^x$ in the energy eigenbasis
\begin{equation}
H_{\text{fb}}(t) =  \frac{\pi}{2}\left(| \eps_{0}(t-\tau)   \rangle\langle \eps_{1}(t-\tau)  | + | \eps_{1}(t-\tau)   \rangle\langle \eps_{0}(t-\tau)  |\right) \,,
\label{eq:feedbackHenergy}
\end{equation}
\item $\sigma^x$ in the computational basis
\begin{equation}
H_{\text{fb}}(t) = H_{\text{fb}} = \frac{\pi}{2} \sigma^x \,,
\label{eq:feedbackHsigmax}
\end{equation}
\item $\sigma^z$ in the computational basis
\begin{equation}
H_{\text{fb}}(t) = H_{\text{fb}} =  \frac{\pi}{2} \sigma^z \,.
\label{eq:feedbackHsigmaz}
\end{equation}
\end{itemize}

These are Pauli operators in different bases. Note that in the decoherence model of the weak coupling limit, any excitation/relaxation takes place in the energy eigen-basis~\cite{ABLZ:12-SI}, thus one would expect that the feedback Hamiltonian $σ^x$ in the energy eigenbasis performs better than $σ^x$ in the computational eigenbasis (see simulations in Fig.~\ref{fig:Hfeedall}). Note also that the significance of the $π/2$ implies that the feedback Hamiltonian acts instantly (pulse/gate-like).

Another advantage of the feedback Hamiltonian formalism is that for delay $τ→0$, one can obtain the feedback master equation that explains the ensemble average of the effect of feedback. The feedback master equation is
\begin{align}
\frac{d}{dt} \rho(t) &= -i \left[H_{\textrm{S}}(t) + H_{\textrm{LS}}(t), \rho(t) \right] \notag\\
&+ \left(
e^{-iH_{\text{fb}}(t)}A_{-\omega}(t)\rho(t) A^\dagger_{-\omega}(t)e^{iH_{\text{fb}}(t)}  - \frac{1}{2} \left\{ A_{-\omega}^{\dagger}(t) A_{-\omega}(t) , \rho(t) \right\} \right) \notag\\
&+ D[A_{0}(t)]\rho(t) + D[A_{\omega}(t)]\rho(t)
\label{Eq:feedbackME_H}
\end{align}
This also allows direct master equation simulation, which can be faster than trajectories simulations for a small number of qubits. Fig.~\ref{fig:Hfeedall} below compares the two approaches of master equation simulation and the quantum feedback trajectories algorithm, over the three feedback bases mentioned in Eqs.~\eqref{eq:feedbackHenergy},~\eqref{eq:feedbackHsigmax} and~\eqref{eq:feedbackHsigmaz}.

% \begin{figure}[!h]
%     \centering
%     \hspace*{\fill}
%     % \begin{subfigure}[t]{0.4\linewidth}
%     % \centering
%     \includegraphics[width=\linewidth]{feed/energysigmaxinstant_crop.pdf}
%     %\includegraphics[width=\linewidth]{energysigmaxinstant_crop_2.pdf}
%     % \caption{}
%     % \label{Fig:2a}
%     % \end{subfigure}
%     % \hfill
%     % \begin{subfigure}[t]{0.4\linewidth}
%     % \centering
%     \includegraphics[width=\linewidth]{feed/sigmaxinstant_crop.pdf}
%     %\includegraphics[width=\linewidth]{sigmaxinstant_crop_2.pdf}
%     % \caption{}
%     % \label{Fig:2b}
%     % \end{subfigure}
%     % \hspace*{\fill}
%     \caption{Linear Schedule. Total annealing time $t_f = 10\mu s$. Instant feedback using (a) $\sigma^{x}$ operator in the instantaneous energy eigenbasis, (b) $\sigma^{x}$ operator in the computational basis. Average over 5000 trajectories. It can also be explained in feedback ME Eq.~\eqref{Eq:feedbackME_H} with feedback Hamiltonian (a) Eq.~\eqref{eq:feedbackHenergy} ($\tau = 0$), (b) Eq.~\eqref{eq:feedbackHsigmax}. QT: quantum trajectories approach. ME: feedback master equation approach.}
%     \label{Fig:2}
% \end{figure}

\begin{figure}[h!]
\centering
% \includegraphics[width=0.65\textwidth]{sigmaxinstant_crop_2.pdf}
\includegraphics[width=0.82\textwidth]{feed/sigmaxinstant.pdf}
\caption{Evolution of ground state population under feedback error correction of different feedback Hamiltonians ($\tau = 0$). The results of quantum jump trajectories (averaged over $5k$ trajectories) of quantum feedback algorithm agree with the feedback master equation solutions. Error bars are $2\sigma$ over $5k$ trajectories.}
\label{fig:Hfeedall}
\end{figure}

% \section{The need for a general QA feedback algorithm}
% Note that Eq.~\eqref{Eq:implicitfeedbackevolution} only describes the effect of feedback. One will need to include the original system evolution into account. One difference between QA and traditional gate-based quantum computing is that QA relies on a time-dependent annealing Hamiltonian. The system Hamiltonian being time-dependent implies there is an additional freedom to annealing time $t_f$, i.e. the annealing schedule ($A(t)$ and $B(t)$, for example in Fig.~\ref{fig: DWavenasaschedule}).
% %As mentioned, one measure of the effectiveness of the quantum error correction is the ground state population at the end of the annealing. Thus one may want
% In QA we wish to maximize the final ground state population, i.e.
% \begin{align*}
% &{\phantom{sssss}}\max_{\boldsymbol{\gamma}}\bra{\epsilon_{0}(t_{f})}\rho_{I}(t_{f})\ket{\epsilon_{0}(t_{f})} \\
% &\text{subject to: } \text{experimental constraints}
% \end{align*}
% \(\boldsymbol{\gamma}\) represents a list of hyperparameters, for example, delay $\tau$, the type of feedback, basis of feedback operator, annealing schedule, etc.
% %~(see Table.~\ref{Tab:hyper}).
% In the linear schedule example of Fig.~\ref{fig:feedbackdeltalinear_new}, only delay $\tau$ varies and other parameters are held constant. However, since there are many hyperparameters in general, the optimization problem would be difficult. Other techniques such as machine learning and deep learning would be useful.

% % \begin{table}[h!]
% % \centering
% % \begin{tabular}{ |c| }
% %  \hline
% %  Hyperparameters \\
% %  \hline
% %  Delay $\tau$  \\
% %  Annealing schedule  \\
% %  Type of feedback  \\
% %  Basis of feedback operators  \\
% %  \hline
% % \end{tabular}
% % \caption{A table of some examples of hyperparameters for QA quantum feedback optimization problem.}
% % \label{Tab:hyper}
% % \end{table}

\chapter{Reverse Annealing: Iterative reverse annealing}
\label{chap:ira}
We want to apply our open-system simulation techniques the study of open-system effects on a particular QA application called reverse annealing (RA). Reverse annealing has been predicted to have better performance than forward annealing in several systems~\cite{perdomo:sombrero,chancellor:reverse,King:2018aa,Ottaviani2018,nishimori:reverse-pspin, nishimori:reverse-pspin-2,Venturelli2019,marshall}.

Reverse annealing can be generally separated into two formalisms: iterated reverse annealing (IRA) and adiabatic reverse annealing (ARA). ARA is similar to standard quantum annealing, but with an extra initialization Hamiltonian. IRA does not have such term but starts at $s=1$ with a random Initial state that is usually an eigenstate of the problem Hamiltonian.

Iterated reverse annealing (IRA) involves (1) annealing backward from a known classical state to a mid-anneal state of quantum superposition, (2) searching for optimum solutions at this mid-anneal point while in the presence of an increased transverse field (quantum state), and then (3) proceeding forward to a new classical state at the end of the anneal.

In this chapter, we study the open-system effects on iterative reverse annealing and perform simulations of IRA in the weak coupling limits to bath, for $p$-spin model with $p=3$. This chapter is based on~\cite{Passarelli2019}. In Chapter.~\ref{chap: IRA_exp}, we perform the IRA simulations with more realistic experimental settings of D-Wave annealers, for $p$-spin model with $p=2$. In Chapter.~\ref{chap: ara}, we study the open-system effects on adiabatic reverse annealing (ARA), for $p$-spin model with $p=3$.

%  Passareli et al. \cite{Passarelli2019} tested the effectiveness of reverse annealing on the $p$-spin models with $p=3$. Numerical calculations using the time-dependent Lindblad equation~\cite{Albash2012, Yip2018} showed that the presence of coupling with the boson field significantly changes the performance of reverse annealing. This result shows that the interaction with the environment can lead to an improvement in performance. Since qubits and the environment should always interact with each other on the hardware, it is likely that this result reflects real physical processes in the D-Wave device. However, it is non-trivial if the coupling model  to the environment of Ref.~\cite{Passarelli2019} is good enough to explain experimental data.

\section{Introduction}
In iterative reverse annealing, the initial state is an eigenstate of the final problem Hamiltonian and the transverse field is cycled rather than strictly decreased as in standard (forward) quantum annealing. We present a numerical study of the reverse quantum annealing protocol applied to the $p$-spin model ($p=3$), including pausing, in an open system setting accounting for dephasing in the energy eigenbasis, which results in thermal relaxation. We consider both independent and collective dephasing and demonstrate that in both cases the open system dynamics substantially enhances the performance of reverse annealing. Namely, including dephasing overcomes the failure of purely closed system reverse annealing to converge to the ground state of the $p$-spin model. We demonstrate that pausing further improves the success probability. The collective dephasing model leads to somewhat better performance than independent dephasing. The protocol we consider corresponds closely to the one implemented in the current generation of commercial quantum annealers, and our results help to explain why recent experiments demonstrated enhanced success probabilities under reverse annealing and pausing.

\section{Iterative reverse annealing: the protocol}
\label{sec:reverse}
\begin{figure*}[t]
\subfigure[]{\includegraphics[width = 0.5\columnwidth]{ra1/schedules}\label{fig:schedules}}
\subfigure[]{\includegraphics[width = 0.5\columnwidth]{ra1/s-of-t}\label{fig:s-of-t}}
\caption{(a) Annealing schedules (in units such that $\hbar=1$) as a function of the annealing fraction $s(t)$, chosen to be similar to the schedules of the D-Wave processors.
(b) Annealing fraction $ s(t) $. The blue solid curve represents standard forward quantum annealing of total annealing time $\tau = 500$ns.
% 	The red dashed curve represents reverse annealing of the same total duration with an inversion point $ (\tinv = 200{\nano\second}, \sinv = 0.6) $.
\end{figure*}

Standard quantum annealing aims at solving optimization problems by employing quantum fluctuations that are slowly decreased to zero, to efficiently explore the solution space~\cite{kadowaki:qa, Santoro, albash:review-aqc}. A system of $ $ qubits is prepared at $ t = 0 $ in the ground state of a transverse field Hamiltonian, \ie, the uniform superposition over the $= 2^$ computational basis states $ {|0⟩, |1⟩, …, |- 1⟩} $.
%The parameter $ \Gamma $ represents the strength of the transverse field.
The magnitude of transverse field
%$ \hamtf $
is then slowly decreased, while the magnitude of the Hamiltonian $ $, encoding the optimization problem, is simultaneously increased. The adiabatic theorem guarantees that if the evolution is long on the timescale set by the inverse of the minimal gap $ = min_t [E_1(t) - E_0(t)] $ between the ground state and the first excited state (we set $ ℏ= 1 $ henceforth), then at the end of the anneal $ t = $ the system populates the target ground state of $ $ with a probability $ $ that approaches unity~\cite{Jansen:07,lidar:102106}. However, any finite sweep rate leads to diabatic Landau-Zener transitions at the avoided crossings~\cite{Joye:LZ}, thus reducing the success probability $ $ of the adiabatic algorithm. Therefore, the output state of quantum annealing is in general a trial solution of the optimization problem, ideally having a large overlap with the target ground state.

Reverse annealing instead aims at refining an already available trial solution~\cite{perdomo:sombrero}. For instance, NP-hard optimization tasks are solved using heuristics, whose output is often an approximation of the true global minimum. The algorithm of reverse annealing is the following.
\begin{enumerate}
\item At $ t = 0 $, the system is prepared in the trial solution state.
\item Quantum fluctuations are increased, until a turning point $ \tinv $ is reached during the dynamics. This ends the reversed part of the dynamics.
\item After the turning point, the dynamics follows the standard quantum annealing schedule; quantum fluctuations are decreased until $ t = \tf $, when the state is eventually measured.
\end{enumerate}
Careful choices of the turning point, and of the initial state, can lead to an improvement in the solution. Moreover, this scheme can also be repeated multiple times, each time starting from the output of the previous stage; hence the terminology of iterated reverse annealing. Assuming that each iteration improves the success probability $ $, this procedure can systematically improve the quality of the solution.

We consider the following time-dependent Hamiltonian, suitable for IRA but not ARA (which includes another term):
\begin{equation}
\label{eq:time-dependent-hamiltonian}
\ham(t) = A\bigl(s(t)\bigr) \hamtf + B\bigl(s(t)\bigr) \hamtarget, \quad  \hamtf = -\sum_{i=1}^{\nspin} \sigma_i^x .
\end{equation}
Since we only consider the IRA protocol, henceforth we simply refer to it as reverse annealing. The function of time $ s(t) $ in Eq.~\eqref{eq:time-dependent-hamiltonian} is the annealing fraction (or dimensionless time) and satisfies $ 0 ≤s(t) ≤1 $ for all $ t $. The two functions $ A(s) $ and $ B(s) $ determine the annealing schedule, which we choose to match the annealing schedule of the D-Wave processors [see Fig.~\ref{fig:schedules}]. They satisfy $ A(0) ≫B(0) $ and $ B(1) ≫A(1) $.

The functional form of $ s(t) $ distinguishes between standard (forward) and reverse annealing. In standard quantum annealing the dimensionless time is defined as $ s(t) = t / $, $ $ being the annealing time. Thus, $ s(t) $ is a monotonic function of $ t $, and is represented in the plane $ (t, s) $ by a straight line going from $ (0, s_0) $ to $ (τ, s_1) $, where $ s_0 = s(0) = 0 $ and $ s_1 = s(τ) = 1 $. This is shown in Fig.~\ref{fig:s-of-t} using a blue solid line, for $ τ= 500$ ns. During standard quantum annealing, quantum fluctuations are very large at $ t = 0 $, and decrease monotonically until $ t = τ$.

In contrast, in reverse annealing $ s_0 = s_1 = 1 $. Starting from $ s = s_0 $, where quantum fluctuations are zero, the annealing fraction is first decreased until it reaches the inversion point, $ s = $, at a time $ t = $. In this first branch, quantum fluctuations are increased. At $ s = $, the annealing fraction is then increased towards $ s = s_1 $, and quantum fluctuations are decreased again to zero. In Fig.~\ref{fig:s-of-t}, we show a typical function $ s(t) $ for a reverse annealing of annealing time $ = 500 $ with an inversion point $ (= 200, = 0.6) $, using a red dashed line.

In general, $ $ and $ $ can be chosen independently of each other. However, in this work we choose the following linear relation, in order to have only one free parameter:
\begin{equation}
\label{eq:inversion-time}
\tinv = \tf (1 - \sinv)\ , \ \ \sinv \ne 0, 1.
\end{equation}
In this way, we have that
\begin{equation}\label{eq:schedule}
s(t) =
\begin{cases}
1-t / \tf & \text{for $ t \le \tinv $},\\
\frac{1 - \sinv}{\tf \sinv} t + \frac{2 \sinv - 1}{\sinv} & \text{for $ t > \tinv $}.
\end{cases}
\end{equation}
%To avoid ill-defined cases, in the following we will only consider $ \sinv \ne 0, 1 $.

Another possible choice would be to fix $ $ so that the two slopes are the same, \ie, $ = / 2 $ for all choices of $ $. This is similar in spirit to what is discussed in Ref.~\cite{nishimori:reverse-pspin-2}.

%	Simulating quantum systems on classical hardware requires an exponentially large amount of computational resources, due to the dimension of the Hilbert space growing exponentially as a function of the number of qubits. Symmetries can often be exploited to restrict to interesting subspaces of the Hilbert space, having smaller dimensions. In this regards, having access to toy models with workable exact solutions and useful symmetries is paramount.
In what follows, we will focus on the fully-connected ferromagnetic $ p $-spin model, a model with a permutationally invariant Hamiltonian and a nontrivial phase diagram, often used as a benchmark for the performance of quantum annealing~\cite{gross:p-spin, derrida:p-spin}.
%We will introduce it in the next Section.

\section{Ferromagnetic $p$-spin model}\label{sec:pspin}

The Hamiltonian of the ferromagnetic $ p $-spin model, in dimensionless units, reads
\begin{equation}\label{eq:pspin}
\hamtarget = -\frac{\nspin}{2}{\left(\frac{1}{\nspin} \sum_{i = 1}^{\nspin} \sigma^z_i\right)}^p,
\end{equation}
with $ p ≥2 $. For even $ p $, there are two degenerate ferromagnetic ground states, whereas for odd $ p $ the ferromagnetic ground state is nondegenerate. For $ p = 2 $, the Hamiltonian of Eq.~\eqref{eq:time-dependent-hamiltonian} is subject to a second-order QPT in the thermodynamic limit at the critical point $ $ (or, equivalently, $ = $), separating a para- and a ferromagnetic phase. For $ p > 2 $, the QPT is first-order~\cite{bapst:p-spin}. The presence of QPTs affects also the finite-size behavior of the system, as the minimal spectral gap $ $, found at $ $ ($ = $), closes as $ ^-1/3 $ for $ p = 2 $ or exponentially in $ $ for $ p > 2 $. $ $ (\ie, $ $) approaches $ $ (\ie, $ $) as $ →∞$. First-order QPTs are especially detrimental for quantum annealing, as the annealing time has to grow exponentially with the system size to compensate the closure of the gap at $ s = $. In the following, we will focus on the case $ p = 3 $.

We can define the total spin operator $𝐒 = (S_x,S_y,S_z)$ and the dimensionless magnetization operators
\begin{equation}
m_{\alpha} = \frac{1}{\nspin} S_{\alpha} ,\quad S_{\alpha} = \sum_{i = 1}^{\nspin} \sigma_i^{\alpha}, \quad \alpha\in\{x,y,z\}
\end{equation}
that allow the rewriting of the time-dependent Hamiltonian of Eq.~\eqref{eq:time-dependent-hamiltonian} as
\begin{equation}
\ham(t) = -\frac{A(t)}{2} S_x - \frac{B(t)N}{2} m_z^p .
\end{equation}
%	 \ham(t) = -\frac{A(t)}{2} S_x - \frac{B(t)}{2 \nspin^{p - 1}} S_z^p .
%with $ \Gamma = 1 $.
Since $[𝐒^2,S_z]=0$ and both are Hermitian, they share an orthonormal basis ${|S,μ_S⟩}$ such that the eigenvalues of $𝐒^2$ are $S(S+1)$ with $S∈{0,1/2,1,…,/2}$ for even $ $ and $S∈{1/2,1,…,/2}$ for odd $ $, and the eigenvalues of $S_z$ are $μ_S ∈{-S,-S+1,…,S}$.
In the subspace with maximum spin $S=/2$, we instead label the basis states as $ |w⟩ ≡|/2 - w⟩ $, with $ w ∈{0, 1, …, } $. These are the eigenstates of $ m_z $ with eigenvalues $ m = 1 - 2 w / $. The target state is the ferromagnetic ground state $ |0⟩ $, \ie, the eigenstate of $ m_z $ with eigenvalue $ m = 1 $.

The Hamiltonian of the $ p $-spin model commutes with $𝐒^2$. Hence sectors differing by $S$ do not become coupled under the dynamics generated by the $ p $-spin model Hamiltonian. Since the ferromagnetic ground state and the initial one belong to the subspace with maximum spin $ S = / 2 $, the interesting dynamics occurs in this subspace, whose dimension scales \emph{linearly} with the number of qubits: $ = + 1 $. This fact enables us to perform numerical calculations with relative large numbers of qubits $$.

In the following, we will start reverse annealing in each of the $ $ excited states $ {|1⟩, …, |⟩ }$ of $ $ in the symmetric sector with $ S = /2 $. The similarity to the ferromagnetic ground state is quantified by the corresponding starting eigenvalue of $ m_z $, denoted $ m_0 $. Note that the $ w $th excited state differs from the ferromagnetic ground state by $ w $ spin flips. Therefore, the initial state and the target solution differ by a fraction $ c = _↑/ = 1 - w /$ of up-aligned qubits. These parameters are also related to the Hamming distance $ dH $,
%\ie, the number of $ -1 $ of the state in the single qubit basis. In fact,
via $ dH = - _↑= (1 - c) $.

\section{Unitary dynamics}
\label{sec:unitary}

% \begin{figure*}[tb]
\begin{figure*}[t]
\subfigure[]{\includegraphics[width = 0.48\columnwidth]{ra1/unitary-tf-100}\label{fig:unitary-tf-100}}
\subfigure[]{\includegraphics[width = 0.48\columnwidth]{ra1/unitary-scaling}\label{fig:unitary-scaling}}
\caption{(a) Success probability in unitary reverse annealing as a function of the inversion point $ \sinv $, for several values of the magnetization of the initial state. The dashed vertical line indicates $ \sinv = \sgap \approx 0.309 $. The annealing time is $ \tf = \SI{100}{\nano\second} $. We sampled the interval $ \sinv \in \rngopen{0}{1} $ using a step size of $ \Delta s = 0.002 $, and repeated the dynamics for each choice of $ \sinv $. (b) Maximum success probability achievable with unitary reverse annealing, as a function of the magnetization of the initial state, for two annealing times: $ \tf = \SI{100}{\nano\second} $ and $ \tf = \SI{1000}{\nano\second} $. The inset zooms in on the region $ m\in\rng{-1.0}{0.0} $. Note the logarithmic scale on the vertical axis.}
\end{figure*}

In this section we study the closed system case of a system of $ = 20 $ qubits, with $ p = 3 $. For our choice of parameters and in terms of the annealing schedules shown in Fig.~\ref{fig:schedules}, the $ p $-spin system has a minimal gap $ ≈2.45 $ at $ ≈0.309 $. The annealing time is $ = 100 $.
%We start without including system-environment interactions.

In Fig.~\ref{fig:unitary-tf-100}, we report the ground state population $ $  at $ t = $, as a function of the inversion point $ $, for several initial states: $ m_0 = 0.9;0.8;0;-1 $. Recall that the target ground state $ |0⟩ $ has $m=1$. We focus on the region $ ∈0.00.5$.

The rightmost part of Fig.~\ref{fig:unitary-tf-100} corresponds to cases in which the anneal is reversed too early, \ie, for $ < $ and $ > $. The system does not cross its quantum critical point, and the success probability is zero. Therefore, no effects on the outcome of the procedure are visible, as the dynamics is slow compared with the minimal inverse level spacing and diabatic transitions are exponentially suppressed. Thus, the system is forced to stay in its initial state, or transition to other excited states. In fact, avoided crossings between pairs of excited eigenstates occur at $ s > $ for this model, and Landau-Zener processes can further excite the $ p $-spin system.

On the other hand, if $ < $ the system crosses the minimal gap twice. Here, the success probability benefits from Landau-Zener processes, inducing transitions towards the ground state. In this region, we also note some non-adiabatic oscillations of the success probability, due to the finite annealing time. These oscillations are more evident for large $ m_0 $. As expected from the adiabatic theorem, they are suppressed for longer annealing times. For instance, we verified that they are no longer visible for $ = 1000 $ (not shown). The sharp rise of the success probability for $ m_0 = 0.9 $ occurs exactly at $ = $. For smaller values of $ m_0 $, the success probability rises more smoothly, as the ground state is reached after a preliminary sequence of Landau-Zener transitions between pairs of excited states, whose corresponding avoided crossings occur at $ s > $. For $ m_0 = 0.8 $, a very small rise of the success probability can still be observed around $ = $. This is due to the fact that during the reverse annealing, the system prepared in the second excited state first encounters an avoided crossing with the first excited state, where part of the population is transferred to the latter, and then the avoided crossing with the ground state, where the system populates its ground state. After reversing the dynamics, the two avoided crossings are encountered again (in the reverse order) and part of the population gets excited, thus reducing the success probability $ $.

As expected, reverse annealing is more effective when the initial state is close to the correct ground state. Moreover, as is also clear from Fig.~\ref{fig:unitary-tf-100}, the inversion time $ $ must be increasingly close to $0$ for decreasing $ m_0 $, in order to obtain a nonzero success probability at $ t = $. This means that almost the entire dynamics is spent in the reverse part of the annealing, and the system is eventually quenched towards $ s = s_1 $ for $ t ≈$. Even so, if the initial state is too far in energy from the correct solution, the success probability of reverse annealing is always close to zero, as evident from the curves for $ m_0 = 0 $ and $ m_0 = -1 $ in Fig.~\ref{fig:unitary-tf-100}.

The maximum success probability decreases rapidly as a function of $ m_0 $. This is clearly seen in Fig.~\ref{fig:unitary-scaling}, where we report the maximum attainable success probability as a function of $ m_0 $, for annealing times $ = 100;1000 $.
%The maximum success probability is also affected by the total annealing time. In fact,
Increasing the annealing time reduces non-adiabaticity and results in a lower success probability, compared with that at the end of a faster reverse anneal. As shown in Fig.~\ref{fig:unitary-scaling}, which zooms in on the region $ m_0 ∈-10 $, this decrease can be of several orders of magnitude for poorly chosen trial solutions.  The influence of the annealing time is less pronounced close to $ m_0 = 1 $, and more evident for intermediate and lower values of $ m_0 $. This is consistent with the adiabatic theorem, since a longer anneal time guarantees that the system will have a higher probability of remaining close to the initial eigenstate it has the largest overlap with (not necessarily the ground state)~\cite{Jansen:07}.

Adopting a conventional quantum annealing procedure, the success probability in the case of $ = 100 $ would be $ P_0 = 0.96 $. This value is larger than any $ P_0 $ achievable using reverse annealing. However, this argument cannot be used to discredit reverse annealing for two main reasons. First, it is clear that, in the analyzed case, we are very close to adiabaticity, where conventional annealing is efficient. Second, as clarified in the next Sections, the role of dissipation may strongly affect this scenario.

The results of this Section are in agreement with those reported in Ref.~\cite{nishimori:reverse-pspin-2}. Namely, as is clear from Fig.~\ref{fig:unitary-scaling}, upon iteration the IRA protocol will only decrease the success probability under unitary, closed system dynamics, unless the initial state was already chosen as the solution of the optimization problem.

\begin{figure*}[t]
\subfigure[]{\includegraphics[width = 0.48 \columnwidth]{ra1/dissipative-tf-100}\label{fig:dissipative-tf-100}}
\subfigure[]{\includegraphics[width = 0.48 \columnwidth]{ra1/dissipative-tf-500}\label{fig:dissipative-tf-500}}
\caption{Success probability in reverse annealing as a function of the inversion point $ \sinv $, for several values of the magnetization of the initial state. The $ p $-spin system is coupled to a collective dephasing bosonic environment as in Eq.~\eqref{eq:system-bath-collective}, and the coupling strength is $ \eta = \num{1e-3} $. The dashed vertical line denotes the time $ \sgap $ of the avoided crossing between the ground state and the first excited state. In (a) the annealing time is $ \tf = \SI{100}{\nano\second} $, in (b) $ \tf = \SI{500}{\nano\second} $. We sampled the interval $ \sinv \in \rngopen{0}{1} $ using a step $ \Delta s = 0.02 $. All other parameters are given in the main text.}
\label{fig:dissipative}
\end{figure*}

\section{Open system dynamics subject to dephasing-induced relaxation}\label{sec:dissipative}

Physical quantum processors always interact with the surrounding environment, which induces decoherence and thermal excitation/relaxation, which in turn impacts the performance of quantum annealing~\cite{childs:robustness,amin_decoherence_2009,Albash:2015nx}. We
%use a description of the bath in terms of harmonic oscillators, and we
assume weak coupling between the qubit system and the environment.
%We formulate the dynamics of the reduced system in terms of a quantum dynamical semigroup.
It can then be shown that the reduced system density matrix evolves according to a quantum master equation in time-dependent Lindblad form~\cite{zanardi:master-equations}, known as the adiabatic master equation:
\begin{equation}\label{eq:lindblad}
\frac{d\rho(t)}{dt} = \iu \bigl[\rho(t), \ham(t) + \ham\ped{LS}(t)\bigr] + \diss\bigl[\rho(t)\bigr].
\end{equation}
In Eq.~\eqref{eq:lindblad}, $ LS(t) $ is a Lamb shift term and $ $ is the dissipator superoperator, which makes the dynamics non-unitary and irreversible. They are expressed in terms of Lindblad operators, inducing dephasing or quantum jumps (pumps and decays) between pairs of adiabatic energy eigenstates~\cite{yip:mcwf}:
\begin{align}
\ham\ped{LS} &= \sum_{ab} L_{ab}^\dagger(t) L_{ab}(t) S_{ab}(\omega_{ab});\\
\diss[\rho(t)] &= \sum_{ab} \gamma_{ab} \biggl( L_{ab}(t) \rho(t) L_{ab}^\dagger(t) \notag\\
&\quad- \frac{1}{2} \left\lbrace L_{ab}^\dagger(t) L_{ab}(t), \rho(t) \right\rbrace \biggr).
\end{align}
The Lindblad operators are determined by the instantaneous eigenbasis of the system Hamiltonian $H(t)$ and the system-bath interaction Hamiltonian $ $. In general, this coupling Hamiltonian involves local operators
%, \ie, ``position'' degrees of freedom both for the qubits in real space (\eg, $ \sigma^z_i $), and for the harmonic oscillators in  momentum space (\eg, $ a_{k, i} + a_{k, i}^\dagger $). This generic form, however,
that break the spin symmetry of the $ p $-spin model, as each qubit is then coupled to its own bath. We study both this independent decoherence model and the collective decoherence model, wherein  %Therefore, we will adopt a specific form for the interaction potential, where
all the qubits are coupled to a collective bath with the same coupling energy $ g $, in order to preserve the spin symmetry. More specifically, we first consider collective dephasing, for which the system-bath coupling Hamiltonian is
\begin{equation}\label{eq:system-bath-collective}
\hamsysbath\api{col} = g S_z
%\sum_k (a_k + a_k^\dagger),
\otimes B ,
\end{equation}
where $B$ is a bath operator [e.g., $B = ∑_k (a_k + a_k^†)$ for an oscillator bath with annihilation operators $a_k$ for the $k$th bosonic mode].
The Lindblad operators are represented in the instantaneous energy eigenbasis of $H(t)$ as~\cite{zanardi:master-equations}:
\begin{equation}\label{eq:Lab(t)}
L_{ab}(t) = \braket{E_a(t)| S_z |E_b(t)} \ket{E_a(t)}\bra{E_b(t)}.
\end{equation}
This represents collective dephasing in the energy eigenbasis, wherein the dephasing process randomizes the relative phase between eigenstates of the system Hamiltonian. Thus, this model does not support phase coherence between energy eigenstates.\footnote{It is worth pointing out a caveat. Namely, the collective dephasing model in general supports decoherence free subspaces (DFSs), \ie, subspaces that evolve unitarily despite the coupling to the bath~\cite{Zanardi:97c,Lidar:1998fk}. For instance, the $ S = 0 $ subspace (for even $ \nspin $) is a DFS of the $ p $-spin model. However, the $p$-spin model is unsuitable for performing quantum annealing inside a DFS, since its Hamiltonian consists of operators that preserve the DFS, so that no dynamics would take place if we were to try to encode a computation using states inside the DFS. Instead, to obtain meaningful dynamics (performing a computation) subject to the collective dephasing model, we would need to add Heisenberg exchange terms to the system Hamiltonian~\cite{Kempe:2001uq}.} As a consequence, thermal relaxation tends to equilibrate the system towards its Gibbs state, with a characteristic timescale set by the inverse of the bath spectral density at the gap frequency~\cite{Albash:2015nx}.

If instead the qubit system is coupled to independent, identical baths, the system-bath coupling operator becomes
\begin{equation}\label{eq:system-bath-independent}
\hamsysbath\api{ind} = g \sum_{i} \sigma_i^z \otimes B_i,
\end{equation}
where, e.g., in the bosonic case $B_i = ∑_k (a_k, i + a_k, i^†)$. Thermal relaxation effects occur here similarly to the collective dephasing case.
However, simulations in this case are more demanding due to the fact that the spin symmetry is broken and that we have $$ times as many Lindblad operators, \ie,
\begin{equation}
L_{ab, i}(t) = \braket{E_a(t)| \sigma_i^z |E_b(t)} \ket{E_a(t)}\bra{E_b(t)}.
\end{equation}
Therefore, in this case we will only investigate reverse annealing starting from the first excited state in the symmetric subspace with maximum spin, \ie, $ |w = 1⟩ $, for $ ∈{3,…,8}$. Moreover, for the particular cases of $ = 7 $ and $ = 8 $, we truncate our system to the lowest $ 29 $ and $ 37 $ eigenstates, respectively, to speed up the numerics. This choice is made since the first three levels of the maximum spin subspace at $ s = 1 $ are spanned by $ ∑_i=0^2 7i = 29 $ (for $ = 7 $) and $ ∑_i=0^2 8i = 37 $ (for $ = 8 $) energy eigenstates. We confirm that this is a good approximation by checking that the total population among these levels is close to $1$ during the reverse annealing when additional levels are included in the simulation.

\begin{figure}[t]
\centering
\includegraphics[width = 0.9 \linewidth]{ra1/individual-vs-collective-nopausing_copy_print.pdf}
\caption{Success probability in reverse annealing as a function of the inversion point $ \sinv $, for $\nspin \in \{3,\dots,8\}$. The initial state is the first excited state of the maximum spin subspace ($ m_0 = 1 - 2/\nspin $). The dashed vertical line denotes the time $ \sgap $ of the avoided crossing between the ground state and the first excited state. The annealing time is $ \tf = \SI{100}{\nano\second} $. We sampled the interval $ \sinv \in (0, 1) $ using a step size of $ \Delta s = 0.005 $.}
\label{fig:withoutpausing}
\end{figure}

\begin{figure}[h!]
\centering
\includegraphics[width=0.5\linewidth]{ra1/peak-fidelity_copy_print.png}
\caption{Maximum success probability achievable with reverse annealing as a function of the number of qubits $ \nspin $, using the collective and the independent dephasing models of Eqs.~\eqref{eq:system-bath-collective} and~\eqref{eq:system-bath-independent}, respectively. The annealing time is $ \tf = \SI{100}{\nano\second} $.}
\label{fig:fiddependenceonn}
\end{figure}

\begin{figure*}[h!]
\subfigure[]{\includegraphics[width = 0.49 \columnwidth]{ra1/dissipative-tf-100-pause-100}	\label{fig:dissipative-tf-100-pause-100}}
\subfigure[]{\includegraphics[width = 0.49 \columnwidth]{ra1/dissipative-tf-100-pause-400}\label{fig:dissipative-tf-100-pause-400}}
\caption{Success probability in paused reverse annealing for the collective dephasing model, as a function of the inversion point $ \sinv $, for several values of the magnetization of the initial state. The coupling strength is $ \eta = \num{1e-3} $. The annealing time is $ \tf = \SI{100}{\nano\second} $. In (a) a pause of duration $ \lpause = \SI{100}{\nano\second} $ is inserted at the inversion point, while in (b) $ \lpause = \SI{400}{\nano\second} $. All other parameters are given in the main text. Dotted lines with empty symbols refer to open system reverse annealing of time $ \tf = \SI{100}{\nano\second} $ and no pauses. The dashed vertical line denotes $ \sinv = \sgap $. The interval $ \sinv \in \rngopen{0}{1} $ is sampled using a step size $ \Delta s = 0.02 $.}
\label{fig:dissipative-tf-100-pause}
\end{figure*}

\begin{figure}[tb]
\includegraphics[width = 0.9 \linewidth]{ra1/individual-vs-collective-pausing_copy_print.pdf}
\caption{Comparison of the success probability in reverse annealing for the collective and independent dephasing models, as a function of the inversion point $ \sinv $, for $\nspin \in\{3,\dots,8\}$. The initial state is the first excited state of the maximum spin subspace ($m_0 = 1 - 2/\nspin$). The dashed vertical line denotes the time $\sgap$ of the avoided crossing between the ground state and the first excited state. The annealing time is $\tf = \SI{100}{\nano\second}$, and a pause of duration $ \lpause = \SI{100}{\nano\second}$ is inserted at the inversion point.}
\label{fig:withpausing}
\end{figure}

The adiabatic master equation in Eq.~\eqref{eq:lindblad} is unraveled using a time-dependent Monte Carlo wavefunction (MCWF) approach~\cite{yip:mcwf}. The advantage of MCWF is that it allows to work with wavefunctions rather than density matrices, thus saving quadratically in the dimension of the objects we need to store for numerical calculations. The tradeoff is that to recover the statistical properties of the density operator we need to average over a large number $ K $ of independent trajectories. For the collective system-bath coupling of Eq.~\eqref{eq:system-bath-collective}, the time evolution operator of each trajectory is generated by the effective non-Hermitian Hamiltonian
\begin{align}\label{eq:non-hermitian-hamiltonian}
\ham\ped{eff}(t) &= \ham(t) + \ham\ped{LS}(t) - \frac{\iu}{2} \sum_{a \ne b} \gamma_{ab} L_{ab}^\dagger(t) L_{ab}(t) \notag \\ & \quad-\gamma_0 \frac{\iu}{2} \sum_{ab}  L_{aa}^\dagger(t) L_{bb}(t),
\end{align}
where $ γ_ab $ and $ γ_0 $ are the rates for jumps and dephasing, respectively. They are related to the temperature and to the spectral density of the bosonic bath,
\begin{equation}\label{eq:spectral-density}
J(\omega) = g^2 \sum_k \delta(\omega -  \omega_k) = 2\pi \eta \omega e^{-\omega / \omegac},
\end{equation}
where $ ω_k $ are the bath eigenfrequencies, $ $ is a high-frequency cutoff and $ η$ is the dimensionless coupling strength. We fix $ = 1 $ and $ η= 1e-3 $. The working temperature is chosen to be $ T = 12.1 = 1.57 $, as in experimental quantum annealing systems~\cite{dwave-site}. Eq.~\eqref{eq:non-hermitian-hamiltonian} is easily extended to the independent dephasing case by including a summation over $ i $ for the Lindblad operators and in the Lamb shift term.

The time evolution operator generated by the non-Hermitian Hamiltonian of Eq.~\eqref{eq:non-hermitian-hamiltonian} is not unitary. Therefore, the norm of the wavefunction decays in time. Whenever the squared norm decreases below a randomly extracted threshold $ r ∈01 $, a quantum jump occurs, projecting the wavefunction on one of the eigenstates of the Hamiltonian $ (t) $.
The adiabatic master equation dynamics is found after averaging over all stochastic trajectories. In this work, we fix $ K = 5000 $ trajectories and consequently find relative Monte Carlo errors $ δO / O $ of the order of $ 1.5 $ over all observables $ O $.

We repeat the simulations we reported in Sec.~\ref{sec:unitary} for $ = 20 $, $ p = 3 $ and $ = 100 $, but now include the role of the environment. We consider both the collective and independent dephasing models.

In Fig.~\ref{fig:dissipative-tf-100}, we show the success probability as a function of the inversion time $ $, for the four initial magnetizations $ m_0 = 0.9;0.8;0;-1 $. Monte Carlo errors are of the order of the point size in all cases and are invisible.

As in the unitary case, if the inversion occurs too early (\ie, for $ ≪$, or, equivalently, $ ≫$), the reverse annealing protocol fails to find the ferromagnetic ground state. In fact, thermal excitations are suppressed, as well as Landau-Zener transitions, due to the large level spacing, compared with the temperature and the inverse of the annealing time. For $ ≈$ ($ ≈$), however, the scenario is drastically different from the unitary case of Fig.~\ref{fig:unitary-tf-100}.

The first difference is that the success probability can be nonzero even if the inversion occurs for $ ≳$, especially for $ m_0 = 0.9 $, where the tail of the curve extends to $ ≈0.75 $. When the instantaneous gap is of the same order of magnitude as the temperature, thermal processes influence reverse annealing even before crossing the minimal gap. Second, for all $ m_0 $ we observe a sudden increase in the success probability around $ ≈$, that eventually brings all curves to an almost flat region at $ < $, where the success probability reaches the large value $ ≈0.957 $. The value of the maximum success probability at the plateau is $ m_0 $-independent within Monte Carlo errors. The time at which the success probability starts to increase with respect to the baseline depends on $ m_0 $. Moreover, the flat region is wider for larger $ m_0 $, although it has a finite width for all $ m_0 $.

These results show that even trial solutions far in Hamming distance from the ferromagnetic ground state can result in a large success probability at the end of a reverse anneal. Moreover, the time window in which inverting the annealing favors the ferromagnetic ordering is relatively large.

We also studied a longer annealing time, $ = 500 $, as shown in Fig.~\ref{fig:dissipative-tf-500}.
%, where we show the success probability $ \pgs $ as a function of the inversion time, for $ m_0 = \text{\numlist{0.9;0.8;0;-1}} $.
Here, we note that the onset of the success probability plateau shifts towards longer values of $ $, compared with the $ = 100 $ case. Therefore, the plateau is wider, and the maximum success probability at the plateau is $ ≈1 $ within Monte Carlo errors for all $ m_0 $ we considered. This is in contrast with the unitary case of Fig.~\ref{fig:unitary-scaling}, where increasing the annealing time had detrimental effects on the algorithm. This evidence supports the notion that the success probability enhancement is due to thermal effects, rather than due to purely unitary quantum dynamics~\cite{passarelli:pspin}. Moreover, the adiabatic theorem for open quantum system guarantees convergence to the steady state of the superoperator generator of the dynamics in the large $$ limit~\cite{Venuti:2015kq,Venuti:2018aa}. This too helps to explain our observations: the steady state of the Davies-Lindblad generator of the open system dynamics we considered here is the Gibbs distribution of the final Hamiltonian, which at sufficiently low temperature relative to the gap is the ferromagnetic ground state. Recall that in our case $ ≈2.45 $ (at $ ≈0.309 $) and $ T = 1.57 $.

Comparing these results with conventional forward annealing in the presence of dissipation, we note that the success probability at the plateau is similar to that of standard quantum annealing for $ η= 1e-3 $ ($ P_0 = 0.98 $). The reason is that collective dephasing favors the ferromagnetic ordering in the $p$-spin model and its induced relaxation increases the success probability compared to the isolated case, in agreement with previous findings~\cite{passarelli:pspin}. We also observe that when the annealing is far from adiabaticity (e.g. $= 1$), reverse annealing becomes a more favorable approach to forward annealing for the $p$-spin model.

We also compare the collective and independent dephasing models of Eqs.~\eqref{eq:system-bath-collective} and~\eqref{eq:system-bath-independent}.
Fig.~\ref{fig:withoutpausing} shows the simulation results for the two models using the adiabatic master equation of Eq.~\eqref{eq:lindblad} for $ ∈{3,…,8} $. As shown in the figure, simulations using the collective dephasing model have larger success probabilities for almost every $ $. This is because, in the independent dephasing model, other states not in the subspace of maximum spin become accessible by thermal excitation or diabatic transition during the reverse anneal. For all of the system sizes we simulated, we had to reverse anneal to a smaller inversion point $ $ for the independent dephasing model to achieve the same success probability as the collective dephasing model. Moreover, the maximum success probability achievable is always smaller for the independent dephasing model. The success probabilities from both models, however, are very similar as $ →0 $, \ie, in the quench limit of the direct part of the evolution.

Figure~\ref{fig:fiddependenceonn} shows how the maximum success probability (over $ $) of both bath models depends on the number of qubits. As $ $ increases, the maximum success probability of the independent dephasing model decreases more rapidly than that of the collective dephasing model. While we can infer that if we modeled independent dephasing for $ = 20 $ we would not observe as large success probabilities as in Fig.~\ref{fig:dissipative}, we stress that  reverse annealing in the independent dephasing model still yields a significantly larger success probability  (for the same $$ values) than the unitary dynamics case described in Section~\ref{sec:unitary}.

\section{Open system dynamics with a pause}\label{sec:pausing}

Quantum annealing in the presence of a low temperature bath can benefit from pauses inserted at certain times during the dynamics~\cite{marshall, passarelli:pausing}. During a pause, $ s(t) = constant $, and the system evolves with a time independent Hamiltonian, subject to dephasing. When a pause is inserted some time after $ $, the environment favors a redistribution of the repopulation according to the Gibbs state at the pause point; at sufficiently low temperature (relative to the gap at this point), this can result in a repopulation of the instantaneous ground state. In this section, we show that pauses at the inversion point can further improve the performance of reverse annealing of the $ p $-spin model.

We repeat the simulations for $ = 20 $, $ p = 3 $ and $ = 100 $, using the collective dephasing model. A pause of duration $ = $ is inserted at $ t = $, so that the total annealing time, including the pause, is $ ' = + = 200 $.

In Fig.~\ref{fig:dissipative-tf-100-pause-100}, we report the success probability as a function of the inversion point, for starting magnetizations $ m_0 = 0.9;0.8;0;-1 $. We compare the paused case with the unpaused case, for which $ = 100 $.
As can be seen in the figure, if the dynamics is reversed too early ($ ≫$), the success probability at the end of the anneal vanishes. The level spacing is large compared with the temperature. The relaxation rate is small and the pause is too short to have impact on the dynamics.

However, the presence of a pause significantly changes the outcome of the annealing around $ ≈$.  In fact, when a pause is inserted at $ ≳$, the success probability reaches $ ≈1 $ for a wide range of inversion points and for all $ m_0 $, within Monte Carlo errors. Here, the ground state is completely repopulated by thermal relaxation. This is in contrast with conventional quantum annealing, where the success probability exhibits a peak as a function of the pausing time, when the pause is inserted about $ 20 $ later than $ $, and then rapidly returns to its baseline value~\cite{marshall, passarelli:pausing}.
In contrast, for $ < $, the effect of the pause is negligible; the solid (with pause) and dotted (no pause) lines in Fig.~\ref{fig:dissipative-tf-100-pause-100} overlap in this region.

We repeated our analysis for a pause duration $ = 400 $, with total annealing time $ ' = 500 $. As shown in Fig.~\ref{fig:dissipative-tf-100-pause-400}, the longer pause duration affects the results only marginally. Comparing with Fig.~\ref{fig:dissipative-tf-100-pause-100}, we note that the qualitative behavior of the curves is the same in the two cases. The pause duration affects mostly the region $ ≳$. A longer pause enhances thermal relaxation, thus the success probability starts to increase from its baseline earlier than for shorter $ $. This results in a wider plateau where the success probability is large, compared with Fig.~\ref{fig:dissipative-tf-100-pause-100}.

Finally, we compare the collective and independent dephasing models while including pausing,
%of Eqs.~\eqref{eq:system-bath-collective} and~\eqref{eq:system-bath-independent},
starting from the first excited state of the maximal spin sector.
%, for $ \nspin = \text{\numlist{3;4;5;6;7;8}} $.
The results are shown in Fig.~\ref{fig:withpausing}, for a pause of duration $ = 100 $ inserted at the inversion point. The collective dephasing model continues to exhibit higher success probabilities than the independent dephasing model, as in the case discussed in the previous section, but the results of the two models coincide when $ < $.
%Interestingly, with pausing the maximum success probabilities of the independent dephasing model can be restored to be as large as those of the collective dephasing model.
Thus, relaxation to the ground state during the pause improves performance for both dephasing models.
%The dependence of the maximum success probabilities on the number of qubits of both models with pausing is shown in Table~\ref{tab:numericaltable}. From the Table we also note
Note that, as $ $ increases, the maximum success probability of the collective dephasing model is achieved at $ > $, while it is achieved at $ < $ in the independent dephasing model. This is in agreement with the $ = 20 $ result shown in Fig.~\ref{fig:dissipative-tf-100}.

%	\footnotesize
%	\centering
%	\caption{Maximum success probability and gap size as a function of the number of qubits ($n$). The parentheses indicate at which $s_{\text{inv}}$ the maximum success probability is achieved.}
%	\begin{tabular}{c c c c c c c c}
%		\toprule
%		\multicolumn{2}{c}{$n$}  & $3$      &$4$ &$5$ &$6$ &$7$&$8$\\ \midrule
%		\multirow{2}[4]{*}{\specialcellc{No pausing: \\ Maximum success probability }}  & Individual bath&\makecell{$0.9335$ \\ $(s_{\text{inv}}=0.040)$} & \makecell{$0.9156$ \\ $(s_{\text{inv}}=0.055)$}  & \makecell{$0.8978$ \\ $(s_{\text{inv}}=0.065)$} & \makecell{$0.8749$ \\ $(s_{\text{inv}}=0.075)$} & \makecell{$0.8487$ \\ $(s_{\text{inv}}=0.085)$} & \makecell{$0.8198$ \\ $(s_{\text{inv}}=0.105)$}\\ \cmidrule{2-8}
%			& Collective bath &  \makecell{$0.9984$ \\ $(s_{\text{inv}}=0.065)$}& \makecell{$0.9976$ \\ $(s_{\text{inv}}=0.095)$} & \makecell{$0.9959$ \\ $(s=0.120)$}& \makecell{$0.9932$ \\ $(s_{\text{inv}}=0.145)$} & \makecell{$0.9882$ \\ $(s_{\text{inv}}=0.160)$} & \makecell{$0.9840$ \\ $(s_{\text{inv}}=0.195)$}\\ \midrule
%		\multirow{2}[4]{*}{\specialcellc{With pausing: \\ Maximum success probability }}  & Individual bath&\makecell{$0.9994$ \\ $(s_{\text{inv}}=0.085)$}  & \makecell{$0.9990$ \\ $(s_{\text{inv}}=0.100)$} & \makecell{$0.9983$ \\ $(s_{\text{inv}}=0.120)$} & \makecell{$0.9973$ \\ $(s_{\text{inv}}=0.150)$} & \makecell{$0.9941$ \\ $(s_{\text{inv}}=0.180)$} & \makecell{$0.9917$ \\ $(s_{\text{inv}}=0.190)$}\\\cmidrule{2-8}
%			&Collective bath&\makecell{$0.9991$ \\ $(s_{\text{inv}}=0.085)$}  & \makecell{$0.9979$ \\ $(s_{\text{inv}}=0.100)$} & \makecell{$0.9960$ \\ $(s_{\text{inv}}=0.120)$} & \makecell{$0.9963$ \\ $(s_{\text{inv}}=0.370)$} & \makecell{$0.9966$ \\ $(s_{\text{inv}}=0.370)$} & \makecell{$0.9972$ \\ $(s_{\text{inv}}=0.370)$}\\ \midrule
%		\multirow{2}[4]{*}{}  & $\Delta$ (GHz) & $7.7551$  & $6.9451$ & $6.3062$ & $5.7842$ & $5.3451$ & $4.9672$\\ \cmidrule{2-8}
%			& $s_{\Delta}$  & $0.3302$   & $0.3210$    & $0.3163$ & $0.3137$ &$0.3120$ &$0.3109$\\  \bottomrule
%	\end{tabular}
%	\label{tab:numericaltable}

\section{Conclusion of this chapter}
\label{sec:conclusions}

Earlier work revealed an intriguing tension between experimental results demonstrating a substantial enhancement in success probabilities for random spin glass instances under reverse annealing compared to standard (forward) annealing~\cite{marshall}, and theoretical results finding that reverse annealing adversely affects performance for the $p$-spin model, in a closed system setting~\cite{nishimori:reverse-pspin-2}. In this work we resolved this tension by performing a numerical study of reverse annealing of the $p$-spin model in an open system setting, where we included dephasing in the instantaneous energy eigenbasis. We found that the associated thermal relaxation results in significant increase in the success probabilities, as long as the inversion point of the reverse annealing protocol is chosen to be close to the avoided crossing point, or before it. Pausing at the inversion point further improves performance.

Since closed-system, unitary dynamics predicts that reverse annealing fails, yet its open system analogue succeeds, it follows that thermal relaxation is the mechanism responsible for the success. Reverse annealing is thus an example of a family of protocols that strictly benefit from thermal effects~\cite{verstraete2009quantum,Venuti:2017aa}.
It may be worth noting that quantum effects are likely to play an important role in thermal relaxation because the success probability is very small for $s_inv$ close to $1$, i.e. when the Hamiltonian stays classical during the anneal, even with a pause.
Whether this can lead to any quantum speedups is an interesting problem worthy of future investigations.

\section{Acknowledgments of this chapter}
This chapter was originally published in~\cite{Passarelli2019}. The simulation data of collective coupling model is provided by Gianluca Passarelli.
The research of KY, DL, and HN is based upon work (partially) supported by the Office of
the Director of National Intelligence (ODNI), Intelligence Advanced
Research Projects Activity (IARPA), via the U.S. Army Research Office
contract W911NF-17-C-0050. The views and conclusions contained herein are
those of the authors and should not be interpreted as necessarily
representing the official policies or endorsements, either expressed or
implied, of the ODNI, IARPA, or the U.S. Government. The U.S. Government
is authorized to reproduce and distribute reprints for Governmental
purposes notwithstanding any copyright annotation thereon. Computation for some of the work described in this chapter was supported by the University of Southern California Center for High-Performance Computing and Communications (hpcc.usc.edu).
\chapter{Experimental iterative reverse annealing}
\label{chap: IRA_exp}
% Iterated reverse annealing (IRA) is a formalism of reverse annealing implemented in current quantum annealers. The system is prepared in a trial solution state, annealed reversely to an inversion point, and then annealed forwardly. It may also be iterated with the last output state as new input.
% We perform IRA experiments on the D-Wave 2000Q device, with a focus on the p-spin problem where p=2. The encoding of problem Hamiltonian is made possible by minor embedding. We examine the dependence of IRA performance on parameters such as problem size, anneal time, inversion points; and also on D-Wave’s functionalities like state reinitialization, number of iterations, pausing. The performance is evaluated in total and partial success probabilities. We can find correct solutions with near certainty. Comparison with standard forward annealing is also made.
% We also perform open-system simulations with dephasing in the instantaneous energy eigen-basis. The experimental total success probabilities are well matched. A difference between independent and collective SB coupling model is reported due to spin-symmetry breaking of the classical input. Simulations also show IRA fails in a closed system, suggesting the importance of thermal relaxation in the convergence to the correct solution.
In this chapter, we continue to study the open-system effects on iterative reverse annealing. We perform the IRA simulations with more realistic experimental settings of D-Wave annealers, for $p$-spin model with $p=2$. We explore both quantum and semiclassical simulation methods, and also compare our simulation results to the experiments.
\section{Introduction}
Iterated reverse annealing (IRA) is a variant of quantum annealing, in which the system is prepared in a trial solution state, reverse-annealed to an inversion point, and then forward-annealed. It may also be iterated with the last output state as the input of the next iteration.
We perform experiments of reverse annealing on the D-Wave 2000Q device, with a focus on the $p$-spin problem with $p=2$. We examine the dependence of the performance on parameters such as problem size, annealing time and inversion points, as well as pausing and the number of iterations. The performance is evaluated in terms of total and partial success probabilities, the latter being the probabilities for each of two degenerate ground states. To explain the experimental results, we perform open-system simulations with decoherence in the instantaneous energy eigenbasis. The experimental total success probabilities  match the simulation data well. We find a difference between the independent and collective dephasing models, which is explained by non-conservation of spin symmetry of classical inputs. Simulations show that iterated reverse annealing fails in a closed system, suggesting the importance of thermal relaxation in the convergence to the correct solution. We also perform semi-classical simulations of spin-vector Monte Carlo to supplement  our understanding of the experimental data.

% Quantum annealing is a quantum metaheuristic to obtain the global optimal solution by using quantum fluctuations to escape from local minima in optimization problems~\cite{kadowaki:qa, Santoro, RevModPhys.80.1061, morita:mathematical-foundation-qa,albash:review-aqc, Hauke2020}.

Recently, D-Wave devices that actually realize quantum annealing on a hardware have been developed. Not only can the D-Wave device be used to solve optimization problems by following quantum annealing  physically, but it is also able to perform physics experiments, for example, spin-glass phase transitions~\cite{Harris2018}, the Kosterlitz-Thouless phase transition~\cite{King:2018aa, King:2019aa}, the Kibble-Zurek mechanism~\cite{Gardas2018}, the Griffiths-McCoy singularity \cite{Nishimura2020}, spin ice \cite{King2020}, and the Shastry-Southernland model \cite{Kairys2020}.% Thus, the quantum annealing device has the potential to be used as a physics experimental machine as well as for optimization problems.\par

As in the case of simulated annealing~\cite{kirkpatrick:sa}, which is a classical counterpart of quantum annealing, one of the simple ways to enhance the performance of quantum annealing is to control the ``annealing schedule''. For conventional quantum annealing~\cite{kadowaki:qa, Santoro}, we choose the ground state of the transverse field term as the initial state and define forward annealing as the temporal change of the system from the transverse field Hamiltonian to the target Hamiltonian. Traditionally, this transition process is supposed to be adiabatic. Diabatic transition processes have also been proposed to overcome the difficulties of adiabatic computing, see e.g., ~\cite{Brady2020,Mbeng2019,Crosson2020}. On the other hand, reverse annealing is defined as the process of choosing an appropriate classical state as the initial state, starting from the target Hamiltonian and mixing it with the transverse field Hamiltonian to return to the original target Hamiltonian~\cite{perdomo:sombrero, nishimori:reverse-pspin, nishimori:reverse-pspin-2}. \par

% Reverse annealing has been predicted to have better performance than forward annealing in several systems~\cite{perdomo:sombrero,Chancellor2017,King:2018aa,Ottaviani2018,nishimori:reverse-pspin, nishimori:reverse-pspin-2,Venturelli2019,Marshall2019}.
% Passareli et al. \cite{Passarelli2019}
In chapter.~\ref{chap:ira} (Ref.~\cite{Passarelli2019}) we tested the effectiveness of reverse annealing on the $p$-spin models with $p=3$. Numerical calculations using the time-dependent Lindblad equation~\cite{ABLZ:12-SI, yip:mcwf} showed that the presence of coupling with the boson field significantly changes the performance of reverse annealing. This result shows that the interaction with the environment can lead to an improvement in performance. Since qubits and the environment should always interact with each other on the hardware, it is likely that this result reflects real physical processes in the D-Wave device. However, it is non-trivial if the coupling model  to the environment of chapter.~\ref{chap:ira} is good enough to explain experimental data.

In chapter.~\ref{chap:ira},  the $p$-spin model with $p=3$  is discussed because it has a first-order phase transition in the thermodynamic limit~\cite{derrida:p-spin, Gross1984, Bapst2012} and thus is a difficult problem for simple quantum annealing.
In the present chapter, we instead choose the $p$-spin model with $p=2$. This is because the $p$-spin model with $p=2$ can be more easily embedded in the hardware graph and thus may avoid possible errors by minor embedding of three-body interactions, and also larger system sizes can be implemented up to about 40. Another reason is the existence of ground-state degeneracy, which leads to additional information as will be described below. Our goal is not to study if a high success probability can be achieved but to understand the mechanism of physical processes working in the device.
\par

In this chapter, we perform reverse annealing on the $p$-spin model with $p=2$ experimentally and numerically. We compare the success probability with numerical results to verify how the performance of reverse annealing on D-Wave 2000Q is enhanced by  coupling to baths of boson field.\par

The structure of this chapter is as follows.  We first define the $p$-spin model with $p=2$ in Sec.~\ref{sec:problem} and describe the sampling details in this Chapter. In Sec.~\ref{sec:exp}, we discuss the sampling results from the D-Wave device. We show and examine the simulation results based on closed system Schr\"{o}dinger equation and quantum adiabatic master equation with two system-bath coupling models in Sec.~\ref{sec:simulations}. In Sec.~\ref{sec:classical}, we perform semi-classical simulations of spin-vector Monte Carlo to further understand the experimental data. The chapter is concluded in Sec.~\ref{section:Conc}.

\section{\label{sec:problem}Problem definition and protocols of experiment}
In this section we describe the problem Hamiltonian and the protocol of reverse annealing we used both in our experiments on the D-Wave device and in numerical simulations.
\subsection{\label{sec:p-spin}$p$-spin model with $p=2$}
The total Hamiltonian of quantum annealing is composed of the driver Hamiltonian $H_D$ and the target Hamiltonian $H_0$, the latter encoding the combinatorial optimization problem represented as an Ising model, as a function of dimensionless time $0≤s(t)≤1$,
\begin{equation}
\label{eq:H}
H(s) = \frac{A(s)}{2} H_{D} + \frac{B(s)}{2}H_{0}.
\end{equation}
Here, $A(s)$ and $B(s)$ are device-dependent annealing schedules. The D-Wave device used in the present experiment has $A(s)$ and $B(s)$ as shown in Fig.~\ref{fig:schedule}(a).

\begin{figure}[h!]
\centering
\includegraphics[width=0.85\columnwidth]{ra2/_paper__DW_2000Q_6_schedule-_RAFA_protocol-eps-converted-to.pdf}%
\caption{(a) Annealing schedule of D-Wave 2000Q with the ``DW\_2000Q\_6`` solver. (b) Forward (FA) and reverse annealing (RA) protocols are defined in Eq.~\eqref{eq:FA}) and Eq.~\eqref{eq:RA}), respectively. The latter incorporates a pause of duration $t_p$. The parameter values for this plot are $t_a=1.2~\mathrm{[\mu s]}$ for forward annealing, and $\tau=1~\mathrm{[\mu s]}$ and $s_\mathrm{{inv}}=0.2$ for reverse annealing. These symbols are defined below Eq.~\eqref{eq:RA}.}
\label{fig:schedule}
\end{figure}

The final values of the schedule functions are $A(s=1)≈0$ and $B(s=1)>0$ so that ideally the ground state of $H_0$ is realized as the final state, irrespective of the type of annealing protocol, i.e., the time dependence of $s(t)$ as discussed in the next subsection.
The driver Hamiltonian $H_D$ is usually chosen as
\begin{equation}\label{eq:QA}
\end{equation}
where $N$ is the number of qubits and $σ_i^x$ is the $x$ component of the Pauli matrix acting on the $i$th qubit.

We study the $p$-spin model,
\begin{equation}
\end{equation}
with $p=2$, in contrast to $p=3$ used in chapter.~\ref{chap:ira}. The main reason for the choice of $p=2$ is that we can study the effects of ground-state degeneracy in $p=2$ in contrast to the non-degenerate ground state of $p=3$. Another technical reason is that two-body interactions can be directly embedded on the D-Wave device, which would reduce possible minor errors coming from reduction of three-body interactions to two-body representations.

\subsection{Reverse annealing protocol}

The traditional protocol of forward  annealing has
\begin{equation}
\label{eq:FA}
s(t)=\frac{t}{t_{a}}\ \ \ t\in[0,t_{a}],
\end{equation}
where $t_a$ is the total annealing time (see Fig.~\ref{fig:schedule}(b), the red full line). The initial state of forward annealing is the ground state of the driver Hamiltonian $H_D$.

In reverse annealing, we start from $s=1$ and decrease $s$ to an intermediate value $s_inv$, pause, then increase $s$ to finish the process at $s=1$. The explicit time dependence $s(t)$ realized on the D-Wave device is
\begin{equation}
\label{eq:RA}
1-\,\displaystyle\frac{\,t}{\tau} & 0 \leq t \leq t_{\rm inv} \\
1-\displaystyle\frac{t_{\rm inv}}{\tau} = s_{\mathrm{inv}} & t_{\rm inv}\le t \leq t_{\rm inv}+t_p \\
2s_{\mathrm{inv}} - 1 + \displaystyle\frac{t}{\tau} & t_{\rm inv} +t_p< t \leq 2t_{\rm inv}+t_p,
\end{cases}
\end{equation}
where $1/τ$ is the annealing rate,  $t_inv$ is the inversion time defined as $t_inv=τ(1-s_inv)$, and $t_p$ is the length of an intermediate pause. Our reverse annealing begins from $s(t=0)=1$ and ends at $s(t=t_a)=1$, where
\begin{equation}
t_a = 2t_{\rm{inv}}+t_p = 2 \tau(1-s_{\mathrm{inv}}) + t_p  ,
\label{eq:t_a}
\end{equation}
passing through the inversion point $s_inv$; see the blue dashed line in Fig.~\ref{fig:schedule}(b). To distinguish between $τ$ and $t_a$, we henceforth refer to the former as the ``annealing time" and the latter as the ``total annealing time". The initial state of reverse annealing is a classical state, usually a candidate solution to the combinatorial optimization problem, i.e, a state that is supposed to be close to the solution.

%\textcolor{blue}{[YB: If the initial state is initialized for each reverse annealing process, it is one-time reverse annealing. On the other hand, if the final state obtained from the previous reverse annealing process is repeatedly adopted as the initial state for the next reverse annealing process, this protocol is iterated reverse annealing.]}\par
We can iteratively repeat the process of reverse annealing with the final state of a process fed into the system as the initial condition of the next iteration. We call this protocol iterated reverse annealing and the denote the iteration number as $r$.
In our experiments for single-run reverse annealing, we constructed 15 instances (embeddings on the Chimera graph of the device) and generated 10 random gauges for each set of parameter values, and 1000 annealing runs were repeated for a given random gauge. Thus, the total number of runs for each set of parameter values is 150,000. In iterated reverse annealing, we constructed 15 instances and generated 90 random gauges for each $r$. Thus, the total number of samples at each $r$ and each $s_inv$ is 1350.\par

Notice that the present protocol of reverse annealing on the D-Wave device is different from the one adopted in Chapter.~\ref{chap:ira} ( Ref.~\cite{Passarelli2019}) for numerical computation. In particular, when $s_inv$ is close to 0, the second part to increase $s$ of their protocol has a very sharp increase of $s$ as a function of $t$, which seems to have lead to a significant drop of success probability for $s_inv≈0$. This is not the case in the present protocol, where the time derivative of $s(t)$ is a constant $±1/τ$ (or 0 in pausing) irrespective of $s_inv$.

\section{Experimental results}
\label{sec:exp}

In this section we report on the performance of reverse annealing for the $p=2$ $p$-spin model on the D-Wave 2000Q device.

\subsection{Dependence on initial conditions}
\label{sec:initcond}

Figure~\ref{fig:RA_20_initial} shows the experimental success probability, i.e., the ground-state probability of the final state for $N=20$ and no pause. The initial condition is specified by the initial value of the normalized total magnetization,
\begin{align}
\label{eq:mo}
m_0=\frac{1}{N}\sum_{i=1}^N \langle \psi_0 |\sigma_i^z |\psi_0\rangle,
\end{align}
where $|ψ_0⟩$ denotes the initial wave function (a classical state). Since the doubly degenerate ground state of the problem Hamiltonian $H_0$ has $±1$ as the normalized magnetization, a value of $m_0$ closer to 1 (or $-1$) represents an initial condition exhibiting higher overlap with the ground state. In Fig.~\ref{fig:RA_20_initial} we present results for a completely unbiased initial condition ($m_0=0$) and initial conditions strongly biased toward the all-up state $σ_i^z=1 ∀i$ ($m_0=0.8,0.9$). The state with $m_0=0$ is the highest excited state, and the state with $m_0=0.9$ is the first excited state.

% \subsection*{Dependence on initial condition}
% Figure~\ref{fig:RA_20_initial} shows the success probability, the ground-state probability at the final state. The parameter values are $N=20$, $t_p=0.0~[\mu s]$, and $\tau=5~[\mu s]$. The initial condition is specified by the initial value of normalized total magnetization,
% \begin{align}
% \label{eq:mo}
%     m_0=\frac{1}{N}\sum_{i=1}^N \langle \psi_0 |\sigma_i^z |\psi_0\rangle,
% \end{align}
% where $|\psi_0\rangle$ denotes the initial wave function (a classical state). Since the ground state of the problem Hamiltonian $H_0$ has $\pm 1$ as the normalized magnetization, a closer value of $m_0$ to 1 (or $-1$) represents a better initial condition.

% Figure~\ref{fig:RA_20_initial}(a) shows the (partial) success probability that the system reaches the perfectly up state ($\sigma_i^z=1,~\forall i$, denoted ``up") and the perfectly down state ($\sigma_i^z=-1,~\forall i$, denoted ``down").
% \begin{figure}[tb]
% \includegraphics[width=0.75\columnwidth]{ra2/[paper]_N20_tau5_total_decomposed.eps}%
% \caption{\label{fig:RA_20_initial}Success probability for different initial conditions $m_0$ in reverse annealing on the D-Wave 2000Q device as a function of the inversion point $s_{\mathrm {inv}}$.% The interval of $s_{\mathrm {inv}}$ is 0.02.
% (a) Partial success probability for $m_0=0.0, 0.8$, and $0.9$. (b) Total success probability for the same set of $m_0$ as in (a). The labels of ``up" and ``down" mean the populations of the all-up state and the all-down state, respectively. The dashed line is the minimum-gap point at $s_{\Delta}\approx 0.36$ for $N=20$.}
% \end{figure}
% When the initial state is $m_0=0.0$, the up and down probabilities are mostly equal for all $s_{\rm {inv}}$ as naturally expected because the initial state has no bias toward one of them. It is observed that both success probabilities are close to $0.5$ for $s_{\rm inv}<0.4$ whereas they are zero for $s_{\rm inv}>0.5$.  The vertical dashed line at $s_{\rm inv}\approx 0.36$ indicates the point where the energy gap between the ground and the first excited states becomes minimum, to be denoted as $s_{\Delta}$. It is likely that the system remains close to its initial state for $s_{\rm inv}>0.5$, where quantum fluctuations are small, and the true final ground state is hard to be reached after the whole process since the initial state with $m_0=0$ is far away from the ideal state with $m_0=\pm 1$. The situation is different if $s_{\rm inv}<0.4$, i.e., if the system goes across the minimum-gap region. The effects of the initial condition disappears and the system becomes quantum-disordered (paramagnetic) by large quantum fluctuations caused by the driver $H_D$. This renders the process similar to the conventional forward quantum annealing, by which two true ground states are reached with almost the same probability 0.5. Notice that the phase transition around $s_{\Delta}$ is of second order in the present problem, and therefore the system finds the true final ground state relatively easily by forward annealing because of the mild, polynomial, closing of the energy gap.

% For the initial states with $m_0=0.8$ and $m_0=0.9$, differences show up between probabilities of the final all-up and all-down states, which reflects the bias of the initial state toward the former (all up) even after the minimum gap region is crossed ($s_{\rm inv}<s_{\Delta}=0.36$). It is remarkable that the success probability for the all-up state is almost 1 around $s_{\rm inv}=0.36$. This may come from relaxation of the system to the ground state as a result of coupling to the environment when the energy gap becomes smallest but quantum fluctuations are not large enough to erase the initial condition.  It is also noticed in Fig.~\ref{fig:RA_20_initial}(a) that the initial condition affects the system properties (i.e., different probabilities for all-up and all-down states) even in the paramagnetic region $s_{\rm inv}<0.36$ where quantum fluctuations are large and all-up and all-down states are expected to have the same probability in equilibrium. The system ``remembers" the initial condition to a certain extent.

% Figure~\ref{fig:RA_20_initial}(b) shows the total success probability, the sum of final up and down probabilities. The resulting curve stays almost flat close to 1 for $s_{\rm inv}<0.4$.  This constant 1 should be a mixture of two effects manifest in Fig.~\ref{fig:RA_20_initial}(a), the genuine effects of reverse annealing around $s_{\rm inv}\approx 0.4$ for $m_0=0.9$ and 0.8 and the effectively-forward-annealing-like behavior for $m_0=0.0$.

% This result in Fig.~\ref{fig:RA_20_initial}(b) is consistent with the tendency found numerically in Ref.~\cite{Passarelli2019} for $p=3$ except for very small values of success probability near $s_{\rm inv}=0$ in their data. This may be explained by the very quick forward annealing in this parameter region in their protocol, keeping the system almost unchanged from the quantum-disordered state at $s=s_{\rm inv}$.
% Aside from this subtlety, our experimental data are consistent with their result, supporting their expectation that the system-environment interactions through spin-boson coupling prompts relaxation of the system toward the ground state around the minimum-gap region.

\subsubsection{$m_0=0$}

Figure~\ref{fig:RA_20_initial}(a) shows the probability that the system reaches the all-up state (denoted ``up") and the all-down state ($σ_i^z=-1, ∀i$, denoted ``down"). We call these the partial success probabilities.

\begin{figure}[h!]
\subfigure[]{\includegraphics[width = 0.4939\columnwidth]{ra2/image_part_001.png}}
\subfigure[]{\includegraphics[width = 0.4939\columnwidth]{ra2/image_part_002.png}}
\caption{Success probability for different initial conditions $m_0$ (magnetization) in reverse annealing on the D-Wave 2000Q device as a function of the inversion point $s_{\mathrm {inv}}$. The dashed line is the minimum-gap point at $s_{\Delta}\approx 0.36$ for $N=20$. When $s_{\mathrm {inv}} < s_{\Delta}$ the reverse direction of the anneal goes through and past the minimum gap point, and then crosses it again during the forward anneal. When $s_{\mathrm {inv}} > s_{\Delta}$ there is no crossing of the minimum gap point.
(a) Partial success probability for $m_0=0, 0.8$, and $0.9$. (b) Total success probability for the same set of $m_0$ as in (a). Here and in all subsequent figures the labels ``up" and ``down" mean the populations of the all-up state and the all-down state, respectively. Error bars denote standard deviations.}
\label{fig:RA_20_initial}
\end{figure}

When the initial condition is $m_0=0$, the up and down probabilities are equal to within the error bars for all $s_inv$, as expected because the initial state is unbiased. Note that both success probabilities are close to $0.5$ for $s_inv<0.4$ whereas they are zero for $s_inv>0.5$.  The vertical dashed line at $s_inv≈0.36$ indicates the minimum-gap point, i.e., the point where the energy gap between the ground and the second excited state becomes minimum, which we denote by $s_Δ$ (the first excited state becomes degenerate with the ground state at the end of the anneal, hence it is the second excited state that is relevant). It is likely that the system remains close to its initial state for $s_inv>0.5$, where the transverse field and the associated quantum fluctuations are small, and the true final ground state with $m_0=±1$ is hard to reach, since the initial state with $m_0=0$ has small overlap with the latter. The situation is different when $s_inv<0.4$, i.e., when the system traverses the minimum-gap region. In this case the system enters the paramagnetic (quantum-disordered) phase dominated by the driver $H_D$, so that the initial condition is effectively erased. This renders the process similar to conventional forward quantum annealing, by which the two true ground states are reached with nearly equal probability $0.5$. Recall that the phase transition around $s_Δ$ is of second order in the present problem, and therefore the system finds the true final ground state relatively easily by forward annealing because of the mild, polynomial, closing of the energy gap, as we discuss in more detail later (see Fig.~\ref{fig:RA_20_size}(c)).
% The overall tendency is similar to Fig.~\ref{fig:RA_20_initial}.

\subsubsection{$m_0=0.8,0.9$: up/down symmetry breaking for $0.3 \lesssim s_{\mathrm{inv}} \lesssim 0.5$}

For the initial states with $m_0=0.8$ and $m_0=0.9$, the experimental results shown in Fig.~\ref{fig:RA_20_initial}(a) reveal significant differences between the probabilities of the final all-up and all-down states. Of course, the initial conditions $m_0=0.8,0.9$ have much larger overlap with the all-up state, and this fact alone suggest a mechanism for breaking the symmetry between the probabilities of the final state being the all-up or all-down state. However,
%there remains a bias toward the all-up state even after the minimum gap region is crossed ($s_{\mathrm{inv}}<s_{\Delta}=0.36$), and
it is remarkable that the success probability for the all-up state is almost 1 in a region $0.3 ≲s_inv≲0.5 $ both to the left and to the right of $s_Δ≈0.36$ (we discuss the additional asymmetry for $s_inv ≲0.3$ below). As we explain in Sec.~\ref{sec:AME}, this behavior is inconsistent with a model of an open quantum system that is weakly coupled to its environment and is thus described by the adiabatic master equation~\cite{ABLZ:12-SI}. However, it is consistent with a simple semiclassical model captured by the spin-vector Monte Carlo algorithm~\cite{shin2014quantum,Albash2020ComparingRM}, as explained in Sec.~\ref{sec:PSP-SVMC}.
%This may be attributed to relaxation of the system to the ground state as a result of coupling to the environment when the energy gap becomes small but quantum fluctuations are not large enough to erase the initial condition. We further discuss the explanation for this effect in Sec.~\ref{sec:simulations}.

\subsubsection{$s_{\mathrm{inv}} \gtrsim 0.5 > s_{\Delta}$: freezing}
\label{sec:freezeout}

Note that for all three initial conditions the success probability eventually vanishes for $s_inv ≳0.5 > s_Δ$. The reason is that in all cases the system is initialized in an excited state, and remains in an excited state at the end of the anneal, since there is no mechanism for thermal relaxation to a lower energy state when $s_inv$ is large. This is a manifestation of the phenomenon of freezing~\cite{ABLZ:12-SI,amin:freeze-out}, i.e., the extreme slowdown of relaxation due to the fact that the system-bath interaction (nearly) commutes with the system Hamiltonian when the transverse field is very small. In addition, the annealing timescale is manifestly too slow for downward diabatic transitions. This is true even in the absence of pausing, which suggests a discontinuity in the derivative of the schedule depicted in Fig.~\ref{fig:schedule}(b). Despite this, the reversal of the anneal direction is apparently too slow in practice to have a non-adiabatic effect, or if diabatic transitions do occur, then they exclusively populate higher excited states.

\subsubsection{$s_{\mathrm{inv}} \lesssim 0.3 < s_{\Delta}$: spin-bath polarization}
\label{sec:spinbathpol}

It is also noteworthy from Fig.~\ref{fig:RA_20_initial}(a) that the initial condition results in different probabilities for the all-up and all-down states even in the paramagnetic region $s_inv ≲0.3 < s_Δ$, where quantum fluctuations are large and the all-up and all-down states are expected to have the same probability in equilibrium. The system ``remembers" the initial condition to a certain extent, an anomaly that may be attributable to spin-bath polarization~\cite{spin-bath-polarization}. Namely, the persistent current flowing in the qubit body during the anneal produces a magnetic field that can partially align or polarize an ensemble of environmental spins local to the qubit wiring, with a much slower relaxation time than the anneal duration. Given the polarized initial condition $m_0=0.8$ or $0.9$, this polarized spin-bath will be aligned with the all-up state even after the system crosses into the paramagnetic phase, thus preventing the system from equilibrating and explaining the observed memory effect.

\subsubsection{Total success probability}

Figure~\ref{fig:RA_20_initial}(b) shows the total success probability, the sum of the final up and down probabilities. The resulting curve stays almost flat and close to $1$ for $s_inv<0.4$.  This constant $1$ is a mixture of the two effects manifest in Fig.~\ref{fig:RA_20_initial}(a), the genuine effects of reverse annealing around $s_inv≈0.4$ for $m_0=0.9$ and 0.8 and the effectively-forward-annealing-like behavior for $m_0=0$. We shall see in Sec.~\ref{sec:simulations-init} that all the salient features seen in Fig.~\ref{fig:RA_20_initial}(b), such as the shift to the left with decreasing $m_0$, are captured by our open system simulations.

The results seen in Fig.~\ref{fig:RA_20_initial}(b) are also consistent with the trend found numerically in Ref.~\cite{Passarelli2019} for the $p=3$ $p$-spin model except for the very small values of success probability near $s_inv=0$ seen in the latter's numerical simulations. This may be explained by the very fast forward annealing in this parameter region in Ref.~\cite{Passarelli2019}, which keeps the system almost unchanged from the quantum-disordered state at $s=s_inv$.
Aside from this subtlety, our experimental data are consistent with the result of Ref.~\cite{Passarelli2019}, supporting the latter's expectation that the system-environment interaction (through spin-boson coupling) prompts relaxation of the system toward the ground state around the minimum-gap region.

\subsection{Dependence on annealing time}
As seen from Eq.~\eqref{eq:t_a}, the total annealing time $t_a$ is linearly dependent on the annealing time $τ$, and proportional to $τ$ when $t_p=0$.
Figure~\ref{fig:RA_20_diff_tau} shows the success probability for different $τ$ with $N=20$, $m_0=0.9$, and no pausing.  The overall trend is similar to Fig.~\ref{fig:RA_20_initial}.

As shown in Fig.~\ref{fig:RA_20_diff_tau}(a), an increase of $τ$ leads to slightly higher all-up success probabilities for $s_inv> s_Δ$ but the other way around for $s_inv < s_Δ$. This can be explained in terms  of an increased relaxation to the ferromagnetic ground state near the minimum energy gap for larger $τ$ and an enhanced relaxation to the paramagnetic ground state for $s_inv<s_Δ$.

\begin{figure}[h!]
\subfigure[]{\includegraphics[width = 0.4939\columnwidth]{ra2/image_part_003.png}}
\subfigure[]{\includegraphics[width = 0.4939\columnwidth]{ra2/image_part_004.png}}
\caption{Success probability for different annealing times $\tau$ in reverse annealing on the D-Wave 2000Q device as a function of $s_{\mathrm{inv}}$.
%The interval of $s_{\mathrm{inv}}$ is 02.
(a) Partial success probability for $\tau=5, 10$, and $20[\mu s]$. (b) Total success probability for the same set of $\tau$ as in (a). The dashed line is the minimum-gap point at $N=20$, $s_{\Delta}\approx 0.36$.}
\label{fig:RA_20_diff_tau}
\end{figure}

The total success probability stays close to $1$ for $s_inv<s_Δ$ and benefits slightly from increasing $τ$ for $s_inv>s_Δ$ (Fig.~\ref{fig:RA_20_diff_tau}(b)).

\subsection{Effects of pausing}

Figure~\ref{fig:RA_20_pause} shows success probabilities for pausing time $t_p∈{0, 2, 10} [μs]$ as a function of $s_inv$  with $N=20$, $m_0=0.9$ and $τ=5 [μs]$.

\begin{figure}[h!]
\subfigure[]{\includegraphics[width = 0.4939\columnwidth]{ra2/pause_image_part_001.jpg}}
\subfigure[]{\includegraphics[width = 0.4939\columnwidth]{ra2/pause_image_part_002.jpg}}
\caption{Success probability for different pausing times $t_p$ in reverse annealing as a function of $s_{\mathrm{inv}}$.
%The interval of $s_{\mathrm{inv}}$ is 02.
(a) Partial success probability for pause $t_p=0, 2$, and $10~[\mu s]$. (b) Total success probability for the same set of $t_p$ as in (a). The dashed line is the minimum-gap point at $N=20$, $s_{\Delta}\sim0.36$.}
\label{fig:RA_20_pause}
\end{figure}

The overall trend is similar to Fig.~\ref{fig:RA_20_diff_tau}, i.e., a longer pause leads to an increased relaxation, but the effects are more significant in the present case. Pausing at $0.5<s_inv<0.6$ greatly improves the success probability from nearly 0 ($t_p= 0 [μs]$) to nearly 1 ($t_p= 2 [μs]$ and $10 [μs]$), implying that pausing in a relatively narrow region slightly past the minimum gap  point is most effective. This is in line with previous findings~\cite{marshall,chen2020pausing,Albash2020ComparingRM}.

Another significant difference between pausing and not pausing is that the partial success probabilities shown in Fig.~\ref{fig:RA_20_pause}(a) are flat for $t_p>0$, showing that the effect responsible for the anomaly disappears under pausing. This is consistent with the spin-bath polarization effect discussed in Sec.~\ref{sec:initcond}, in that the pause provides the time needed for this polarization to relax to equilibrium, on a timescale of a few $μs$.

\begin{figure}[h!]
\subfigure[]{\includegraphics[width = 0.4939\columnwidth]{ra2/differentN_image_part_001.jpg}}
\subfigure[]{\includegraphics[width = 0.4939\columnwidth]{ra2/differentN_image_part_002.jpg}}
\subfigure[]{\includegraphics[width=0.4739\linewidth]{ra2/gaprevealplot_3.pdf}}
\caption{Success probability for different system sizes $N$ in reverse annealing  as a function of $s_{\mathrm{inv}}$.
(a) Partial success probability for system size $N=4, 8, 12, 16, 20$, and $24$. (b) Total success probability for the same set of $N$ as in (a). (c) The position of the minimum gap $\Delta$ for $4 \le N \le 22$, along with the value of $s_{\mathrm{inv}}$ at which the success probability $=.5$ for each $N$. Inset: the minimum gap for each $N$.}
\label{fig:RA_20_size}
\end{figure}

\subsection{Size dependence}

% \DL{About Fig.~\ref{fig:RA_20_size}: Bottom panel was an appendix figure. Let's have this panel also display the value of $s_{\mathrm{inv}}$ at which the success probability $=.5$ for each $N$; let's plot the gap as the inset (as a log-log plot to highlight the poly($N$)-dependence) and $s_\Delta$ as the main plot. Let's also replot this panel so the font and size are consistent with the two panels above it. It needs to be labeled (c).}

Figures~\ref{fig:RA_20_size}(a) and (b), respectively, show the partial and total success probability for different system sizes $N$ as a function of $s_inv$ with $τ=5 [μs]$. The initial state for each $N$ is the first excited state with $m_0<1.0$ (closest to the ground state for which $m_0=1.0$), i.e., a state with one spin flipped. Although we observe a slight non-monotonicity with $N$ in the region $s_inv≥s_Δ$ in Fig.~\ref{fig:RA_20_size}(a), the overall trend is clear, namely, the drop-off to zero success probability occurs at smaller $s_inv$ as $N$ is increased. This is consistent with the reduction in $s_Δ$ as a function of $N$, shown in Fig.~\ref{fig:RA_20_size}(c), which is tracked by the value of $s_inv$ at which the success probability equals $0.5$.

Additionally, the asymmetry between the all-up and all-down states for $s_inv<s_Δ$ is enhanced with increasing $N$. This is consistent with the decreasing gap seen in Fig.~\ref{fig:RA_20_size}(c), which should lead to slower relaxation in thermal equilibrium: at inverse temperature $β$, the probability of an excitation is $exp(-βΔ)$ times the probability of relaxation between the same pair of states, separated by a gap $Δ$. I.e., the asymmetry induced by the spin-bath polarization effect is expected to take longer to dissipate for larger $N$ due to the slower relaxation of the system.

Finally, the overlap of data for $N=20$ and $24$
%around $s_{\mathrm{inv}}=0.5$
suggests that large-$N$ effects have already converged at these sizes, i.e., that $N∼20$ is  sufficiently large to infer the behavior at larger $N$. Also, even the smallest systems with $N=4$ and $8$ already share qualitative features with larger systems.

\subsection{Effects of iteration}

We next study the effects of iteration on reverse annealing, i.e., how the success probability depends on the number of iterations $r$. Figure~\ref{fig:IRA} shows the result for $r=1, 10, 25$ and $50$ at $s_inv$ with $N=20$, $m_0=0.9$, $t_p=0 [μs]$, and $τ=5 [μs]$.
In Fig.~\ref{fig:IRA}, the success probabilities improve in the region $s_inv≥s_Δ$ as the number of iterations $r$ increases, regardless of the initial state $m_0$. The relaxation to the low-energy state due to coupling to the environment should be successively induced by iterations, and the occupation of the ground state should increase. The results in Fig.~\ref{fig:IRA} confirm this expectation.

\begin{figure}[h!]
\centering
\includegraphics[width=0.95\columnwidth]{ra2/_paper__reinitialize_N20_tau5-eps-converted-to.pdf}%
\caption{Success probability for different number of iterations $r$ in iterated reverse annealing  as a function of $s_{\mathrm{inv}}$. (a) and (b) are partial and total success probabilities for number of iterations $r=1, 10, 25$, and $50$ with the initial state $m_0=0.9$ (the first excited state), respectively. (c) and (d) are partial and total success probabilities with the initial state $m_0=0$ (the highest excited state). The dashed line is the minimum-gap point at $N=20$, $s_{\Delta}\sim0.36$.}
\label{fig:IRA}
\end{figure}

Comparing Figs.~\ref{fig:IRA}(a) and (b) with $m_0=0.9$ and Figs.~\ref{fig:IRA}(c) and (d) with $m_0=0$, we can see that the initial state with $m_0=0.9$ has a larger improvement in success probability with fewer iterations $r$ in the region where $s_inv ≥s_Δ$. The initial state with $m_0 = 0$ deviates greatly from the ground state, and there are many excited states between it and the ground state. Therefore, when the system state transitions from the initial state of $m_0 = 0$ to other states by relaxation due to coupling to the environment, those excited states become candidates. Many iterations are expected to be necessary to obtain a significant improvement in the success probability. On the other hand, there are only a few excited states between the initial state $m_0≲1.0$ and the ground state. Therefore, it is expected that only a few iterations will be sufficient to make the transition to the ground state. Our results support this picture of improved performance under iterated reverse annealing.

\section{Numerical simulations}
\label{sec:simulations}

We next present closed and open-systems simulations and compare the results with the data presented in the previous section. We choose relatively small system sizes $N=4$ and $8$ to facilitate numerical computations. We note that while we can thus simulate the initial condition $m_0=0$ (equal number of spins up and down), we are limited to $m_0=0.5$ for $N=4$ and $m_0=0.75$ for $N=8$ (a single spin down). However, our experimental data do not show a strong size dependence, as discussed in the previous section, so the qualitative comparison our numerical results will enable should suffice to draw relevant physical conclusions.

For simplicity, we focus on reverse annealing without pausing (i.e., $t_p=0$) in this section. We expect that in general pausing will lead to an overall success probability  in the results of open-system simulation.

\subsection{Closed system model}
An analysis of the closed system case, while being unrealistic due to the absence of thermal effects, is instrumental in isolating the effect and importance of diabatic transitions in explaining our experimental results.

The time evolution of reverse annealing in a closed system is described as follows.
The state after a single iteration (cycle) is
\begin{equation}
\ket{\psi(2t_{\text{inv}})} = U(2t_{\text{inv}},0)\ket{\psi(0)}\,,
\end{equation}
where $|ψ(0)⟩$ is the initial state and
\begin{equation}
U(2t_{\text{inv}}, 0) = {\cal T}\exp\left[-i\int_0^{2t_{\text{inv}}}H(t')dt'\right]
\end{equation}
is the unitary time-evolution operator and $T$ denotes forward time ordering.
At the end of $r$ cycles, the final state $|ψ(2rt_inv)⟩$ can be expressed as:
\begin{equation}
\ket{\psi(2rt_{\text{inv}})} = U(2rt_{\text{inv}}, 0)\ket{\psi(0)}\,,
\end{equation}
\begin{equation}
U(2rt_{\text{inv}}, 0) = \prod_{i=0}^{r-1}U(2(i+1)t_{\text{inv}}, 2it_{\text{inv}}) \,.
\end{equation}
Here, the final state of the $r$th cycle is the initial state of the next cycle. This condition is shared by the experiments in the previous section. The solution states are doubly degenerate, and the total success probability at the end of $r$ cycles is
\begin{equation}
p(r)=|\langle\psi(2rt_{\text{inv}})\ket{\text{up}}|^2 + |\langle\psi(2rt_{\text{inv}})|\ket{\text{down}}|^2 \,,
\end{equation}
where $|up⟩ = |↑⟩^⊗N$ and $|down⟩ = |↓⟩^⊗N$.
For higher efficiency without loss of accuracy, we rotate the state vector into the instantaneous energy eigenbasis representation at each time step in our numerical simulations.

\subsubsection{Dependence on system size and annealing time}

We initialize the state to have a single spin down, i.e., as the computational basis state $|0001⟩$ ($|0⟩ ≡|↑⟩$, $|1⟩ ≡|↓⟩ (m_0=0.5)$) for $N=4$ and $|0000001⟩ (m_0=0.75)$ for $N=8$, respectively.
Note that in our simulations these are not exact eigenstates of $H(1)$, due to a small residual transverse field at $s=1$, as in the D-Wave annealing schedule shown in Fig.~\ref{fig:schedule}.
Note further that our chosen computational basis initial states are not eigenstates of $S^2 = ∑_α∈{x,y,z}(S^α)^2$, where $S^α=1/2∑_i=1^Nσ_i^α$. Therefore the computational basis states do not lie in the subspace of any fixed value of $S$, where $S$ represents the total spin quantum number of $N$ qubits.