Quantum annealing with twisted fields

Quantum annealing (QA) [1–5], which has attracted considerable attention in recent decades, was originally used as a method for solving combinatorial optimization problems. It is well known that such problems can be mapped to finding the ground states of Ising Hamiltonians [6–8]. It is known that the problem to obtain the ground state of the Ising spin glass models is NP-hard [9]. In QA, we prepare a ground state of a Hamiltonian of a transverse field, which is called a driving Hamiltonian; then, we let the system evolve with a time-dependent Hamiltonian to change from the transverse field to the problem Hamiltonian. As long as an adiabatic condition is satisfied and the Hamiltonian contains matrix elements to induce transitions between the initial state and the target state, the ground state of the problem Hamiltonian can be prepared with high probability after the dynamics. By measuring the ground state, we can obtain the solutions of combinatorial optimization problems [10, 11].

By choosing a special driving Hamiltonian and problem Hamiltonian, we can simulate an arbitrary quantum circuit by an adiabatic evolution, and this is called an adiabatic quantum computation [12]. However, to realize this, many ancillary qubits are required to generate a complicated state called a history state, which is in stark contrast to the original QA [1–5]. The focus of the present paper is the original one that does not require either such ancillary qubits or a history state.

Furthermore, many studies have attempted to use QA for quantum chemical calculations. The Hamiltonian of molecules is written in the second quantized form. There are methods for transforming the Hamiltonian of the second quantized form into that of qubit forms [13–18]. In quantum chemistry, high accuracy is required in the calculation of molecular energy, which is referred to as chemical accuracy. Knowledge of energy with chemical accuracy is essential for predicting the chemical reactions of molecules [19].

There are several other applications of QA. For example, QA has been applied to clustering [7, 20]. Further, a method that uses QA to perform calculations for topological data analysis has been reported [21]. In addition, QA has been investigated for solving the shortest-vector problem, which is a candidate for postquantum cryptography [22].

D-Wave Systems Inc. has realized QA machines composed of thousands of qubits [23]. Superconducting flux qubits have been employed in these machines. Many experimental demonstrations of QA have been performed using these machines, including machine learning, and graph coloring problems [24–27].

Although considerable effort has been devoted toward identifying useful applications of QA, many obstacles are yet to be overcome for using QA to solve practical problems. In particular, there are two major obstacles: decoherence and non-adiabatic transitions [28]. To suppress decoherence, we need to implement QA with a shorter time schedule. However, as the annealing time becomes shorter, more non-adiabatic transitions will occur during QA. This trade-off makes it difficult to use QA for solving practical problems.

A Hamiltonian is called Ising type when it consists of σ^(z), which is similar to the classical Ising model. When the problem Hamiltonian is of the Ising type, we randomly obtain the ground state of the problem Hamiltonian after QA by performing measurements in the computational basis. The density matrix after QA is expressed as ρ = ∑_j p_j |E_j ⟩⟨E_j |, where p_j denotes the population and |E_j ⟩ denotes an eigenvector of the problem Hamiltonian. In this case, the probability of obtaining the ground state depends on the population of the ground state. Thus, if the ground-state population is finite, we can obtain the ground state by increasing the number of trials.

By contrast, if we choose a problem Hamiltonian with off-diagonal matrix elements as the problem Hamiltonian, we cannot estimate the energy of the Hamiltonian by performing measurements in the computational basis. In this case, we need to estimate the energy of the Hamiltonian by performing the Pauli measurements. We consider the case in which the Hamiltonian is composed of the summation of the products of the Pauli matrices, such as H = ∑_i c_i σ_i , where σ_i is the Pauli product and c_i is the coefficient. We obtain the expectation value of each term by performing measurements of the Pauli products in the quantum states after QA, and we take the sum of the expectation values of all the terms. The energy of the Hamiltonian for a density matrix ρ is given by ⟨H⟩ = Tr [Hρ] = ∑_n p_n E_n , and we obtain the expectation values of the Hamiltonian ⟨H⟩ by performing a large number of measurements. Hence, if there is a population in the excited state, the energy measured experimentally is different from the ground-state energy. Therefore, we need to prepare a density matrix close to the exact ground state in order to determine the ground-state energy with high accuracy.

In this paper, we propose a method for suppressing both decoherence and non-adiabatic transitions by using inhomogeneous twist operators that change the angles of the transverse fields during QA. We define inhomogeneous twist operators that rotate the direction of the transverse fields of the driving Hamiltonian, and we also define the twist parameters that correspond to the rotation angle at each qubit. Further, we apply these operators to the driving Hamiltonian. Our method is especially useful especially when the problem Hamiltonian contains off-diagonal terms.

By using this twisted driving Hamiltonian, we can implement QA for the given twist parameters and measure the energy of the state after QA. To minimize the energy, we update the parameters by using, e.g., gradient descent methods, and we perform QA again with different twist parameters. By repeating these processes, we can obtain the ground-state energy, which is lower than that obtained using conventional QA. Through numerical simulations, we demonstrate that our approach suppresses the effect of decoherence and non-adiabatic transitions in QA for some problem Hamiltonians.

The remainder of this paper is organized as follows. Sections 2 and 3 review QA and gradient descent, respectively. Section 4 introduces our scheme with twist operators. Section 5 describes numerical simulations conducted to evaluate the performance of our scheme and shows that, for some problem Hamiltonians, the ground-state energy obtained using our scheme is more accurate than that obtained using the conventional scheme. Finally, section 6 concludes the paper.

Here, we review QA for the ground-state search [3]. We choose the driving Hamiltonian H_D as the transverse field (i.e., ${H}_{\mathrm{D}}=-{\sum }_{i=1}^{L}{\hat{\sigma }}_{i}^{x}$). The total Hamiltonian is described as follows:

$$\begin{equation}H(t)=\left(1-\frac{t}{T}\right){H}_{\mathrm{D}}+\frac{t}{T}{H}_{\mathrm{P}},\end{equation} \tag{ 1 }$$

where T is the annealing time, H_D is the driving Hamiltonian, and H_P is the problem Hamiltonian. We prepare the ground state of the driving Hamiltonian. The driving Hamiltonian is adiabatically changed into the problem Hamiltonian. If the dynamics is adiabatic, the adiabatic theorem guarantees that we can obtain the ground state of the problem Hamiltonian with high probability.

The accuracy of QA is degraded by many types of noise. The relevant types are environmental decoherence and non-adiabatic transitions [28–35]. For the suppression of non-adiabatic transitions, it is necessary to implement QA with a longer time schedule. However, as the annealing time becomes longer, the decoherence effects become more significant during QA. Owing to this trade-off, it is not straightforward to solve practical problems using QA.

Several methods for suppressing decoherence and non-adiabatic transitions have been investigated. Susa et al proposed a way to accelerate the annealing process using an inhomogeneous driving Hamiltonian for a specific case [36, 37]. It is known that 'non-stoquastic' Hamiltonians with negative off-diagonal matrix elements improve the performance of QA for some specific problem Hamiltonians [38–41]. However, it is worth mentioning that the non-stoquastic Hamiltonian does not always increases the energy gap for general problem Hamiltonians [42, 43].

Direct estimation of the energy gap between the ground state and the first excited state by QA has been proposed [44, 45], and this method is robust against non-adiabatic transitions. In addition, considerable effort has been devoted toward suppressing the effect of environmental noise. Error correction of QA has been investigated to suppress decoherence [46–49]. In addition, the idea of using a decoherence-free subspace for QA has been proposed [33, 50]. Spin lock techniques can be adopted to use long-lived qubits for QA [51–53]. Furthermore, several methods using non-adiabatic transitions and quenching for efficient QA have been studied [54–60]. Other approaches have also been proposed to suppress non-adiabatic transitions and decoherence by using variational methods [61–63]. It is worth noting that such variational methods have been adopted to find a ground state of the Hamiltonian by using variational algorithms with near-term intermediate-scale quantum devices [64, 65].

Here, we review the basic gradient descent method. The gradient descent method is a method that searches for the lowest value of the cost function by the gradient. It consists of four steps. First, we determine the learning rate that indicates by how much we can change the parameters when we update them. Second, we derive the gradient of the cost function. When we cannot obtain the gradient of the cost function analytically, we need to use numerical differentiation. Third, we update the parameters using the gradient and learning rate as follows:

$$\begin{equation}{a}^{(1)}={\left.{a}^{(0)}-\alpha \frac{\partial f(a)}{\partial a}\right\vert }_{a={a}^{(0)}},\end{equation} \tag{ 2 }$$

where f is the cost function, α is the learning rate, and a⁽⁰⁾ and a⁽¹⁾ are the initial and updated parameters, respectively. Fourth, we repeat the second and third steps until the cost function converges to a specific value. It is worth mentioning that the gradient descent does not always provide a global minimum but may provide a local minimum of the cost function.

Here, we introduce our scheme for using twist operators with QA. By deforming the driving Hamiltonian with an inhomogeneous twist operator, we aim to obtain the ground-state energy with higher accuracy than that obtained using conventional QA. The inhomogeneous twist operator consists of several parameters, which we refer to as twist parameters. To find the optimal twist parameters for an efficient ground-state search, we use the so-called variational methods, where we adaptively update the parameters according to the measurement results after QA.

4.1. Inhomogeneous twist operator

Let us consider the problem Hamiltonian composed of L qubits. Let ${\sigma }_{j}^{x}$, ${\sigma }_{j}^{y}$, and ${\sigma }_{j}^{z}$ denote the standard Pauli matrices at the jth site. Then, the inhomogeneous twist operator is defined by

$$\begin{equation}{U}_{\text{twist}}({\theta }_{1},\dots ,{\theta }_{L}){:=}\prod\limits _{j=1}^{L}\mathrm{exp}[\mathrm{i}{\theta }_{j}{\sigma }_{j}^{y}],\end{equation} \tag{ 3 }$$

where L is the number of qubits and ${\left\{{\theta }_{j}\right\}}_{j=1}^{L}$ are the twist parameters. We define a single rotational operator of the y-axis at the jth site as

$$\begin{equation}{U}_{j}^{y}{:=}\mathrm{exp}[\mathrm{i}{\theta }_{j}{\sigma }_{j}^{y}].\end{equation} \tag{ 4 }$$

In addition, we define the inhomogeneous twist operator as

$$\begin{align}\hfill {U}_{\text{twist}}({\theta }_{1},\dots ,{\theta }_{L}){:=}& \prod\limits _{j=1}^{L}{U}_{j}^{y}\hfill \\ \hfill & =\mathrm{exp}\left[\mathrm{i}\sum\limits _{j=1}^{L}{\theta }_{j}{\sigma }_{j}^{y}\right].\hfill \end{align} \tag{ 5 }$$

Further, we deform the driving Hamiltonian with the inhomogeneous twist operator. The key idea of our scheme is to use the deformed driving Hamiltonian for QA as follows:

$$\begin{equation}{H}^{(\mathrm{t}\mathrm{w}\mathrm{i}\mathrm{s}\mathrm{t})}(t,{\theta }_{1},\dots ,{\theta }_{L})=\left(1-\frac{t}{T}\right){U}_{\text{twist}}^{-1}({\theta }_{1},\dots ,{\theta }_{L}){H}_{\mathrm{D}}{U}_{\text{twist}}({\theta }_{1},\dots ,{\theta }_{L})+\frac{t}{T}{H}_{\mathrm{P}}.\end{equation} \tag{ 6 }$$

The energy spectrum of this total Hamiltonian is changed by the twist parameters. When we apply the twist operator, the ground state of the driver Hamiltonian is changed from a state pointing along x direction in the Bloch sphere to a state pointing along a direction rotated by θ_j in the x–y plane. Therefore, we can realize the inhomogeneous twist operator in this paper using the transverse magnetic field and the longitudinal magnetic field. In our scheme, we need an independent control of the transverse magnetic field and the longitudinal magnetic field. The current QA machine [23] can generate both the transverse magnetic field and the longitudinal magnetic field. However, it is still difficult to realize independent control of the transverse and the longitudinal magnetic field in the current device. On the other hand, recently, a new architecture of QA to use a capacitively shunted flux qubit was suggested [53, 66], and it is possible to control the transverse and the longitudinal magnetic field independently by adjusting the amplitude and detuning of the microwave fields. So our proposal will be realized by using such an architecture.

4.2. Numerical differentiation

In this subsection, we introduce the numerical differentiation in detail. The numerical differentiation is a method of approximate differentiation when the form of a specific function is not known. Several variants of numerical differentiation exist. Here, we use the following formalism. We define the numerical differentiation of the function f by x as

$$\begin{equation}\frac{\mathrm{d}f(x)}{\mathrm{d}x}\enspace : =\frac{f(x+{\epsilon})-f(x-{\epsilon})}{2{\epsilon}},\end{equation} \tag{ 7 }$$

where is a small constant. If the variance increase, this definition can be used in the same way when there is more than one variable.

4.3. Variational QA with gradient descent

Here, we explain how to apply gradient descent with our scheme. Let E^(ann)(T, θ₁, ..., θ_L ) denote the energy that we measure after QA using the deformed annealing Hamiltonian (6). Our scheme consists of three steps. First, to obtain the derivative of E^(ann)(T, θ₁, ..., θ_L ) with respect to ${\left\{{\theta }_{j}\right\}}_{j=1}^{L}$ for a given T, we set the inhomogeneous twist operators with some twist parameters, perform QA with the deformed Hamiltonian, and measure the energy of the state after QA. It is worth noting that we cannot analytically obtain the derivative of E^(ann)(T, θ₁, ..., θ_L ); hence, we use numerical differentiation. Second, we update the twist parameters on the basis of the results of step 1. Third, we repeat steps 1 and 2 until the energy converges to a finite value. The entire procedure is summarized in algorithm 1. Let us define N as the number of measurements to know the energy of the state after the QA. If the number of measurements required, our protocol requires 2KLN measurements where K denotes the number of steps to update the parameters. The number of measurements required for our protocol is given as

$$\begin{equation}N\propto \frac{KL}{{{\epsilon}}^{2}},\end{equation} \tag{ 8 }$$

where is a constant that appears in the numerical differentiation. We explain the details in appendix B.

Algorithm 1. Calculate optimal ${\left\{{\theta }_{j}\right\}}_{j=1}^{L}$ using gradient descent.

1: Initial states ${\left\{{\theta }_{j}\right\}}_{j=1}^{L}$

2: for k(step number) do

3:  for j(site number) do

4:   $\frac{\partial {E}^{(\mathrm{a}\mathrm{n}\mathrm{n})}({\theta }_{1},\dots ,{\theta }_{L})}{\partial {\theta }_{j}}{\leftarrow}=\frac{{E}^{(\mathrm{a}\mathrm{n}\mathrm{n})}({\theta }_{1},\dots {\theta }_{j}-{\epsilon},\dots ,{\theta }_{L})-{E}^{(\mathrm{a}\mathrm{n}\mathrm{n})}({\theta }_{1},\dots {\theta }_{j}+{\epsilon},\dots ,{\theta }_{L})}{2{\epsilon}}$ (numerical differentiation with the annealing machine)

5:   ${\theta }_{n}^{(k+1)}{\leftarrow}{\theta }_{n}^{(k)}-\alpha {\left.\frac{\partial {E}^{(\mathrm{a}\mathrm{n}\mathrm{n})}({\theta }_{1},\dots ,{\theta }_{L})}{\partial {\theta }_{j}}\right\vert }_{{\theta }_{n}={\theta }_{n}^{(k)}}$ (update twist parameter)

6:  end for

7: end for

It is worth noting that by performing the scheme shown in algorithm 1 for several values of the annealing time, we can find the optimized annealing time that minimizes E^(ann). This optimized annealing time is denoted by T^(opt). We choose suitable learning rate so that the results should sufficiently converge with a reasonable number of iterations.

This section describes numerical simulations conducted to evaluate the performance of our scheme. In particular, to account for decoherence, we employ the Lindblad master equation. For the problem Hamiltonians, we consider two examples: a XYZ with transverse field model and a deformed spin star model. Moreover, in the appendix, we show numerical results when we choose a hydrogen molecule as the problem Hamiltonian.

5.1. Lindblad master equation

In this subsection, we introduce the Lindblad master equation to consider the decoherence during QA. The Lindblad master equation that we use in this paper is given by

$$\begin{equation}\frac{\mathrm{d}\rho (t)}{\mathrm{d}t}=-\mathrm{i}[H(t),\rho (t)]+\sum\limits _{n}\gamma [{\sigma }_{n}^{(k)}\rho (t){\sigma }_{n}^{(k)}-\rho (t)],\end{equation} \tag{ 9 }$$

where ${\sigma }_{j}^{(k)}(k=x,y,z)$ denote the Lindblad operators acting at site j, γ denotes the decoherence rate, and ρ(t) is the density matrix of the quantum state at time t. It is worth mentioning that this master equation was used as a phenomenological model to describe the decoherence during the QA with the superconducting flux qubits [67]. To treat the decoherence more properly, we should derive the master equation from the first principle without (or with less) approximations as explained in [12, 68]. However, such an analysis is beyond the scope of our paper, and we adopt the widely used model where the Lindblad operator is described by a simple Pauli operator [67].

We solve the Lindblad master equation using QuTiP [69, 70]. Throughout this paper, we choose ${\sigma }_{j}^{z}$ as the Lindblad operator because this type of noise is considered as the main source of decoherence for the qubits used in QA [67].

5.2. XYZ model with transverse field

In this subsection, we discuss the numerical results obtained using a XYZ model with the transverse field as the problem Hamiltonian. In our numerical simulations, the dynamics of the state strongly depends on the annealing time.

First, we introduce the XYZ model. The Hamiltonian of the XYZ is given by

$$\begin{equation}H=\sum\limits _{j=1}^{L-1}\left({\sigma }_{j}^{(x)}{\sigma }_{j+1}^{(x)}+{J}_{1}{\sigma }_{j}^{(y)}{\sigma }_{j+1}^{(y)}+{J}_{2}{\sigma }_{j}^{(z)}{\sigma }_{j+1}^{(z)}+{\Gamma}{\sigma }_{j}^{(x)}\right),\end{equation} \tag{ 10 }$$

where the coefficients J₁ and J₂ are anisotropic parameters and the boundary condition is open. It is known that in the case of Γ = 0 this model is integrable system but in the case of Γ ≠ 0 this model is non-integrable system [71–74]. We choose the parameters as (J₁, J₂, Γ) = (0.4, 0.7, 1.0), and the site number as L = 3. Throughout of our paper, the unit of the parameter in the Hamiltonian is GHz.

We conduct numerical simulations to quantify the performance of our scheme, where we use inhomogeneous twist parameters for the driving Hamiltonian to lower the energy in a variational manner. For comparison, we also perform numerical simulations with conventional QA, where the driving Hamiltonian is chosen as the transverse field. For both cases, we optimize the annealing time. However, the way to optimize the annealing time in our scheme is different from that in the conventional scheme. More specifically, in our scheme, we firstly find an optimal twist angle by using the gradient descent for a given annealing time, and we secondly repeat the first step for many choice of the annealing time, and we finally find the optimal annealing time from the list obtained in the second step. On the other hand, in the conventional scheme, by fixing the twist angle as 0, the annealing time is optimized. We set γ = 10⁻⁴ and α = 0.1. Furthermore, the number of steps is 200. The coherence time of the superconducting qubit used by the device of D-wave is less than 1/γ = 100 ns [75], while the coupling strength is as large as J = 2π × 5 GHz [76]. So the ratio between the coupling strength and decay rate is less than J/γ = 3 × 10³, which is smaller than that adopted with our numerical simulations. On the other hand, the coherence time of the C-shunt flux qubits can be a hundred micro second and the coupling strength will be around J = 2π × 15 MHz [77]. We can achieve J/γ ≃ 10³ with the C-shunt flux qubits, and so our scheme can be realized with C-shunt flux qubits where spin-lock technique is adopted [53].

We plot the energy spectrum of the annealing Hamiltonian for each scheme in figure 1. We will discuss the interpretation of these energy diagrams later. In this paper, we define an estimation error as the difference between the true ground state energy and the measured energy after QA. If our scheme provides a smaller estimation error than the conventional scheme, it is considered to be more accurate than the conventional scheme. In figure 2, we plot the estimation error against the number of variational steps and annealing time using either our variational scheme with the optimal variational parameters or the conventional approach, and we find that the estimation error of our scheme is two order of magnitude smaller than that of the conventional scheme. In addition, we find that the optimal annealing time for minimizing the energy in our scheme is shorter than that in the conventional approach.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Let |⟨φ_drive|φ_prob⟩| denote the overlap, where |φ_drive⟩ denotes the ground state of the driving Hamiltonian and |φ_prob⟩ denotes the ground state of the problem Hamiltonian. In figure 3(a), we plot the overlap between the initial state of twisted QA and the ground state of the problem Hamiltonian. This shows that the ground state of the problem Hamiltonian has a finite distance from that of the driving Hamiltonian in our scheme.

Zoom In Zoom Out Reset image size

To investigate the effect of the decoherence, we use the so-called purity, which is known as a measure for quantifying the effect of decoherence. It is defined by

$$\begin{equation}P=\mathrm{Tr}({\rho }^{2}),\end{equation} \tag{ 11 }$$

where ρ is a density matrix. For a pure state, the purity becomes 1, whereas it becomes exponentially small against the number of qubits for a completely mixed state. We plot the purity to quantify the effect of decoherence (see figure 3(b)). In conventional QA, as we increase the annealing time, the purity decreases owing to decoherence. By contrast, in our scheme, the decoherence effect on the purity is smaller. It is worth mentioning that our method is useful not only for suppressing the decoherence effect but also for making the energy gap larger. We plot the purity after we optimize the twist angle, which makes the energy gap larger. So, actually, we observe the combinational effect of both the suppression of the decoherence and non-adiabatic transitions in the figure 3(b).

Furthermore, we analyzed a so-called adiabatic condition. If the following quantity A_j is much smaller than 1 for all j, the adiabatic condition is satisfied [11, 28–30]:

$$\begin{equation}{A}_{j}=\frac{\vert \langle {E}_{j}\vert \frac{\partial H}{\partial t}\vert {E}_{0}\rangle \vert }{{({E}_{j}-{E}_{0})}^{2}},\end{equation} \tag{ 12 }$$

where the numerator denotes the transition matrix elements of the derivative of the Hamiltonian from the ground state to the jth excited state and the denominator denotes the energy gap between the ground state and the jth excited state. In addition, in figure 1, we plot an energy diagram for the conventional QA, and also plot an energy of a state during the QA with a sufficiently long T without decoherence. We can see that a crossing occurs twice via QA, when we perform a ground state search in the conventional scheme. Actually, the system during the QA is in the first excited state for a time around t/T = 0.3. This suggests that there is a symmetry [78, 79] in the Hamiltonian to induce the level-crossing [66, 80, 81].

Thus, to check whether the adiabatic condition is satisfied, we need to investigate the energy gap and the transition matrices between the first excited state and the jth excited states (j ⩾ 2) when we consider a time around t/T = 0.3 in the conventional scheme. Actually, |E₂ − E₁| for a time around t/T = 0.3 in the conventional scheme becomes much smaller than |E_j − E₀| for j = 1, 2, 3, 4 in our scheme, as shown in figures 4 and 5. Also, we plot the transition matrix that corresponds to the numerator of A_j for conventional QA and our method in figures 6 and 7, respectively. In figures 6(a) and (b), it seems that $\vert \langle {E}_{j}\vert \frac{\partial H}{\partial t}\vert {E}_{0}\rangle \vert$ for j = 1, 2, 3, 4 in conventional scheme is smaller than those of our scheme. On the other hand, we plot the transition matrix between the first excited state and jth excited states (j ⩾ 2) in the conventional scheme, as shown in figure 7. This clearly shows that the transition matrix $\vert \langle {E}_{2}\vert \frac{\partial H}{\partial t}\vert {E}_{1}\rangle \vert$ in the conventional scheme for a time around t/T = 0.3 is larger than $\vert \langle {E}_{j}\vert \frac{\partial H}{\partial t}\vert {E}_{0}\rangle \vert$ for j = 1, 2, 3, 4 in our scheme.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

We investigate how the estimation error and overlap scales with the number of qubits for the XYZ model in figure 8. The results show that, as we increase the number of qubits, the overlap with our scheme becomes smaller and converges to that with the conventional scheme. Also, we observe that, as we increase the number of qubits, the difference between the estimation error with our scheme and that with the conventional scheme becomes larger, which manifests the superiority of our scheme for a larger system. Importantly, these results demonstrates that, even if the overlap is small, our scheme can give a significant improvement. Therefore, our method should work even for larger systems.

Zoom In Zoom Out Reset image size

Therefore, in the case of the XYZ model with transverse field, our method is advantageous in all aspects, namely the effects of decoherence, the energy gap, and the transition matrix of the derivative of the annealing Hamiltonian.

5.3. Deformed spin star model

In this subsection, we consider a deformed spin star model as the problem Hamiltonian where L satellite qubits are coupled with a central qubit. The deformed spin star model Hamiltonian is given by

$$\begin{equation}H=\omega {\hat{\sigma }}_{0}^{z}+{\omega }_{1}{\hat{J}}^{z}+J({\hat{\sigma }}_{0}^{+}{\hat{J}}^{-}+{\hat{\sigma }}_{0}^{-}{\hat{J}}^{+}),\end{equation} \tag{ 13 }$$

where ${\hat{J}}^{+}\equiv {\sum }_{j=1}^{L}{\mathrm{e}}^{2\pi \frac{j}{L}}{\sigma }_{j}^{+}$ and ${\hat{J}}^{-}\equiv {\sum }_{j=1}^{L}{\mathrm{e}}^{-2\pi \frac{j}{L}}{\sigma }_{j}^{-}$. This model has been studied to represent a hybrid system composed of a superconducting flux qubit and nitrogen-vacancy centers in a diamond lattice [82–87]. It is known that the ground state of the deformed spin star model with ω = ω₁ is

$$\begin{equation}\vert {W}_{\theta }\rangle =\frac{1}{L}\sum\limits _{j=1}^{L}{\sigma }_{j}^{+}\,{\mathrm{e}}^{ij\theta }\vert {\downarrow}{\downarrow}{\downarrow}\dots {\downarrow}\rangle .\end{equation} \tag{ 14 }$$

As this state is highly entangled, an overlap with a product state created by the driving Hamiltonian in our scheme cannot be large.

We set γ = 10⁻⁴, α = 0.001, J = 15, ω = ω₁ = 1, and L = 3. Further, the number of steps is 500. We plot the energy spectrum of the annealing Hamiltonian at each time in figure 9. From this graph, the energy spectrum does not seem to have changed significantly; however, we observe that the energy gap between the ground state and the excited states in our scheme is larger than that in the conventional scheme, as shown in figure 12(a).

Zoom In Zoom Out Reset image size

In figure 10, we plot the estimation error against the number of variational steps and the annealing time when we adopt either twisted QA or conventional QA. From these plots, we find that our scheme improves the accuracy by an order of magnitude compared to the conventional scheme. For the annealing time with 80 < T < 400, the performance of our (conventional) scheme is getting worse (better) as we increase the annealing time. This suggests that the decoherence is more relevant than the non-adiabatic transition for our scheme while the non-adiabatic transition is more relevant than the decoherence for the conventional scheme.

Zoom In Zoom Out Reset image size

In figure 11(a), we plot the overlap between the initial state of twisted QA and the ground state of the problem Hamiltonian. This shows that the ground state of the deformed spin star model is highly entangled, and the overlap between the initial ground state of the twisted driving Hamiltonian and the ground state of the problem Hamiltonian cannot be large regardless of the twist parameters. We investigate the effects of decoherence and non-adiabatic transitions. First, we plot the purity to show the effect of decoherence in figure 11(b). We can see that the purity of our scheme is slightly higher than that of the conventional scheme. However, as the difference in purity between them is small, we conclude that our scheme cannot significantly suppress the effect of decoherence in this case.

Zoom In Zoom Out Reset image size

Furthermore, we show the energy gap between the ground state and the first excited state in figure 12(a). The smallest energy gap in the conventional method is twice smaller than that in our method. In the case of the energy gap between the ground state and the second excited state, the results obtained for the conventional scheme and our scheme are nearly the same.

Zoom In Zoom Out Reset image size

In addition, we consider the transition matrix elements of the derivative of the Hamiltonian from the ground state to the jth excited state. The transition matrices between the ground state and the first and second excited states are plotted for conventional QA and our method in figure 13. The difference between the conventional method and our method is rather small.

Zoom In Zoom Out Reset image size

From these results, we conclude that the improvement in the accuracy of our scheme arises from the increase in the energy gap owing to the twisting operations.

Finally, we compare the energy errors of the twisted scheme with that of conventional scheme, when the satellite spin number (i.e. L − 1) is larger. When the system size is large, it is difficult to calculate the Lindblad master equation due to the computational cost. So we solve the time-dependent Schrödinger equation without the decoherence for a fixed annealing time. We choose T = 10 (ns) as the annealing time and choose L − 1 = 2, 3, 4, 5 as the satellite spin numbers. We plot the overlap between the initial state of twisted QA and the ground state of the problem Hamiltonian against the number of the satellite spins in figure 14. These plots show that the overlap becomes smaller as the number of satellites is increased. Also, we plot the estimation error against the number of qubits. On the other hand, we see that the difference in estimation error of energy between the conventional and twisted QA cases increases when the number of satellites number increases. Thus, it can be seen that our method is advantageous even when the number of satellites with small overlap is large. Therefore, this figure shows that our method is effective even when the number of sites increases.

Zoom In Zoom Out Reset image size

In this paper, we proposed a variational method for determining a suitable driving Hamiltonian for QA. We employed a parameterized twist operator to change the driving Hamiltonian, where each spin was rotated with some angle characterized by the twist parameters. Starting from the conventional (transverse-field) driving Hamiltonian, we updated the parameters to minimize the energy of the state after QA until the energy converged to a certain value. For the XYZ model, we show that our method is advantageous for the suppressing decoherence, increasing the energy gap, and decreasing the transition matrix elements. On the other hand, for the deformed spin star model, our method succeeds to make the energy gap larger, and we confirm that the advantage persists even for a larger system.

In summary, our QA approach of using variational methods with twist operators can estimate the ground state energy of the problem Hamiltonian with higher accuracy than conventional QA.

This work was supported by MEXT's Leading Initiative for Excellent Young Researchers and JST PRESTO (Grant No. JPMJPR1919), Japan. This paper is partly based on the results obtained from a project, JPNP16007, commissioned by the New Energy and Industrial Technology Development Organization (NEDO), Japan.

While preparing our manuscript, we became aware of a related work that also uses the change in the driving Hamiltonian in a variational way with QA [88].

In this appendix, we discuss the numerical results obtained using a hydrogen molecule as the problem Hamiltonian. In our numerical simulations, the dynamics of the state strongly depends on the annealing time.

First, we introduce the hydrogen molecule. The Hamiltonian of the hydrogen molecule is described by the second quantized form. We map the Hamiltonian with the second quantized form to a spin Hamiltonian by the Jordan–Wigner transformation. The Hamiltonian of the hydrogen molecule is given by

$$\begin{align}\hfill H& ={h}_{0}I+{h}_{1}{\hat{\sigma }}_{0}^{z}+{h}_{2}{\hat{\sigma }}_{1}^{z}+{h}_{3}{\hat{\sigma }}_{2}^{z}+{h}_{4}{\hat{\sigma }}_{3}^{z}+{h}_{5}{\hat{\sigma }}_{0}^{z}{\hat{\sigma }}_{1}^{z}+{h}_{6}{\hat{\sigma }}_{0}^{z}{\hat{\sigma }}_{2}^{z}+{h}_{7}{\hat{\sigma }}_{1}^{z}{\hat{\sigma }}_{2}^{z}+{h}_{8}{\hat{\sigma }}_{0}^{z}{\hat{\sigma }}_{3}^{z}+{h}_{9}{\hat{\sigma }}_{1}^{z}{\hat{\sigma }}_{3}^{z}\hfill \\ \hfill & \quad +{h}_{10}{\hat{\sigma }}_{2}^{z}{\hat{\sigma }}_{3}^{z}+{h}_{11}{\hat{\sigma }}_{0}^{y}{\hat{\sigma }}_{1}^{y}{\hat{\sigma }}_{2}^{x}{\hat{\sigma }}_{3}^{x}+{h}_{12}{\hat{\sigma }}_{0}^{x}{\hat{\sigma }}_{1}^{y}{\hat{\sigma }}_{2}^{y}{\hat{\sigma }}_{3}^{x}+{h}_{13}{\hat{\sigma }}_{0}^{y}{\hat{\sigma }}_{1}^{x}{\hat{\sigma }}_{2}^{x}{\hat{\sigma }}_{3}^{y}+{h}_{14}{\hat{\sigma }}_{0}^{x}{\hat{\sigma }}_{1}^{x}{\hat{\sigma }}_{2}^{y}{\hat{\sigma }}_{3}^{y},\hfill \end{align} \tag{ A.1 }$$

where we use the STO-3G basis and Jordan–Wigner transformation [89]. The coefficients h₀, h₁, ..., h₁₄ of the Hamiltonian in equation (A.1) depend on the interatomic distance. We obtain these coefficients of the Hamiltonian expressed by the spin in equation (A.1) using OpenFermion for each interatomic distance [90]. We choose the interatomic distance as 0.74 Å, and the coefficients of the Hamiltonian with this interatomic distance are listed in table A1. It is known that the ground state of the Hamiltonian of the hydrogen molecule is very close to the separable state using the Hartree–Fock approximation.

In this case, as will be discussed later, the twist parameters tend to be chosen such that the ground states of the driving Hamiltonian and problem Hamiltonian are very close. We conduct numerical simulations to quantify the performance of our scheme, where we use inhomogeneous twist parameters for the driving Hamiltonian to lower the energy in a variational manner. For comparison, we also perform numerical simulations with conventional QA, where the driving Hamiltonian is chosen as the transverse field. We set γ = 10⁻⁴ and α = 0.05. Furthermore, the number of steps is 200. We plot the energy spectrum of the annealing Hamiltonian for each scheme in figure A1. We observe that the twist operations increase the energy gap between the ground state energy and the first excited state energy. In figure A2, we plot the estimation error against the number of variational steps and annealing time using either our variational scheme with the optimal variational parameters or the conventional approach, and we find that the estimation error of our scheme is one order of magnitude smaller than that of the conventional scheme. In addition, we find that the optimal annealing time for minimizing the energy in our scheme is shorter than that in the conventional approach.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

When we set a short annealing time, our variational method tends to choose a driving Hamiltonian whose ground state has a large overlap with the target ground state of the problem Hamiltonian. By contrast, when we set a long annealing time, our method tends to choose a driving Hamiltonian whose ground state is robust against decoherence.

In figure A3(a), we plot the overlap between the initial state of twisted QA and the ground state of the problem Hamiltonian. From figure A3(a), we see that the overlap becomes especially large for a short annealing time. We plot the purity to quantify the effect of decoherence (see figure A3(b)). In conventional QA, as we increase the annealing time, the purity decreases owing to decoherence. By contrast, in our scheme, the decoherence effect is negligible. This is probably because both the initial state of the driving Hamiltonian in our scheme and the ground state of the problem Hamiltonian are nearly eigenstates of ${\hat{\sigma }}_{z}$, which is robust against ${\hat{\sigma }}_{z}$ noise.

Zoom In Zoom Out Reset image size

We plot the energy gap between the ground state and the first excited state (see figure A4(a)). We can see that our method enlarges the energy gap between the ground state and the first excited state. Similarly, the energy gap between the ground state and the second excited state is increased in our scheme, as shown in figure A4(b). In addition, we plot the transition matrix that corresponds to the numerator of A_j for conventional QA and our method in figures A5(a) and (b), respectively. We find that the transition matrix in our scheme is around two orders of magnitude smaller than that in the conventional one.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Therefore, in the case of the hydrogen molecule, our method is advantageous in all aspects, namely the effects of decoherence, the energy gap, and the transition matrix of the derivative of the annealing Hamiltonian.

In this section, we discuss the number of measurements required for the our method. At first, we decompose the problem Hamiltonian into the local operators as follows.

$$\begin{equation}H=\sum\limits _{j}^{M}{c}_{j}{\hat{a}}_{j},\end{equation} \tag{ B.1 }$$

where M denote the number of local operator, ${\hat{a}}_{j}$ denote local operator, and c_j denote the coefficient of the local operator ${\hat{a}}_{j}$. We can estimate the energy via the measurements of the local operators. The standard deviation of the energy value obtained by measurement is

$$\begin{equation}\delta {H}_{\mathrm{E}}\propto \frac{1}{\sqrt{N}},\end{equation} \tag{ B.2 }$$

where N denote the number of measurements to know the energy for a given state. Thus, if we measure the Hamiltonian, we obtain

$$\begin{equation}\langle H\rangle +\delta {H}_{\mathrm{E}}.\end{equation} \tag{ B.3 }$$

In addition, we obtain the numerical differentiation of E^(ann)(θ₁, ..., θ_L ) with respect to θ_j with measurement error as follows.

$$\begin{align}\frac{\mathrm{\partial }{E}^{(\mathrm{a}\mathrm{n}\mathrm{n})}({\theta }_{1},\dots ,{\theta }_{L})}{\mathrm{\partial }{\theta }_{j}}& =\frac{{E}^{(\mathrm{a}\mathrm{n}\mathrm{n})}({\theta }_{1},\dots {\theta }_{j}-{\epsilon},\dots ,{\theta }_{L})-{E}^{(\mathrm{a}\mathrm{n}\mathrm{n})}({\theta }_{1},\dots {\theta }_{j}+{\epsilon},\dots ,{\theta }_{L})}{2{\epsilon}}\\ & \quad +\delta \left(\frac{\mathrm{\partial }{E}^{(\mathrm{a}\mathrm{n}\mathrm{n})}({\theta }_{1},\dots ,{\theta }_{L})}{\mathrm{\partial }{\theta }_{j}}\right),\end{align}$$

where, $\delta (\frac{\partial {E}^{(\mathrm{a}\mathrm{n}\mathrm{n})}({\theta }_{1},\dots ,{\theta }_{L})}{\partial {\theta }_{j}})$ denotes the standard deviation of the numerical differentiation of E^(ann)(θ₁, ..., θ_L ). Here, the standard deviation of this value is

$$\begin{equation*}\delta \left(\frac{\partial {E}^{(\mathrm{a}\mathrm{n}\mathrm{n})}({\theta }_{1},\dots ,{\theta }_{L})}{\partial {\theta }_{j}}\right)\propto \frac{1}{{\epsilon}}\sqrt{\frac{1}{N}}.\end{equation*}$$

Since the numerical derivatives are performed L times per step, and this is done in K steps, the total standard deviation is

$$\begin{equation}\delta \left(\sum\limits _{k=1}^{K}\sum\limits _{j=1}^{L}\frac{\partial {E}^{(\mathrm{a}\mathrm{n}\mathrm{n})}({\theta }_{1},\dots ,{\theta }_{L})}{\partial {\theta }_{j}}{\vert }_{{\theta }_{n}={\theta }_{n}^{(k)}}\right)\propto \frac{1}{{\epsilon}}\sqrt{\frac{KL}{N}},\end{equation} \tag{ B.4 }$$

where these are assumed to be i.i.d distributed and we use the central limit theorem. Thus, we need to choose the number of measurements N to satisfy

$$\begin{equation}\frac{1}{{\epsilon}}\sqrt{\frac{KL}{N}}\ll 1.\end{equation} \tag{ B.5 }$$

Therefore, the number of measurements required for our protocol is

$$\begin{equation}N\propto \frac{KL}{{{\epsilon}}^{2}}.\end{equation} \tag{ B.6 }$$