Phase transitions in the frustrated Ising ladder with stoquastic and nonstoquastic catalysts

The role of nonstoquasticity in the field of quantum annealing and adiabatic quantum computing is an actively debated topic. We study a strongly-frustrated quasi-one-dimensional quantum Ising model on a two-leg ladder to elucidate how a first-order phase transition with a topological origin is affected by interactions of the $\pm XX$-type. Such interactions are sometimes known as stoquastic (negative sign) and nonstoquastic (positive sign)"catalysts". Carrying out a symmetry-preserving real-space renormalization group analysis and extensive density-matrix renormalization group computations, we show that the phase diagrams obtained by these two methods are in qualitative agreement with each other and reveal that the first-order quantum phase transition of a topological nature remains stable against the introduction of both $XX$-type catalysts. This is the first study of the effects of nonstoquasticity on a first-order phase transition between topologically distinct phases. Our results indicate that nonstoquastic catalysts are generally insufficient for removing topological obstacles in quantum annealing and adiabatic quantum computing.


I. INTRODUCTION
Quantum annealing exploits quantum-mechanical fluctuations to solve combinatorial optimization problems [1][2][3][4][5][6][7][8]. In a typical formulation, the combinatorial optimization problem one wants to solve is expressed as an Ising model, represented in terms of the z components of the Pauli matrices, and quantum fluctuations are introduced as a transverse field (sometimes called the X term), the sum of the x components of the Pauli matrices over all sites. One of the bottlenecks of quantum annealing under unitary Schrödinger dynamics is a firstorder quantum phase transition as a function of the relative weight of the Ising Hamiltonian and the X term, which exists in almost all cases of interest. At a firstorder quantum phase transition, the energy gap between the ground state and the first excited state closes exponentially as a function of the system size except for a few very limited cases [9][10][11]. This leads to exponentially long computation times within the framework of * Present address: Japan Digital Design, Inc., 3-3-5 Nihonbashihongokucho, Chuo, Tokyo 103-0021, Japan; mtcamkxacdiaiki@gmail.com adiabatic time evolution, according to the adiabatic theorem of quantum mechanics [12,13]. 1 There has been a large amount of effort to circumvent this difficulty including, but not limited to, the introduction of nonstoquastic catalysts [14][15][16][17][18][19][20][21][22][23][24], inhomogeneous field driving [25][26][27][28][29][30][31], reverse annealing [32][33][34][35][36][37][38][39][40][41][42][43], pausing [38,[44][45][46], and diabatic quantum annealing [47]. In the present paper, we focus on the approach of nonstoquastic catalysts. A Hamiltonian is called stoquastic if there is a choice of a local (tensor-product) basis in which the Hamiltonian off-diagonal matrix elements are all real and nonpositive [48]. Otherwise the Hamiltonian is called nonstoquastic, and the inevitable positive or complex offdiagonal matrix elements of the Hamiltonian lead to the infamous sign problem [49,50]. We note that even when the Hamiltonian is stoquastic, but it is presented in a form in which this stoquasticity is unapparent, the problem of deciding whether there exists a local transformation to a basis that "cures the sign problem" (makes it stoquastic) by making all of the Hamiltonian matrix elements real and nonpositive, is NP-complete [51,52].
Many interesting quantum models are stoquastic, e.g., the transverse-field Ising model of conventional quantum annealing, the Bose-Hubbard model, and particles subject to a position-dependent potential, which includes superconducting flux circuits if we associate flux with position and charge with momentum [48,53]. The partitionfunction decomposition of a stoquastic Hamiltonian leads to a sum of nonnegative weights (i.e., the absence of sign problem, which means that the path-integral quantum Monte Carlo algorithm [54] can typically be used efficiently for sampling tasks). Moreover, the ground state of a stoquastic Hamiltonian has only nonnegative amplitudes by the Perron-Frobenius theorem. These two observations are often cited as reasons that stoquastic Hamiltonians should be considered to be less computationally powerful than their nonstoquastic counterparts. In addition, it is well known that nonstoquastic adiabatic quantum computation is universal [55][56][57], while universality in the stoquastic case requires excited states [58]. From a computational complexity perspective, the class StoqMA associated with the task of estimating the ground-state energy of stoquastic local Hamiltonians is known to be no harder than QMA and no easier than MA (and hence no easier than NP) [59] (see Ref. [7] for more background).
However, recently evidence has been mounting that nonstoquasticity does not have a clear-cut computational benefit even in the ground state setting. For example, very recently examples were found for which evolution in the ground state of a stoquastic Hamiltonian can solve a problem superpolynomially [60] or even subexponentially [61] faster than is possible classically. In addition, there is both theoretical and numerical evidence that adiabatic paths based on nonstoquastic Hamiltonians generically have smaller gaps between the ground state and the first excited state, with the implication that they are less useful than stoquastic Hamiltonians for adiabatic quantum optimization [24]. In this work we contribute further evidence that nonstoquasticity is not necessarily useful for efficiently solving a problem.
The setting for our work is the conventional transversefield Ising model used in quantum annealing, but with added antiferromagnetic XX interactions (XX terms with positive coefficients). Such terms, which can turn a stoquastic Hamiltonian into a nonstoquastic one, are sometimes called nonstoquastic "catalysts" when they are turned on only at intermediate times (typically as s(1 − s), where s ∈ [0, 1] is the dimensionless time along the anneal) [7,21], a terminology inspired by its use in entanglement theory [62]. 2 One of the outstanding features of the introduction of 2 The first use of a nonstoquastic catalyst in adiabatic quantum optimization dates back to Ref. [14], which referred to it as an "extra piece of the Hamiltonian." In that work the catalyst did an XX-type nonstoquastic catalyst lies in the universality and QMA-completeness of the resulting Hamiltonian, assuming independent control of each term in the Hamiltonian [63]. Our focus is, however, on another aspect, namely, that a certain set of nonstoquastic catalysts of the XX type is known to reduce the order of quantum phase transitions from first to second [15,16,18,[20][21][22]. This means that adiabatic evolution converges to the ground state of the final Hamiltonian in quantum annealing with an exponential speedup relative to the stoquastic case. However, it does not guarantee a quantum speedup relative to classical solution methods, and in addition examples are known where XX-type nonstoquastic catalysts do not lead to performance improvements [19,22,24,25]. All of these studies concern problems with relatively simple phase transitions without a drastic change in the topological structure of thermodynamic phases. The frustrated ladder model introduced by Laumann et al. [11] is an exceptional case in that it undergoes a first-order quantum phase transition between topologically different phases despite the simplicity of the problem, which is defined on a quasi-one-dimensional two-leg ladder with nearest-neighbor interactions. The model is closely related to a dimer problem defined on a dual ladder, through which the topological aspect of its phase transition can be understood intuitively, as will be detailed in Sec. II. For this reason, we study the effects of nonstoquastic XX catalysts on the phase transition of the frustrated Ising ladder, and the stoquastic case is also treated for completeness.
We employ theoretical and numerical tools, i.e., the real-space renormalization group (RG) method and the density-matrix renormalization group (DMRG) method, to study the structure of the phase diagram with and without XX catalysts of both signs (stoquastic and nonstoquastic). We find that the phase diagrams obtained by the two methods are qualitatively consistent with each other and conclude that the XX-type catalysts, both stoquastic and nonstoquastic, keep the order of the phase transition intact. This is the first example where the role of stoquastic and nonstoquastic catalysts is revealed systematically in a low-dimensional system that exhibits a first-order phase transition with a change of the topological structure.
This paper is organized as follows: In the next section we formulate the problem and describe the known properties of the model system. In Sec. III we analyze the problem via the real-space RG method. The results of extensive numerical computations are explained in Sec. IV and are compared with those from the real-space RG method. We conclude in Sec. V. Additional technical details are given in the Appendices.
not have the specific XX form but was a random Hermitian (hence nonlocal) matrix.

II. MODEL
We consider the frustrated Ising ladder with transverse fields and XX interactions. This is a generalization of the model proposed and analyzed in Ref. [11], where the transverse field was uniformly applied to all sites and no XX interactions were taken into account. This section describes the definition of the model and its basic properties in the case without XX interactions, and is largely a recapitulation of Ref. [11].
A. Definition of the model As depicted in Fig. 1, the system is composed of qubits (spin-1/2 particles) located on sites of a two-leg ladder. The ladder has a top row (t) and a bottom row (b). Nearest-neighbor interactions are of magnitude K > 0 with ferromagnetic (solid lines) or antiferromagnetic (dashed lines) signs. Local longitudinal fields applied to the top and bottom rows have magnitudes K and U/2 > 0, respectively, and are oppositely directed. With transverse fields (X terms) of magnitude Γ t or Γ b and transverse interactions (XX terms) with magnitude Ξ tt , Ξ tb , or Ξ bb , the Hamiltonian is written aŝ whereX a,i ,Ŷ a,i , andẐ a,i are the Pauli operators at sites i = 1, . . . , L on row a = t, b. We assume the length L to be even and impose the periodic boundary conditionŝ X a,L+1 =X a,1 ,Ŷ a,L+1 =Ŷ a,1 ,Ẑ a,L+1 =Ẑ a,1 . (2) The classical part ofĤ (the terms involving Z operators) is highly frustrated due to the competition between positive and negative interactions as well as between interactions and longitudinal fields. Let us discuss the stoquasticity of the HamiltonianĤ. As is easily checked,Ĥ is stoquastic for Ξ tt , Ξ bb , Ξ tb ≥ 0. Additionally, when Ξ tt < 0, Ξ bb < 0, or Ξ tb < 0, there are cases where a local curing transformation makes the Hamiltonian stoquastic. For example, if Γ t = Γ b = 0 and Ξ tt , Ξ bb , Ξ tb ≤ 0, consider the transformation obtained by conjugating some qubits byẐ a,i , an operation under whichX a,i → −X a,i . Then, the following is a curing transformation: for odd i, conjugate the qubits on the top row and for even i, conjugate the bottom row. This transformation is equivalent to flipping the signs of Ξ tt , Ξ bb , and Ξ tb when It can be shown [64] thatĤ remains nonstoquastic under single-qubit Clifford transformations if the following set of conditions is satisfied: Since the general problem of deciding whether a local curing transformation exists is NP-complete even for singlequbit Clifford transformations [51,52], we do not consider here the nonstoquasticity ofĤ in more general cases than the conditions given by Eq. (3). In the following, we refer toĤ as nonstoquastic if there is no curing transformation that is a product of single-qubit Clifford unitaries, and as stoquastic otherwise.
B. Phase diagram for the case without XX terms Laumann et al. [11] studied the phase diagram in the case of the uniform transverse field Γ t = Γ b = Γ and no XX interactions Ξ tt = Ξ bb = Ξ tb = 0 using numerical diagonalization of small-size systems and perturbation from the large-K and small-K limits. We confirm their findings in Fig. 2, which shows that a first-order transition exists as a function of Γ/U with K/U fixed to a large value, while a second-order transition appears for small K/U . The values of Γ/U that minimize the energy gap between the ground state and the first excited state (which we refer to henceforth as "minimum gap are the locations of the minimum energy gap for fixed values of K/U . Here, we set the system size to L = 10. In both of (a) and (b), the first-and second-order phase boundaries predicted by perturbation theory are shown as solid and dotted lines, respectively. The first-order transition line for K/U 1 is Γ/U ≈ 1/c + U/(4Kc 3 ) (c ≈ 0.6) and the second-order transition points") for L = 10 and fixed values of K/U are indicated in Fig. 2 by black squares. According to perturbation theory, the first-order transition line for K/U 1 is Γ/U ≈ 1/c + U/(4Kc 3 ) (c ≈ 0.6) and the second-order transition line for K/U 1 is Γ/U ≈ K/U [11]. These two lines meet at (Γ/U, K/U ) ≈ (2.2, 2.2).
Note that the perturbation theory prediction agrees with the location of the minimum gap points for K/U 2.2 or K/U 1.5, but the agreement breaks down for 1.5 K/U 2.2, where the numerically computed locations of the minimum gap deviate from the perturbative second-order transition line. A jump in the minimum gap points is observed at K/U ≈ 1.5 between Γ/U ≈ 1.5 and 1.9. We provide an explanation of this phenomenon in terms of the appearance of a "double-well" in the energy gap as a function of Γ/U at the critical value K/U ≈ 1.5, which is the origin of the observed discontinuity. See Appendix A for additional details. Figure 2 also shows the staggered magnetization of the top row and the magnetization of the bottom row for L = 10, where · · · denotes the expectation value in the ground state. We see that m t and m b have discontinuities as functions of Γ/U for large K/U , which indicate the existence of the first-order transition. On the other hand, m t and m b change continuously around the secondorder transition, which occurs as Γ/U is decreased with K/U fixed to a small value. The "symmetric" and "staggered" phases shown in Fig. 2 are associated with columnar and staggered configurations, which are defined as follows: • Columnar configurations are product states of local eigenstates ofẐ a,i in which all the bottom-row spins are up and there are no nearest-neighbor (consecutive) down spins on the top row.
• Staggered configurations are product states of local eigenstates ofẐ a,i in which all the bottom-row spins are down and nearest-neighbor top-row spins are antiparallel (antiferromagnetically ordered).
When K is large, the two phases can be characterized by perturbation theory in K −1 . The leading-order part of the state in the symmetric phase is a superposition of columnar configurations (e.g., the red arrows in the top panel of Fig. 3), while that in the staggered phase has staggered configurations (red arrows in the middle and bottom panels of Fig. 3). The name "symmetric phase" means that all columnar spin configurations have comparable amplitudes with the wave function exhibiting translational invariance, and this invariance is broken in the staggered phase. The system is in the symmetric phase for large Γ/U and is in the staggered phase for small Γ/U . Reference [11] also verified that the energy gaps at the first-order transition points decay exponentially as the length L grows and that the decay rate is proportional to ln(K/Γ). This can be intuitively understood from the fact that the transition from a columnar configuration to a staggered configuration requires flipping all the bottom spins as shown in Fig. 3.
As reviewed in Appendix B 1, the ground states of the nonperturbative Hamiltonian proportional to K, are the columnar and staggered configurations. In the presence of the perturbative terms with the coefficients U and Γ, the degeneracy ofĤ (0) is lifted as described above.
The number of columnar configurations is the sum of two Fibonacci numbers F L−1 + F L−3 , which is exponentially large in the length L. Here, F L is defined by the recurrence relation F L = F L−1 +F L−2 and the initial values F 1 = 2 and F 2 = 3. On the other hand, there are only two staggered configurations. These basic properties are confirmed in Appendix B 1.
In the opposite limit U → ∞, the bottom spins are fixed to down, and the top row of the original Ising ladder (1) with Γ t = Γ b = Γ and Ξ tt = Ξ bb = Ξ tb = 0 is effectively an antiferromagnetic Ising chain in a transverse field, L i=1 (KẐ t,iẐt,i+1 − ΓX t,i ). This results in a second-order transition approximately at Γ = K for small K/U , as shown in Fig. 2. The symmetric and staggered phases of the ladder correspond to the quantum paramagnetic and antiferromagnetic phases on the top row, respectively. C. Dimer model on the dual lattice A significant property of the frustrated Ising ladder is that in the limit of large frustration K → ∞ the model is equivalent to the quantum dimer model on a two-leg ladder, as shown in Ref. [11]. We can find this equivalence by locating dimers on the dual lattice according to the types and positions of frustration exhibited by the columnar and staggered configurations. The dual lattice is a two-leg ladder whose position is shifted horizontally and vertically by half of the unit length from the Ising ladder (see Fig. 3). Types of frustration in the columnar and staggered configurations are classified as (i) a down spin on the top row (opposite to the longitudinal field), (ii) a horizontally aligned ferromagnetic pair on the top row (opposite to the antiferromagnetic interaction), and (iii) a vertically aligned antiferromagnetic pair (opposite  Fig. 1 (namely, the black circles, the black solid lines, and the black dashed lines are longitudinal fields, ferromagnetic interactions, and antiferromagnetic interactions, respectively, of equal magnitude K), while the dual lattice for the latter model is indicated by the dotted lines. The rightmost points are identified with the leftmost points for each lattice, which is subject to periodic boundary conditions. The red arrows are spins and the blue thick line segments are dimers. The topological number is w = 0, +1, −1 from top to bottom, the first being a columnar configuration and the second and the third being staggered ones.
to the ferromagnetic interaction). In the dimer picture these three types become (i) a top horizontal dimer, (ii) a vertical dimer, and (iii) a bottom horizontal dimer, respectively. In this manner the columnar and staggered configurations in the frustrated Ising ladder are mapped one-to-one onto hardcore dimer coverings on the dual two-leg ladder.
Dimer coverings on a two-leg ladder are classified into three topological sectors: the columnar sector w = 0 and the two staggered sectors w = ±1, where w denotes the difference in the number of dimers between the top and bottom rows on an arbitrarily given unit square (plaquette). The fact that one cannot transform a dimer covering into another dimer covering with a different w by a series of local movements of dimers allows us to regard w as a topological number.
From this mapping of the columnar and staggered configurations in the frustrated Ising ladder onto dimer coverings, it follows that the Hamiltonian of the frustrated Ising ladder (1) in the limit K → ∞ is equivalent to the following Hamiltonian of the quantum dimer model: where a constant energy shift was ignored. The summations | , , and are performed over vertical lines |, plaquettes , and pairs of two neighboring plaquettes on the dual lattice, respectively. The line segments in the kets and bras denote dimers that are located in the summation "variables" |, , and . We call the limit K → ∞ the dimer limit. Note that the terms with the coefficients Γ b , Ξ bb , and Ξ tb do not appear in the dimer model, because the bottom spins should exhibit complete ferromagnetic order in the limit K → ∞.
Since the only nonvanishing matrix elements ofĤ dimer are those between states in the columnar sector, it naturally follows that there is a strict (not avoided) energylevel crossing between the columnar and staggered sectors in the quantum dimer model on a two-leg ladder. Indeed, the HamiltonianĤ dimer has a strict level crossing at Γ t /U ≈ 1/0.6 in the Ξ tt = 0 case (vanishing XX interactions on the top row) according to numerical diagonalization results for small-size systems [11]. Although for the frustrated Ising ladder with large but finite K the level crossing turns into an avoided crossing [due to nonvanishing transition probabilities to defect states with energy penalties O(K)], the numerical consequence that the energy gap at the first-order transition decays exponentially at least for Ξ tt = 0 reflects the topological nature relating to the quantum dimer model [11].

III. REAL-SPACE RG ANALYSIS
We perform the real-space RG analysis of the frustrated Ising ladder in the limit of large frustration K → ∞, namely the dimer limit, in Sec. III A. Subsequently we analyze the limit of small frustration U → ∞, in Sec. III B.

A. Large frustration limit
We analyze the zero-temperature phase transition of the frustrated Ising ladder in the dimer limit K U, Γ a , Ξ aa using the real-space RG method. Technical details are delegated to Appendix B.
Although the standard real-space RG transformation amounts to separating the Hamiltonian into intrablock and interblock terms and projecting the Hilbert space onto a low-energy space of the intrablock Hamiltonian [65], this procedure will fail to find a low-energy subspace due to the presence of the strong interblock interactions with magnitude K. To circumvent this problem, we employ a real-space RG method in which the projector onto a low-energy space is variationally determined such that the dimer structure is preserved at each RG step. Here, the dimer structure means that the columnar and staggered configurations remain the ground states of the leading-order Hamiltonian in K −1 . Since the columnar and staggered configurations will be the low-energy states with a small energy splitting near the transition, such a variational ansatz may enable us to extract critical properties of the entire system at zero temperature.

Effective Hamiltonian
To write down the RG equations, let us define the generalized Hamiltonian that appears in our RG analysis: where nondimer configurations |σ indicate product states of local eigenstates ofẐ a,i that are neither columnar nor staggered configurations. The coefficients 0 < K σ = O(K) stand for energy penalties on nondimer configurations σ. We denoted Γ t = Γ and Ξ tt = Ξ for notational simplicity. The longitudinal field on the top row V is produced after a step of the RG transformation. The operatorÔ has a magnitude comparable to U , Γ, V , and Ξ, but is irrelevant in the sense of RG theory. The couplings Γ b , Ξ bb , and Ξ tb are included inÔ. For a more detailed version of the generalized Hamiltonian, see Eq. (B14) (although the operators with the coefficients Γ, V , and Ξ are slightly different from Eq. (8), the differences are irrelevant). The bare Hamiltonian (1) can be obtained by setting V = 0 except for a constant energy difference proportional to K, as described in Appendix B 1.

RG equations
Let us perform the RG transformation repeatedly. We denote the coupling constants that have been renormalized l times by U (l), Γ(l), V (l), and Ξ(l). The bare couplings are U (0), Γ(0), V (0) = 0, and Ξ(0). Our RG transformation has the scaling factor b = 3 (note that the scaling factor b should be an odd number to keep the antiferromagnetic order in the staggered phase through the RG transformation) and the system size scales as is the length of the original ladder. As derived in Appendix B, the renormalized coupling constants are determined by the RG equations Here, α 1 (l), α 2 (l), z(l), β 1 (l), β 2 (l) ∈ R are the variational parameters that minimize the function under the constraint where ϕ = (1 + √ 5)/2 is the golden ratio. This minimization is equivalent to minimizing the trace of the renormalized Hamiltonian in the subspace spanned by the columnar and staggered configurations at each RG step.
Note that the function (10) and the constraint (11) are invariant under the transformations (α 1 , α 2 , z) → (−α 1 , −α 2 , −z) and (z, β 1 , β 2 ) → (−z, −β 1 , −β 2 ). Each of these transformations is equivalent to changing the phase of a variational state included in the projector onto the coarse-grained space and thus leaves the projector invariant. Therefore, the optimal set of the variational parameters α 1 (l), α 2 (l), z(l), β 1 (l), β 2 (l) ∈ R is four-fold degenerate and we choose one of the solutions. The choice among the four solutions does not affect the critical properties of the system since the multiplication of (α 1 (l), α 2 (l), z(l)) or (z(l), β 1 (l), β 2 (l)) by −1 only flips the sign of the renormalized transverse field Γ(l + 1), which corresponds to the gauge transformation Let us describe important features of the RG equations. Since Eqs. (9b)-(9d) do not contain U (l) explicitly and the variational parameters minimizing f Γ(l),V (l),Ξ(l) subject to the constraint (11) can be regarded as functions of Γ(l), V (l), and Ξ(l), the renormalized coupling constants Γ(l + 1), V (l + 1), and Ξ(l + 1) are expressed as functions of Γ(l), V (l), and Ξ(l), not including U (l). In addition, it is convenient to write the coupling constant U (l) as which yields the initial value U (0) = 0 and the recurrence relation Now Eqs. (9b)-(9d) and (13) are independent of U (l) as well as U (0). Since U (0) = V (0) = 0, it follows by mathematical induction that U (l), Γ(l), V (l), and Ξ(l) are determined only by the "time" l and the initial values Γ(0) and Ξ(0). Although the penalty constants K σ for nondimer configurations σ are also renormalized, their specific values are unimportant for our analysis of the zero-temperature phase transition in the dimer limit. The only essential point is that all the penalty constants are positive and proportional to K. This means that any spin configuration that is not mapped onto a dimer covering does not contribute to the leading-order part of the ground state.

RG flow and fixed points
We show the coupling constants U (l), Γ(l), V (l), and Ξ(l) as functions of l for the initial values (Γ(0), Ξ(0)) = (1, −1), (1, 0), (1, 1) in Fig. 4 [see the next paragraph on U (l) and its initial value U (0)]. For the behavior of the coupling constants that have other initial values and the behavior of the variational parameters, see Figs. 14-16 in Appendix B 7. Note that multiplying all the bare couplings by a positive constant c leaves the variational parameters at each l unchanged and multiplies every renormalized coupling by c. This operation corresponds to a rescaling of the vertical axis of each graph in Fig. 4.
We find that U (l) and V (l) converge to finite positive values and Γ(l) and Ξ(l) vanish in the limit l → ∞. The coupling constant U (l) asymptotically behaves as where U Γ(0),Ξ(0) := lim l→∞ U (l) is the limiting value of U (l) determined by the bare couplings Γ(0) and Ξ(0) [recall that V (0) = 0]. Figure 4 shows the behavior of U (l) for several values of U (0) including U Γ(0),Ξ(0) . It turns out that the limit of U (l) is This indicates that the fixed points of the present RG transformation are (U, Γ, V, Ξ) = (0, 0, V, 0), (±∞, 0, V, 0) with V being an arbitrary positive value. Indeed, we can confirm that these points are fixed points by finding that (11) and substituting these values into the RG equations (9).
To gain a clearer understanding of the behavior of the coupling constants, we depict the RG flow diagrams in Fig. 5. We see that the coupling constants approach the following values in the limit l → ∞: where lim l→∞ V (l) is a finite positive value. The fixed points (U, Γ, V, Ξ) = (±∞, 0, V, 0) (V > 0) are stable (i.e., any coupling constant is irrelevant around these fixed points 3 ) and these correspond to the two phases: On the other hand, the fixed point Around this fixed point, U is relevant and the other coupling constants are irrelevant. 3 Recalling that the scaling factor of our RG transformation is b = 3 and U (l) is asymptotically proportional to 3 l as shown in Eq. (14), we find that the scaling dimension of U around the unsta-

Phase diagram
We now obtain the phase diagram of the frustrated Ising ladder in the dimer limit. In the following, we denote the bare couplings by U , Γ, and Ξ [not U (0), Γ(0), and Ξ(0)] for notational simplicity. Although there are three couplings U , Γ, and Ξ, it suffices to plot the phase diagram in the Ξ/U -Γ/U plane because dividing these couplings by U > 0 does not change the phase [note, as already remarked above, that U cΓ,cΞ = cU Γ,Ξ for any c > 0 and thus lim l→∞ U (l) is unchanged under multiplication of every bare coupling by c]. We can draw the phase diagram by marking the critical points  We can deduce from scaling theory [65], which is closely related to RG theory, that the phase boundary shown in Fig. 6 is of first order. The fact that the scaling dimension of the longitudinal field on the bottom row y U is equal to the spatial dimensionality of the system d = 1 indicates that the phase transition is of first order, because the correlation function on the bottom row at the critical point does not decay: (17) where · · · U =U Γ,Ξ denotes the expectation value at the critical point U = U Γ,Ξ and the approximation is a consequence of scaling theory [65]. This behavior of the correlation function leads to the value of the anomalous dimension η, one of the critical exponents defined by G bi,bj ∼ |i − j| 2−d−η : Recall that in general the scaling law for a quantum system is G bi,bj ∼ |i − j| −2(d+z−y U ) with z being the dynamic critical exponent [66] (which should not be confused with the variational parameter z) instead of Eq. (17). However, we find that z = 0 for the present system, because under a sufficiently large number of RG transformations the system becomes (almost) classical, by elimination of the transverse field Γ and the XX interaction Ξ.
Consequently, Fig. 6 demonstrates that there is a firstorder phase boundary that completely separates the staggered phase from the symmetric phase on the Ξ/U -Γ/U plane. In other words, the first-order phase transition encountered during quantum annealing cannot be removed by stoquastic or nonstoquastic XX catalysts.
Note that whether the bare Hamiltonian (1) is nonstoquastic depends not only on Ξ (= Ξ tt ) but also on Γ (= Γ t ), Γ b , Ξ bb , and Ξ tb [recall that the Hamiltonian is nonstoquastic if the fields and interactions satisfy Eq. (3), whose third condition |U/2|, |Ξ tb | < K holds in the dimer limit K U, Γ a , Ξ aa ]. For example, consider the case of sgn Γ b = sgn Γ and sgn Ξ bb , sgn Ξ tb ∈ {0, sgn Ξ}, where the sign function is defined as sgn x = x/|x| for x = 0 and sgn x = 0 for x = 0. Then, the Hamiltonian is stoquastic if Γ = 0 ∨ Ξ ≥ 0 and nonstoquastic if Γ > 0 ∧ Ξ < 0. Stoquasticity for Γ = 0 ∧ Ξ < 0 follows from the curing transformation that flips the signs of Ξ tt , Ξ bb , and Ξ tb , as pointed out in Sec. II A. The existence of this transformation indicates that the critical points on Γ = 0, U Γ=0,Ξ , are invariant under a change of the sign of Ξ.
It is noteworthy that the transition point derived from our real-space RG method in the absence of XX interactions, Γ/U = Γ/U Γ,Ξ=0 = 1.9314900, is not far from that obtained by numerical diagonalization of the quantum dimer model (whose Hamiltonian is given by Eq. (7) with Ξ tt = 0), Γ/U ≈ 1/0.6 = 1.66 · · · [11], in spite of the fact that our RG analysis is not an exact method.

B. Small frustration limit of the Ising chain
In this section we step back from the full Ising ladder and analyze the phase transition in the limit of small frustration U K, Γ a , Ξ aa at zero temperature. Since the bottom spins are fixed to down in this limit, we obtain an antiferromagnetic Ising chain with a transverse field and XX interactions as an effective model: Our purpose in this section is to reveal how the XX interactions affect the second-order transition that appears in the case of no XX interaction and large U (see Fig. 2 or Ref. [11]). Readers who are interested only in the effects of the XX interactions on the first-order transition can proceed to the next section on the DMRG calculations. Similar analyses were conducted by Langari [67,68]. He carried out the real-space RG analysis of the XXZ chain in a magnetic field [67], whose Hamiltonian is equivalent to Eq. (19) with i KŶ t,i ,Ŷ t,i+1 added, and derived the zero-temperature phase diagram. In addition, he obtained the zero-temperature phase diagram of the model given by the Hamiltonian (19) for Ξ tt ≤ 0 using the real-space RG method [68].
To compare the effects of ferromagnetic and antiferromagnetic XX interactions, we perform the real-space RG analysis for positive and negative Ξ tt . Instead of the Hamiltonian (19), we consider the ferromagnetic Ising chain with a transverse field and XX interactions, which is obtained by the gauge transformation for every odd i (we omitted the subscript t in Γ t , Ξ tt , and the Pauli operators for notational simplicity). The real-space RG analysis of the model (20) proceeds in the standard manner [65], which is detailed in Appendix C. We first separate the Hamiltonian into intrablock and interblock Hamiltonians after partitioning the chain every two sites, which is valid when Ξ is not a negative large value (if Ξ < 0 and |Ξ| is large, block partitioning every odd number of sites will be needed to retain the antiferromagnetic order in the x direction). Then, we diagonalize the intrablock Hamiltonian and project the Hilbert space onto the two-dimensional low-energy subspace in each block. Our real-space RG transformation is slightly different from that in Ref. [68], in the sense that the intrablock Hamiltonian in our transformation includes the magnetic field only at the left site in each block while that in Ref. [68] includes the fields at both sites. An advantage of our scheme is that the critical point Γ = Γ c and the critical exponent ν for Ξ = 0 coincide with the exact results Γ c = K and ν = 1, where ν is defined by the correlation length L corr ∼ |Γ − Γ c | −ν near the critical point [65].
Equation (21) indicates that the XX interaction is irrelevant in the limit U → ∞ unless the bare magnitude |Ξ| is large. This consequence is consistent with the result in Ref. [68], which predicts that the antiferromagnetic XX interaction gradually disappears when its bare magnitude is not large. It is also known for the model (19)  with the antiferromagnetic Y Y interactions of the magnitude K that Ξ is irrelevant for |Ξ/K| ≤ 1 regardless of the sign of Ξ [67].
The critical points are determined by the equation for the bare couplings γ and ξ. Note that all the critical points but (γ, ξ) = (0, 1) are inside the region |ξ| < 1 + γ 2 , in which our analysis can be applied. Since the RG equations (21) imply that the antiferromagnetic Ising chain with the transverse field and the XX interactions belongs to the same universality class as the transverse-field Ising chain, there will be secondorder transitions at the critical points (22). We show the phase diagram in Fig. 7, which indicates that the ±XXtype catalysts do not eliminate the second-order phase transition encountered during quantum annealing. The phase diagram is in qualitative agreement with the diagrams obtained by Langari [67,68], though the model in Ref. [68] covers only the case of the antiferromagnetic XX interactions and that in Ref. [67] includes the Y Y interactions that have the same sign and magnitude as the ZZ interactions.

IV. DMRG CALCULATION
We next employ the DMRG method [69][70][71][72] to obtain the phase diagram with a moderate coupling magnitude of K/U in Sec. IV A, which supplements the real-space RG analysis. In addition, we compute the energy gap of the finite-size system with a nonstoquastic XX catalyst in Sec. IV B.

A. Phase diagram
We perform the DMRG calculations for the frustrated Ising ladder with the transverse fields Γ t = Γ b = Γ and the XX interactions Ξ tt = Ξ and Ξ bb = Ξ tb = 0 in Eq. (1).
To avoid trapping of the calculations in one of energy local minima, we start with several sweeps by taking not only the ground state but also the first excited state as the target states in the finite DMRG procedure, followed by the single target DMRG procedure to obtain the convergence of the ground-state energy. The truncation number m of the density-matrix eigenvalues in the DMRG calculations is set as large as m = 200 to achieve a maximal truncation error less than 10 −12 .
We first show the ground-state phase diagram for K/U = 5 and L = 20 in Fig. 8. 4 Since the third condition of Eq. (3) is satisfied, the Hamiltonian is stoquastic for Γ = 0 ∨ Ξ ≥ 0 and nonstoquastic for Γ > 0 ∧ Ξ < 0. It is significant that the phase diagram obtained by the DMRG method has a similar shape to that for K → ∞ in the thermodynamic limit L → ∞ predicted by the realspace RG method (see Fig. 6 for the latter diagram). As explained below, the phase boundary between the staggered phase and the symmetric phase in Fig. 8 is characterized by first-order transitions.
In order to determine the phase boundary in Fig. 8, we calculate the staggered order parameter for spins on the top row defined as where · · · denotes an average over the ground-state wave function. Figure 9(a) shows the results for K/U = 5 and L = 20 as a function of Γ/U with several representative values of Ξ/U . When Ξ/U = −0.4, the staggered order parameter S remains approximately 1 for Γ/U up to ∼ 2.2, indicating that the ground state is in the staggered phase, and then suddenly decreases to almost zero (S ≈ 0.08) as Γ/U increases further, which corresponds to a first-order transition to the symmetric phase. A similar behavior continues until Ξ/U ≈ −2.65 with varying Γ/U as shown in Fig. 9(a), e.g., for Ξ/U = −2.6, where the first-order transition occurs at Γ/U ≈ 1.1. On the other hand, when Ξ/U = −3.0, the staggered order parameter S gradually increases from S ≈ 0.08 with Γ/U , implying that the ground state remains in the symmetric phase. Indeed, as shown in Fig. 9(b), the first-order transition from the staggered phase to the symmetric phase occurs at Ξ/U ≈ ±2.65 when Γ = 0.
To support these results, we calculate the first derivatives of the ground-state energy E 0 = Ĥ with respect to Γ and Ξ since a first-order transition is characterized by a point where ∂E 0 /∂Γ or ∂E 0 /∂Ξ is discontinuous. According to the Hellmann-Feynman theorem [73], the first derivatives of the ground-state energy E 0 with respect to Γ and Ξ are calculated as As shown in Fig. 10, the first derivatives of the groundstate energy E 0 with respect to Γ and Ξ exhibit discontinuities exactly at the points where the staggered order parameter S changes abruptly in Fig. 9, confirming the first-order nature of the transition.

B. Energy gap
We next obtain the energy gap between the ground state and the first excited state of the finite-size system with a nonstoquastic XX catalyst using the DMRG method. We adopt the multi-target DMRG procedure with the ground state as well as the two lowest excited states as the target states, taking the truncation number m as large as 2,000. We show in Fig. 11 the finitesize scaling of the minimum gap ∆(L) through the firstorder transition driven by varying Γ/U at K/U = 5 for Ξ/U = 0, −1, −2. We find that the gap ∆(L) decreases exponentially with L, i.e., ln ∆(L) ∼ −αL, for all three values of Ξ/U (for Ξ/U = 0 and L ≤ 12, we reproduce the previously reported results in Ref. [11]). It is noteworthy that the absolute value of the slope, α, in ln ∆(L) tends to be smaller in the presence of negative Ξ. From the point of view of quantum annealing and optimization this implies a quantitative, if not qualitative, improvement by the nonstoquastic XX catalyst, since a larger gap implies a faster time-to-solution by the adiabatic theorem [6,7].
We now discuss why the nonstoquastic XX catalyst reduces the decay rate α of the minimum gap ∆(L). Consider the sum of the terms with the coefficients U , Γ, and Ξ in the HamiltonianĤ as a perturbation. It may be expected that the ground state and the first excited state at the transition are superpositions of columnar and staggered configurations and ln ∆(L) is roughly proportional to L ln(Γ/K), because all the bottom spins need to be flipped by the bottom transverse field Γ b = Γ to move from the symmetric phase to the staggered phase (note that the present XX interactions Ξ are applied only along the top row and cannot flip the bottom spins). Indeed, log 10 (1.85/5) = −0.43 and log 10 (2.85/5) = −0.24 are close to the slopes of the lines for Ξ/U = 0, −2 in Fig. 11, respectively, where Γ/U = 1.85 is the transition point for Ξ/U = 0 and Γ/U = 2.85 is for Ξ/U = −2 when K/U = 5 (see Fig. 8). This is consistent with the fact that the decay rate α is proportional to ln(K/Γ) when K/U is varied in the absence of the XX catalyst [11]. Our argument suggests that the nonstoquastic XX catalyst on the top row with an appropriate magnitude Ξ increases the transverse field Γ at the transition and reduces the exponential decay rate of the gap in the present system. Although we have not carried out a numerical calculation of the gap of the model with the stoquastic XX catalyst due to the associated heavy computational cost, it is implied that the decay rate α of the minimum gap ∆(L) becomes larger for Ξ > 0 (as well as for Ξ/U −2.5) than for Ξ = 0, because Γ/U at the transition is smaller in the former case, as shown in Fig. 8.
Similar arguments were presented for a geometrically local Ising model on two connected rings [21]. There it was shown by numerical diagonalization that the nonstoquastic XX catalyst makes the minimum gap larger than in the case without the XX catalyst and in the case with the stoquastic XX catalyst. It was argued that one of the possible reasons for the softening of the avoided level crossing with the nonstoquastic catalyst is that the driver with the nonstoquastic XX catalyst causes the level crossing earlier in quantum annealing (which means a larger transverse field) than with the stoquastic XX catalyst, which is a consequence of perturbation theory.

V. SUMMARY AND CONCLUSIONS
We have studied the effects of stoquastic and nonstoquastic catalysts, implemented via ferromagnetic and antiferromagnetic XX interactions, in the setting of the frustrated Ising ladder. This model -without the XX interactions -is known to have a first-order phase transition with an exponentially decaying energy gap, which is characterized by a change in the topology of dimer configurations in the limit of strong frustration K → ∞ [11]. We have formulated a real-space RG transformation such that the symmetry of the problem is preserved and used it to obtain the phase diagram in the presence of XX interactions of both signs, stoquastic and nonstoquastic. The result shows that the first-order transition persists in the presence of XX interactions of moderate magnitude. The transition point obtained by the real-space RG method in the case without XX interactions is close to the value obtained by numerical diagonalization [11]. This is surprising, given that the real-space RG approach involves a number of uncontrolled approximations. In addition, we applied the real-space RG method to the case with small frustration and found that the second-order transition persists under the influence of XX interactions of moderate magnitude. We next performed extensive numerical computations by the DMRG method for a large but finite value of K in order to directly study the behavior of various physical quantities. The results for the order parameter and derivatives of the ground-state energy clearly indicate the existence of first-order phase transitions in the presence of stoquastic or nonstoquastic XX interactions, which is consistent with the conclusion from the real-space RG study. The structure of the phase diagram qualitatively and even semi-quantitatively resembles the one obtained by the real-space RG analysis. Our DMRG results furthermore confirm that the energy gap decays exponentially as a function of system size at the first-order transition points, both with or without XX interactions on the top row of the ladder. However, the nonstoquastic XX catalyst reduces the decay rate of the gap, as we showed numerically. We pointed out that this softening of the exponential decay of the gap may be due to an increase of the transverse field at the transition in the presence of a nonstoquastic catalyst.
A first-order phase transition is a sudden change of the system state between very different phases, e.g., between water and ice, and is unlikely to be induced or reduced by a series of gradual local changes. In the case of spin systems, the latter gradual change is exemplified by the introduction of XX interactions, which change the state of the system in the computational basis by flipping only pairs of spins simultaneously. Nevertheless, there exist examples in which nonstoquastic XX interactions change the order of a phase transition from first to second [15,16,18,[20][21][22], meaning a drastic reduction of the "strength" of the phase transition by nonstoquastic XX interactions, although counterexamples abound [24].
The first-order transition in the problem of the frustrated Ising ladder without XX interactions belongs to a more stable class in the sense that the two phases are separated by different topological structures in the limit of strong frustration. These topological structures are generated by the columnar and staggered configurations, which have magnetizations with opposite signs on the bottom row of the ladder. In the presence of frustration of infinite magnitude, the addition of any local operators of finite magnitude does not allow transitions between the topologically distinct states. It is therefore expected that the first-order transition is stable against the introduction of the XX interactions of either sign.
From the perspective of quantum annealing, the persistence of the first-order transition means that the computational complexity of the problem, exponential in the system size, remains intact under the introduction of stoquastic or nonstoquastic XX catalysts as long as the system evolves under adiabatic unitary dynamics. It is an important and interesting open question to study whether any of these conclusions are modified under nonunitary (open-system) dynamics or diabatic evolution, both of which take place in real quantum devices [47]. This appendix displays additional data to supplement the phase diagram of the Ising ladder without XX terms in Fig. 2. As mentioned in Sec. II B, it appears that the minimum gap point (namely, the location at which the energy gap between the ground state and the first excited state is minimized for a fixed K/U ) has a discontinuity at K/U ≈ 1.5 between Γ/U ≈ 1.5 and 1.9. To examine whether this phenomenon is a finite-size effect, we plot the minimum gap points for L = 4, 6, 8, 10 in Fig. 12. We find near convergence already for L = 10, which suggests that the discontinuity of the minimum gap point location is not a finite size effect.
An explanation is provided in Fig. 13, which shows the energy gap ∆E as a function of Γ/U for L = 10 and K/U = 1.25, 1.49, 1.5, 1.75. It turns out that the energy gap for K/U ≈ 1.5 takes a "double-well" form, which is the origin of the discontinuity of the minimum gap point. The energy gap ∆E is minimized in the left "well" for K/U 1.5, while ∆E is minimized in the right "well" for K/U 1.5. We leave a detailed understanding of this phenomenon as a future topic of research.  at different values of K/U . The inset is a magnification of the lines for K/U = 1.49, 1.5, which shows that the global minimum of ∆E moves from the left "well" to the right "well" as K/U changes.

Appendix B: Derivation of RG equations in the limit of large frustration
We derive the RG equations of the frustrated Ising ladder in the limit of large frustration (i.e., the dimer limit), which were used in Sec. III A. Appendix B 1 gives a review of basic properties of the model. In Appendix B 2, we define a Hamiltonian that appears in the RG analysis. In Appendices B 3-B 7, we explain the way to construct the real-space RG transformation including the variational ansatz in detail and write down the RG equations. We also calculate the renormalized couplings and show their behavior.

Preliminary analysis
To prepare for the real-space RG analysis, we demonstrate why the columnar and staggered configurations are the low-energy states of the frustrated Ising ladder (1) in the dimer limit K → ∞ and derive the number of low-energy states, as indicated in Ref. [11].
Since the overall energy scale does not change the statistical properties in the zero-temperature limit, we consider the dimensionless Hamiltonianĥ =Ĥ/K: where u = U/K, γ a = Γ a /K, and ξ aa = Ξ aa /K are sufficiently small dimensionless parameters. We assume that u, γ a , and ξ aa are of the same order δ 1. We separate the Hamiltonianĥ Using translational invariance yields iẐ t,i = i (Ẑ t,i +Ẑ t,i+1 )/2 and iẐ t,iẐb,i = i (Ẑ t,iẐb,i +Ẑ t,i+1Ẑb,i+1 )/2. We can thus rewrite the zeroth-order Hamiltonianĥ (0) aŝ Next, let us rewrite the same Hamiltonian as a linear combination of projectors: .
Here, each state vector containing an array of arrows means a product state and the positions in the array correspond to those in the ladder: We denoted the normalized eigenstates ofẐ a,i with the eigenvalues +1 and −1 by |↑ a,i and |↓ a,i , respectively.

Generalized Hamiltonian
The bare dimensionless Hamiltonian is given by Eqs. (B2a), (B6), and (B2c). Now we define a more general Hamiltonian, because an RG transformation will produce interactions not included in the bare Hamiltonian: . The length L is even and position L + 1 is identified with 1 due to the assumption of periodic boundary conditions. We denoted an arbitrary operator of O(δ 1 ) byô and the projection operator onto (1) will disappear after an RG transformation and is irrelevant in the sense of RG theory. The bare Hamiltonian can be obtained by assigning the coefficients in Eq. (B6) to k(σ i , σ i+1 ) and setting γ = γ t , v = 0, and ξ = ξ tt .
In Eq. (B14b), we denoted configurations of two spins at each position i by As a first step of an RG transformation, we partition the ladder into L = L/b blocks of spins, where b is the scaling factor. We project the Hilbert space onto a four-dimensional space for each block. The projector is written aŝ where |σ = L I=1 |σ I I for σ = (σ 1 , . . . , σ L ) and { |σ I I } σ I ∈H1 is a set of four orthonormal states in the Ith block constituted by the 2b sites (a, i) (a = t, b; i = b(I −1)+1, . . . , bI). In the present system, the scaling factor b should be an odd number to prevent the RG transformation from breaking the antiferromagnetic order in the staggered phase (an RG transformation with even b would turn antiferromagnetic order into ferromagnetic order). We take b = 3 in the following.
We expand the states |σ I I (σ I ∈ H 1 ) up to first order in powers of δ: These states satisfy the orthonormality which yields the constraint at each order: The block-product state has the zeroth-and first-order parts

Variational ansatz for the projector
The renormalized Hamiltonian is given byˆ h =QĥQ. Now we need to construct the projectorQ. In the standard real-space RG method, the Hilbert space would be projected onto a low-energy space of the intrablock Hamiltonian after separating the Hamiltonian into intrablock and interblock Hamiltonians [65]. However, we do not adopt this method, because the interblock interactions are strong and diagonalization of the intrablock Hamiltonian will not yield a low-energy space of the entire system correctly (recall that the interblock operators are of O(δ 0 )). In order to find a low-energy space not of an intrablock Hamiltonian but of the whole Hamiltonianĥ, we propose a new RG procedure in which the projectorQ is variationally determined.
Let us construct the zeroth-order projectorQ (0) = σ∈H L |σ . It is sufficient to take into account only the penalty Hamiltonianĥ (0) to determine the forms of |σ I I . The construction proceeds as follows: . This assumption is motivated by the expectation that the dimer structure will be kept around the fixed point that dominates the phase transition as the RG transformation is applied many times (we need to explore a neighborhood of a fixed point to derive critical properties). The dimer structure means that the set of low-energy states has a one-to-one correspondence with the set of dimer coverings on a two-leg ladder. We can preserve the dimer structure by constructing a zeroth-order renormalized Hamiltonian that imposes a penalty on |σ I σ I+1 (0) I+1 for (σ I , σ I+1 ) ∈ H 2 \ D 2 and no penalty for from paying a penalty, ↑ ↑ ↑ ↓ (0) ↓ ↓ (0) where the parameters α 1 , α 2 , z, β 1 , β 2 , β 3 ∈ C satisfy the normalization conditions: Equation (B23) can be regarded as a variational ansatz for the projectorQ (0) . We will determine the variational parameters α 1 , α 2 , z, β 1 , β 2 , and β 3 in Appendix B 6.