Dynamics of Order Parameters of Non-stoquastic Hamiltonians in the Adaptive Quantum Monte Carlo Method

We derive macroscopically deterministic flow equations with regard to the order parameters of the ferromagnetic $p$-spin model with infinite-range interactions. The $p$-spin model has a first-order phase transition for $p>2$. In the case of $p\geq5$ ,the $p$-spin model with anti-ferromagnetic XX interaction has a second-order phase transition in a certain region. In this case, however, the model becomes a non-stoqustic Hamiltonian, resulting in a negative sign problem. To simulate the $p$-spin model with anti-ferromagnetic XX interaction, we utilize the adaptive quantum Monte Carlo method. By using this method, we can regard the effect of the anti-ferromagnetic XX interaction as fluctuations of the transverse magnetic field. A previous study derived deterministic flow equations of the order parameters in the quantum Monte Carlo method. In this study, we derive macroscopically deterministic flow equations for the magnetization and transverse magnetization from the master equation in the adaptive quantum Monte Carlo method. Under the Suzuki-Trotter decomposition, we consider the Glauber-type stochastic process. We solve these differential equations by using the Runge-Kutta method and verify that these results are consistent with the saddle-point solution of mean-field theory. Finally, we analyze the stability of the equilibrium solutions obtained by the differential equations.


I. INTRODUCTION
Finding the best solution in combinatorial optimization problems is computationally intractable. Nevertheless, efficient solutions to such problems have been studied in various fields. Quantum annealing (QA) stochastically solves combinatorial optimization problems with the aid of quantum fluctuations [1,2]. To do so, the cost function is regarded as the physical energy of the system, and the minimizer corresponds to the ground state of the physical system.
In conventional QA, the Hamiltonian is described by H = H 0 + H 1 . The symbol H 0 denotes the target Hamiltonian where we want to solve the ground state. This Hamiltonian consists of z components of Pauli matrices (σ z 1 , · · · , σ z N ), where N is the total number of spins. We add the quantum driver Hamiltonian H 1 as the quantum fluctuation, which is written Here, Γ denotes the strength of the transverse magnetic field, and σ x i is the x component of the Pauli matrix at site i. In QA, we initially set the strength of the transverse magnetic field very high such that the system explores a wide range of the state space to obtain the ground state. We then gradually decrease the strength of the transverse magnetic field. Following the Schrödinger equation, the initial ground state adiabatically evolves into a nontrivial final ground state, which is the solution to the combinatorial optimization problem.
The theoretically sufficient condition to obtain the ground state in QA is assured by the quantum adiabatic theorem [3]. The total evolutionary time τ of the Schrödinger equation to obtain the ground state depends on the minimum energy gap ∆ min from the ground state: τ ≫ ∆ −2 min . In a system with a second-order phase transition, the minimum energy gap polynomially decays with the system size as ∆ min ∝ N −α . Thus, * shunta.arai.s6@dc.tohoku.ac.jp the total evolutionary time polynomially increases: τ ∝ N 2α . In this case, QA efficiently solves the problems. On the other hand, when the system has a first-order phase transition, the minimum energy gap exponentially decays with the system size as ∆ min ∝ exp(−αN). Since the total evolutionary time exponentially increases such that τ ∝ exp(2αN), QA cannot be performed efficiently.
Seki and Nishimori proposed that QA with antiferromagnetic XX interaction can avoid the first-order phase transition for the ferromagnetic p-spin model with infiniterange interactions [4,5]. This entails an exponential efficiency speedup of conventional QA, because the second-order phase transition has a minimum energy gap that polynomially decreases as a function of the system size. A non-stoquastic Hamiltonian like the model with anti-ferromagnetic XX interaction cannot be simulated with the standard quantum Monte Carlo method (QMC) because the non-stoquastic Hamiltonian has positive values in off-diagonal elements in the computational basis to diagonalize the z component of the Pauli matrix and has a negative sign problem [6][7][8][9].
However, a method has been proposed to avoid the negative sign problem involved in a particular class of non-stoquastic Hamiltonians [10]. This method is called the adaptive quantum Monte Carlo method (AQMC). The AQMC treats the effect of the anti-ferromagnetic XX interaction as the fluctuation of the transverse magnetic field by using the delta function and its Fourier integral transformation. We can calculate various physical quantities of the non-stoquastic Hamiltonian by estimating the transverse magnetization and changing the corresponding transverse magnetic field obtained by a saddle-point solution.
In this paper, we focus on the dynamics of the AQMC. To simulate QA, we often utilize the QMC, which is mainly designed for sampling from a Boltzmann distribution. Although the dynamics of QMC differ from those of QA [11], it has been found that some aspects of the dynamics of QA can be expressed by QMC [12,13]. Therefore it is useful to consider the dynamics of QMC or AQMC for QA with and without a non-stoquastic Hamiltonian.
We analyze the dynamics of the order parameters with anti-ferromagnetic XX interaction. For cases with a stoquastic Hamiltonian , Inoue analytically derived macroscopically deterministic flow equations of the order parameter, for example longitudinal magnetization in infinite-range quantum spin systems [14][15][16]. The differential equations with respect to the macroscopic order parameter are obtained from the master equation by considering the transition probability of the Glauber-type dynamics of microscopic states under the Suzuki-Trotter decomposition.
Following this approach, we derive the macroscopically deterministic flow equations of order parameters with antiferromagnetic XX interaction in AQMC. The adaptive transverse magnetic field is changed by the transverse magnetization in AQMC. Therefore we newly introduce the dynamics of transverse magnetization and derive differential equations for magnetization and transverse magnetization. We compare the non-trivial behavior of the dynamics of order parameters with and without anti-ferromagnetic XX interaction.
It is useful to establish a way of simulating a class of nonstoquastic Hamiltonians with a classical computer in order to validate a new quantum annealer. To date, conventional QA with a transverse magnetic field has been implemented in the D-Wave machine [17][18][19]. However, a quantum annealer for non-stoquastic Hamiltonians is being developed. Analyzing the dynamics of order parameters with non-stoquastic Hamiltonians will help us to verify the performance of this new quantum annealer for non-stoquastic Hamiltonians.
The remainder of this paper is organized as follows. In Sec. II, we show the algorithm for AQMC. In Sec. III, we derive the macroscopically deterministic flow equations with respect to a non-stoquastic Hamiltonian from the master equation. In Sec. IV, we analyze the stability of the solutions obtained by the macroscopically deterministic flow equations. In Sec. V, we show the numerical results of the differential equations of order parameters. Finally, in Sec. VI, we summarize our results and discuss future research directions.

II. ADAPTIVE QUANTUM MONTE CARLO METHOD
In this paper, we treat the ferromagnetic p-spin model with infinite-range interactions as the target Hamiltonian which is written as (1) We add the quantum driver Hamiltonian as where γ is the strength of the XX interaction. The partition function of the total Hamiltonian is given by where β is the inverse temperature. We employ the Suzuki-Trotter decomposition to divide the total Hamiltonian into two parts [6]. After that, we introduce the delta function as We utilize the Fourier integral representation of the delta function and introduce the auxiliary variablem xk on the Trotter slice k. We obtain the partition function as where M is the Trotter number. We have dropped a trivial coefficient 1/2π in the above expression. This partition function (5) is equivalent to the partition function of the transverse field Ising model. Furthermore, we use static approximatioñ m xk =m x and m xk = m x to simplify the problem. Finally, the partition function is written as where we define Here, we regard σ z ik as the classical spin σ ik ∈ {−1, +1}. In the thermodynamic limit, we may take the saddle-point in the integral. The saddle-point is evaluated bym x = Γ−γm x . Thus, the instantaneous transverse magnetic field is determined by the transverse magnetization m x . To estimate the transverse magnetization, we consider the conditional probability distribution as where Z(m x ) is the partition function of the effective spin model defined by the conditional probability distribution. The transverse magnetization is written as where the bracket denotes the expectation with respect to the weight of the conditional probability distribution P(σ|m x ). We can realize the non-stoquastic Hamiltonian by estimating the transverse magnetization in the standard QMC method and updating the instantaneous transverse magnetic fieldm x according to the saddle-point solution. The steps for doing so are given in Algorithm 1.

Algorithm 1 adaptive quantum Monte Carlo method
Set the adaptive transverse magnetic fieldm x while the physical quantities converge do 1: Perform the standard QMC for (10 )  In this section, we introduce the dynamics of the p-spin model with XX interaction in AQMC. We consider a local field at site i on the k-th Trotter slice as Here, we neglect the term less than O(1). The effective Hamiltonian is a classical system under the Suzuki-Trotter decomposition. Therefore, the dynamics of AQMC can be written as a Glauber-type stochastic process whose transition probability is given by The master equation for the probability p t ({σ k }), which represents the probability of a macroscopic state including the M-Trotter slices {σ k } ≡ (σ 1 , · · · , σ M ), σ k = (σ 1k , · · · , σ Nk ) at time t is written as follows: where we define the probability of a macroscopic state on the k-th Trotter slice as p t (σ k ) and a spin flip operator We impose the periodic boundary conditions σ 1 = σ M+1 . Next, we introduce a probability distribution P t ({m k } , {m xk }) of the macroscopic order parameters as where we define K = −(1 − tanh 2 (βm x /M))/(2 tanh(βm x /M)), the set of the longitudinal magnetization as {m k } ≡ (m 1 , · · · , m M ), and the set of the transverse magnetization as {m xk } ≡ (m x1 , · · · , m xM ). The notation σ k is written as We regard m k as the magnetization on the Trotter slice k, and m xk as the transverse magnetization between the Trotter slice k and k + 1. Taking the derivative of Eq. (18) with respect to t, we can obtain differential equations as From Eq. (20), we can obtain these deterministic flow equations as Here, we use the saddle-point solutionm x = Γ − γm x . The derivations of Eqs. (21) and (22) are described in Appendix A.

IV. STABILITY ANALYSIS OF THE EQUILIBRIUM SOLUTIONS
In this section, we derive the stability of the solutions obtained by the deterministic flow equations [20][21][22]. To simplify the problem, we consider the zero-temperature limit β → ∞. We can rewrite the deterministic flow equations (21) and (22) as We assume the existence of the equilibrium solutions (m * , m * x ). These solutions satisfy f (m * , m * x ) = 0 and g(m * , m * x ) = 0. We consider infinitesimal increments of m and m x around the equilibrium solutions as The Taylor expansions for f (m, m x ) and g(m, m x ) around the equilibrium solutions yield and From Eqs. (25)-(28), the time differential of infinitesimal increments are written as follows: where we define the Jacobian matrix as To determine the stability of the equilibrium solutions, we consider the characteristic polynomial of the Jacobian matrix whose eigenvalues are λ 1 and λ 2 . By evaluating the trace of the Jacobian matrix trJ = λ 1 + λ 2 , the determinant detJ = λ 1 · λ 2 , and whether the eigenvalues are complex, we can investigate the stability of the equilibrium solutions. The condition with real roots of eigenvalues is (trJ) 2 > 4 · detJ. This condition determines whether the equilibrium solutions have oscillations.

V. EXPERIMENTAL RESULTS
In this section, we describe our numerical experiments. We numerically solve differential equations (21) and (22) using the Runge-Kutta method with a sufficiently low-temperature T and inverse temperature β = 1/T = 100. The ferromagnetic p-spin model has a first-order phase transition for p > 2 [23]. We set the parameters γ = 0 or γ = 18. In the case of γ = 18, the p-spin model has a second-order phase transition for p > 4 following the phase diagram in [4]. In this study, we consider p = 5 and compare the dynamics of order parameters with and without anti-ferromagnetic XX interaction.
First, we show the dynamics of order parameters without anti-ferromagnetic XX interaction γ = 0 in Figs. 1 and 2. The original model has a first-order phase transition. We set three different conditions: Γ = 0.5 in the ferromagnetic phase, Γ = 2.5 in the paramagnetic phase, and Γ = 1.3 in the coexisting phase between the critical point Γ c and the spinodal point Γ sp . These figures indicate that the dynamics exponentially converge to each steady state depending on the initial condition.
We also consider the pseudo free energy to evaluate the equilibrium solutions. By using mean-field theory, the pseudo free energy is written as Standardly, with mean-field theory, we utilize the saddle-point conditions ∂F/∂m = 0 and ∂F/∂m x = 0 form andm x , respectively. Strictly speaking, however, these conditions need not be used to estimate the pseudo free energy precisely (31). Therefore, we utilize the saddle-point conditions ∂F/∂m = 0 and ∂F/∂m x = 0. These conditions lead to Even if we utilize these saddle-point conditions, we can ultimately obtain the saddle-point equations as These saddle-point equations, obtained in the standard manner for mean-field theory [4], are consistent with dm/dt = 0 and dm x /dt = 0.
To show the validity of the solution obtained from the master equation, we plot the pseudo free energy (31) with respect to the function of longitudinal magnetization m in Fig.3 Here, we can utilize the equation m 2 + m 2 x = 1 in the thermodynamic limit N → ∞. According to the initial values, each equilibrium solution from the master equation converges to the minimum values. From Fig.3(a), the pseudo free energy has two different stable values m ≃ 1 and m ≃ 0. The solution m ≃ 0 is the metastable state in the ferromagnetic phase. Therefore, we can see that the equilibrium solutions in Fig.1(a) converge to two different values, m ≃ 1 and m ≃ 0, according to the different initial values. In the coexisting phase, a similar phenomenon obtains.
We plot the equilibrium solutions from the master equation and the exact solutions from the saddle-point equations in Fig. 4. We find that these order parameters change discontinuously. The longitudinal magnetization is the multivalued function with respect to the strength of the transverse magnetic field. After the spinodal point, the ferromagnetic stable state m > 0 appears. The dotted line in Fig. 4 denotes the critical point where the pseudo free energy takes the same value. From the viewpoint of dynamics, we can confirm that this model has a first-order phase transition.
The dynamics of these order parameters with antiferromagnetic XX interaction γ = 18 are shown Figs. 5 and 6. We set three different conditions: Γ = 5 in the ferromagnetic phase, Γ = 15 in the critical phase, and Γ = 25 in the paramagnetic phase. These figures indicate that the dynamics exponentially converge to the only steady state.
We plot the pseudo free energy in Fig. 7. Finally, we investigate the stability of the equilibrium solutions obtained by the deterministic flow equations. We numerically compute the trace of the Jacobian matrix, the determinant, and the condition with real roots of eigenvalues by utilizing the equilibrium solutions (m * , m * x ) of Eqs. (23) and (24). We plot the behavior of the trace of the Jacobian matrix, the determinant, and the condition with real roots of eigenvalues with and without anti-ferromagnetic XX interaction in Fig. 9. Figure 9(a) shows that the equilibrium solutions are stable nodes when the equilibrium solutions are minimum values of the free energy because (trJ) 2 > 4 · detJ, detJ > 0, and trJ < 0 hold. We can see the differences between Fig. 9(a) and Fig. 9(b). When the equilibrium solutions are not necessarily minimum values of the free energy, the equilibrium solutions are saddled and unstable near the spinodal point. Then, (trJ) 2 > 4 · detJ, detJ < 0, and trJ < 0 are established.
In the case of a non-stoquasitc Hamiltonian γ = 18 shown in Fig. 9(c), the equilibrium solutions are stable nodes in the ferromagnetic phase and the paramagnetic phase. Near the point Γ = γm x , the equilibrium solutions are saddled and unstable. Subsequently, the equilibrium solutions reach stability. As the strength of the transverse magnetic field decreases, the equilibrium solutions become a stable spiral and are asymptotically stable because (trJ) 2 < 4 · detJ, detJ > 0 and trJ < 0 hold . In this region, there is an oscillation because the eigenvalues of the Jacobian matrix have complex values. Thus, strong quantum fluctuations like the anti-ferromagnetic XX interaction affect the stability of the equilibrium solutions.

VI. CONCLUSION
We derived macroscopically deterministic flow equations of order parameters from the Glauber-type master equation under the Suzuki-Trotter decomposition. In AQMC, the adaptive transverse magnetic field is governed by transverse magnetization. By changing the adaptive transverse magnetic field in accordance with the saddle-point solution, we can obtain the dynamics of order parameters of the p-spin model with and without anti-ferromagnetic XX interaction. We found that the equilibrium solutions obtained by the deterministic flow equations are identical to the saddle-point solutions obtained by mean-field theory. We can obtain the behavior of order parameters until the system is equilibrated. Without anti-ferromagnetic XX interaction, the metastable state appeared near the spinodal point because the original model has a first-order phase transition. In the ferromagnetic phase and the coexisting phase, the equilibrium solutions converged to different values depending on the initial conditions, owing to the existence of the metastable solution. By adding the anti-ferromagnetic XX interaction, the metastable solution vanished. Therefore, the equilibrium solutions converged to the saddle-point solution in all phases. We confirmed that the equilibrium solutions minimize the free energy. Finally, we investigated the stability of the equilibrium solutions under a zero-temperature limit. When the equilibrium solutions include metastable solutions in the case without antiferromagnetic XX interaction, the stability of the equilibrium solutions are unstable. For the model with anti-ferromagnetic XX interaction, the stability of the equilibrium solutions significantly changed compared to the original model. We found that these strong quantum fluctuations have an impact on the stability of the equilibrium solutions.
This approach of using the master equation is useful for the dynamics of QA, which is inhomogeneous the Markovian stochastic process where the strength of the transverse magnetic field is time-dependent. This helps us to investigate not only conventional QA but also QA with a non-stoquastic Hamiltonian, inhomogeneous driving of the transverse field, and reverse annealing [24][25][26]. In our future work, we will analyze the dynamics of these new types of QA. We here drive the differential equations (21) and (22). To derive the deterministic flow equations of order parameters, we utilize the following assumptions: where we define the effective single local field βΦ(k) ≡ (βp/M)m p−1 k + B/2(σ k−1 + σ k+1 ) on the k-th Trotter slice and the average · · · \σ k ≡ lim M→∞ σ\σ k (· · · )p(σ 1 , . . . , σ k−1 , σ k+1 , . . . , σ M ) Here, we write the summation with respect to all sites except for the Trotter slice k as We can rewrite tanh(βΦ(k)) and σ k+1 tanh(βΦ(k)) as σ k σ k+1 p(σ k |σ 1 , . . . , σ k−1 ). (A6) In a manner similar to a previous study [15], we can obtain the expectation of tanh(βΦ(k)) and that of Kσ k+1 tanh(βΦ(k)) In order to derive a compact representation of the differential equations, we substitute P t ({m k } , {m xk }) = M k=1 δ(m k − m k (t))δ(m xk − m xk (t)) into Eq. (A9) and carry out the integral with respect to k m k and k m xk after multiplying itself m k . Finally, we can obtain the differential equations for each Trotter slices k as dm k dt = −m k + σ k path .
In order to derive a compact representation of the differential equation, we utilize the static approximation m k = m, m xk = m x . Under this approximation, we inverse the procedure of the Suzuki-Trotter decomposition:  23) and (24). The horizontal axis denotes the strength of the transverse magnetic field. In Figs. 9(a) and 9(b) without anti-ferromagnetic XX interaction γ = 0, the equilibrium solutions do and do not include the metastable solutions, respectively. In Fig. 9(c), we consider the case with anti-ferromagnetic XX interaction γ = 18.
We can regard σ k path = lim M→∞ M −1 k σ k path as We substitute Eq. (A12) for Eq. (A10) and obtain the deterministic equation (21). For m xk , we similarly consider the flow equation as