Counterdiabatic Hamiltonians for multistate Landau-Zener problem

We study the Landau-Zener transitions generalized to multistate systems. Based on the work by Sinitsyn et al. [Phys. Rev. Lett. 120, 190402 (2018)], we introduce the auxiliary Hamiltonians that are interpreted as the counterdiabatic terms. We find that the counterdiabatic Hamiltonians satisfy the zero curvature condition. The general structures of the auxiliary Hamiltonians are studied in detail and the time-evolution operator is evaluated by using a deformation of the integration contour and asymptotic forms of the auxiliary Hamiltonians. For several spin models with transverse field, we calculate the transition probability between the initial and final ground states and find that the method is useful to study nonadiabatic regime.


Introduction
The Landau-Zener (LZ) formula is an analytic solution for a dynamical two-level quantum system [1,2]. The system Hamiltonian has a linear-time dependence and the time evolution of the state is exactly solvable. By examining together with some generalized results [3,4], we now have a firm picture of nonadiabatic transitions between two instantaneous eigenstates in two-level time-dependent Hamiltonians. This result has a great influence on various fields of physics and we can find a vast amount of applications in literature.
Since the energy-level crossing occurs in general quantum systems, we expect that a similar picture holds for multistate systems where the time-dependent Hamiltonian takes the form with arbitrary operatorsẐ andX. Under this Hamiltonian with some minor conditions, the transition probability from the ground state at t → −∞ to the highest state at t → ∞ was conjectured by Brundobler and Elser (BE) [5]. The following studies revealed that the formula is indeed correct and some extensions are possible [6][7][8][9][10][11][12]. The formula can be interpreted as a natural generalization of the original LZ formula. The same type of the Hamiltonian is used in the method of quantum annealing [13][14][15][16][17]. To obtain the optimized solution, we treat Ising Hamiltonians with a transverse field whose time dependence is linear in the standard protocol. Starting from the ground state of the Hamiltonian at t = −∞ where the transverse-field termX is dominant, we consider the time evolution until the Ising termẐ becomes dominant at t = 0. Although the standard LZ protocol considers the time evolution until t goes to infinity, we treat the same Schrödinger equation and expect that the LZ picture is still applicable [18]. The field is changed very slowly and the adiabatic picture is applied except around the crossing point. Although this adiabatic approximation gives a reasonable picture, the dynamical description is required to calculate the transition probability, which will be useful to estimate the efficiency of the protocol.
In the optimization problem, we are basically interested in the probability that the system remains the ground state. Unfortunately, the BE formula cannot be applied to this situation. As we discuss below, the description of the transition from the ground state to the highest state requires few information about the system in contrast to describing the transitions to the other states.
The standard form of the two-level LZ Hamiltonian is written asĤ(t) = J 2σ z + γt 2σ x wherê σ z andσ x are Pauli matrices. Then, it can be shown that the time-evolution operator from t = −T to t = T is decomposed asÛ at T → ∞.Û ±T are diagonal inσ x basis and are rapidly-oscillating as functions of T . We are not interested in these phase contributions and the transition probabilities from eigenstates at t = −T → −∞ to eigenstates at t = T → ∞ are determined byÛ 0 only.Û 0 is independent of T and is given byÛ where This is the LZ formula. This decomposition of the unitary operator looks reasonable since the nonadiabatic transition only occurs around t = 0. The time domain is decomposed into adiabatic domain |t| ≫ J/γ, which is described byÛ ±T , and diabatic domain |t| < J/γ, described byÛ 0 . Then, to calculate the transition probability we may focus only onÛ 0 . This operator is independent of T and is analytically tractable. It is a fundamental question whether such a decomposition is generally possible.
Recently, an interesting approach was proposed in Ref. [19]. The Hamiltonian contains several parameters and they are treated as dynamical variables. Using the introduced auxiliary Hamiltonians, we can deform the one-dimensional path −T → T of the time evolution to a multi-dimensional one in parameter space. The auxiliary Hamiltonian plays a role of the generator for the corresponding parameter. Several solvable models were treated to construct the auxiliary Hamiltonians and to demonstrate that the method is useful to describe the time evolution [19][20][21].
In this paper, we use the method of the auxiliary Hamiltonians to demonstrate the separation of the unitary operator as in Eq. (2). We also show that the auxiliary Hamiltonians are closely related to the counterdiabatic terms discussed in the method of shortcuts to adiabaticity (STA) [22][23][24][25][26][27]. By using the auxiliary Hamiltonians, we evaluate the time-evolution operator for several spin models with transverse field to calculate the transition probabilities.
The organization of this paper is as follows. In Sec. 2, we study the general structure of time-dependent Hamiltonian. Combining the idea of Ref. [19] with the result from STA, we find a general form of the auxiliary Hamiltonians. The method is applied to the Landau-Zener Hamiltonian. In Sec. 3, we demonstrate the decomposition of the time-evolution operator in Eq. (2). Based on the general results, we numerically calculate the transition probability for several spin models in Sec. 4. The final section 5 is devoted to discussions and conclusions.

Auxiliary Hamiltonian for parameter generation
We consider a general time-dependent HamiltonianĤ(t, ω) with a set of parameters ω = (ω 1 , ω 2 , . . . ). The Schrödinger equation for the state |ψ(t, ω) is written as In Ref. [19], the derivative with respect to one of the parameters was considered to define the auxiliary Hamiltonian as The auxiliary HamiltonianĤ j (t, ω) plays the role of generator for the corresponding parameter ω j .Ĥ andĤ i satisfy the consistency conditions The equations are interpreted as the zero curvature condition in noncommutative geometry.
It is generally a difficult task to find the auxiliary Hamiltonians for a givenĤ(t, ω). In the previous studies [19][20][21], the condition that the both sides of the above two equations are respectively equal to zero were examined to treat several solvable models. To examine more general cases, we show a relation betweenĤ andĤ j by using the method of STA.
A general time-dependent Hamiltonian is written aŝ where λ(t, ω) = (λ 1 (t, ω), λ 2 (t, ω), . . . ) is a set of time-dependent parameters. Each term of the Hamiltonian is given respectively bŷ ǫ n (λ) is real and {|n(λ) } represents a set of orthonormal basis. The solution of the Schrödinger equation is given by the adiabatic state ofĤ 0 (λ): where c n is a constant determined from the initial condition at t = t 0 and θ n is represented by the sum of dynamical phase and geometric phase: We assume that the initial parameter λ(t 0 , ω) = λ 0 is independent of ω. This separation of the Hamiltonian means that the state is written in terms of the eigenstates of the first term and the second term prevents nonadiabatic transitions. Although it is nontrivial to find this separation for a given total Hamiltonian, the separation is always possible in principle [28]. We also note that the form of λ(t, ω) cannot be obtained straightforwardly from the total Hamiltonian since the total Hamiltonian is dependent not only on λ(t, ω) but also on the time derivative of λ(t, ω). When the state is evolved under this Hamiltonian, the average energy at each time is determined byĤ 0 and the state is changed by the counterdiabatic termξ. This means that the Hamiltonian is separated according to its roles: energy and generator.ξ j corresponds to the generator for the parameter ω j , which implies a relation to the auxiliary HamiltonianĤ j . In fact, it is straightforward to show that the counterdiabatic terms satisfy the zero curvature condition This is obtained from the spectral representation in Eq. (11). We also obtain the form ofĤ i by taking the derivative of Eq. (12) with respect to ω j aŝ whereĤ We note that the structure of the Hamiltonians is insensitive to the choice of ǫ n and is mainly determined by the counterdiabatic terms. The simplest choice is to set ǫ n = 0 as we treat below. Then, the auxiliary Hamiltonians are completely equivalent to the counterdiabatic terms.

Landau-Zener Hamiltonian
In the following we focus on the Hamiltonian whereẐ andX are arbitrary operators. We take one of the parameters γ to be positive without losing generality. We also assume that the spectrum of the operatorX is nondegenerate so that we can apply the BE formula to this system [5]. WhenẐ andX are chosen as two of the Pauli operators, we obtain the standard LZ Hamiltonian. Here our motivation is to setẐ as an Ising model andX as the transverse field. We consider the case where the Hamiltonian is represented by a real symmetric matrix. Then, the form of the Hamiltonian in Eq. (9) implies thatĤ 0 = 0. This is because the Hamiltonian cannot be real symmetric when we considerĤ 0 andξ at the same time. Even ifĤ 0 is real symmetric, the definition of the counterdiabatic term in Eq. (11) inevitably introduce complex components to the Hamiltonian. We do not consider the case whereẐ andX commute with each other since this case describes a trivial situation. We also exclude the trivial contribution coming from a term proportional to the identity matrix. Then, the only possible choice is to setĤ 0 = 0 and to take the counterdiabatic term real symmetric. By restricting ourselves to the real symmetric Hamiltonian, we can identify what the corresponding parameter λ(t, ω) is. The Schrödinger equation is written as This equation clearly shows that the state is given by a function of √ γt and Jt. On the other hand, Eq. (12) shows that the state is a function of λ whenĤ 0 = 0. Combining these observations, we can identify the parameter λ as The Hamiltonian is decomposed aŝ We need two counterdiabatic terms. They are related to each other aŝ By using Eq. (14), we obtain Several remarks are made below: (i).ξ 1 andξ 2 cannot be directly obtained from the original Hamiltonian. These new operators are required when we consider the deformation of the path between initial and final Hamiltonians as we show in the following. (ii). Even ifĤ is real symmetric,ξ 1,2 , andĤ 1,2 , generally have complex elements. Previous studies only treated the case where the auxiliary Hamiltonians are real symmetric [19][20][21]. (iii). The auxiliary Hamiltonians in the present example are given byĤ j = tξ j with j = 1 and 2. (iv).Ĥ 1 and H 2 are not independent. We find this nontrivial property by using STA. It is convenient to rewrite the differential equation (22) as a function of λ = λ 1 and g = λ 2 /λ 1 . We putξ Then, we obtain whereẐ The differential equation is with respect to λ and not to g. This can be understood from the Schrödinger equation in Eq. (18) where g appears as a time-independent parameter. We also note that the unitary transformation represented by the exponential factor in Eq. The equation can be formally expressed aŝ with the use of the path-ordered unitary operator To avoid the divergence ofξ 2 (λ) for λ → 0, we require the conditionξ(0, g) = 0, which determines the constant of integration. We can also write a form expanded in terms of g: 3 Decomposition of the time-evolution operator

Path deformation
By introducing the auxiliary Hamiltonians, we can consider the path deformation of the time evolution. The unitary time-evolution operator represented by the time-ordered productŝ is equivalent to the path-ordered operator under the zero curvature condition in Eqs. (7) and (8). Here, we write ω µ = (t, ω 1 , ω 2 ), H µ = (Ĥ,Ĥ 1 ,Ĥ 2 ), and C represents a path in (t, ω)-space. We fix the initial point (t, ω) = (−T, √ γ, J) and final point (T, √ γ, J) and take the limit of T → ∞.
The path C can be chosen arbitrarily. Here we take the path where each segment represents the straight line, as we show in Fig. 1. The advantage of using this path is that the time variable takes a large constant value at the first and last segments, |t| = T , and the operator from the second segment gives a trivial contribution. We havê The second exponential factor comes from the second segment and gives the identity operator. By using Eq. (23) and using the propertyξ 2 (λ 1 , λ 2 ) =ξ 2 (−λ 1 , λ 2 ), we can writê where δ = J/ √ γ. Then, we only need to know the asymptotic behavior ofξ(λ, g) at λ → ∞.

Asymptotic form
As we can understand from Eq. (26),ξ diverges linearly in λ at most. The linear divergence is obtained whenẐ has diagonal elements inX basis. We also generally find thatξ has a logarithmic-divergent term. We writeẐ =Ẑ 0 +Ẑ 1 whereẐ 0 represents the diagonal part andẐ 1 offdiagonal part. At zeroth order in g, the diagonal part ofξ is linear in λ. The jth element is given by The logarithmic divergence comes from the diagonal part at the first order in g. The jth element is given by where X j represents the eigenvalue ofX. The estimation of the integral is given in Appendix A. It should be noticed that these divergent terms are obtained from the dynamical phase calculated in the adiabatic approximation [8] φ j (T ) = T dt ǫ j (t).
The instantaneous energy ǫ j (t) is evaluated perturbatively and we obtain up to second order in J as Then, we have The first term corresponds to the phase in Eq. (23), the second to Eq. (34), and the third to Eq. (35). The key observation in Ref. [8] is that these are only divergent contributions to the phase at T → ∞. In Appendices A and B, we show that the same property holds forξ. The divergent terms inξ come from the diagonal part as we estimated above. In the following, we show that they can be represented as a phase in the time-evolution operator and do not give any contribution to the transition probabilities.
Using the arbitrariness ofφ, we can remove the divergent contribution from the diagonal part ofξ asξ where the subscript 0 denotes the diagonal part and 1 denotes the offdiagonal part. It is possible to chooseφ such that the conditionξ 0 = 0 is satisfied. The addition of phase in ξ 1 does not induce any additional divergence. Thus, we conclude that the time-evolution operator is decomposed aŝ as is expected in Eq. (2). To calculate transition probabilities between the eigenstates of X, we may focus onÛ 0 which is only dependent on δ = J/ √ γ and contains no divergent contributions.

Transition probability
It is still a difficult problem to find the explicit form ofξ in general systems. Even if we can findξ, we have to manage the path-ordered integral in Eq. (43). Although these calculations look a formidable task, the BE formula implies that the unitary operator is well-approximated byξ at zeroth order in g:Û The formula treats the transition probability from the lowest state at t → −∞ to the highest state at t → ∞ (and vice versa). Their states are represented by the same eigenstate ofX and are denoted by |0 . The transition probability is evaluated by the second-order cumulant expansion as where X 0 represents the eigenvalue ofX for |0 . This is exactly the BE formula [5].
In most of the applications we are interested in the transition probability between the ground states at t → −∞ and at t → ∞. Below, for several spin models, we study the transition probabilities by using the truncation approximation ofξ(∞, g) keeping the lowerorder terms in g. The approximation is valid for nonadiabatic regime where J < √ γ and can be a complementary method of the adiabatic approximation.

Two-body interacting spin model
First, we study the two-body interacting Ising model with the transverse field: S = (Ŝ x ,Ŝ y ,Ŝ z ) represents the total spin of the spin-1/2 systemŜ = 1 2 N i=1σ i where N represents the number of 1/2-spin. The Hamiltonian commutes withŜ 2 and we consider the block with the quantum numberŜ The ground state and highest state are contained in this block and the Hamiltonian is represented by a (N + 1) × (N + 1) matrix. We start the time evolution from the ground state of the initial Hamiltonian and examine the transition probabilities at large T . In this model, the BE formula is evaluated as We note that this Hamiltonian has diagonal elements ofẐ inX basis. In the present model, the decomposition ofẐ =Ẑ 0 +Ẑ 1 is written aŝ To examine the effect by the diagonal term, we also treat a modified Hamiltonian with the diagonal part ofẐ removed: We note that the BE formula in Eq. (48) is unchanged by this modification. We use two kinds of approximation. The first one is to keep only the zeroth-order term in Eq. (28). The time-evolution operator is approximated as Eq. (44). The second approximation is to use some ansatz taking into account the higher-order corrections. We put whereξ 0 represents the zeroth-order term as in Eq. (44) andŜ ± =Ŝ y ±iŜ z . The operator form of the double-commutator term is taken from the second term of the expansion in Eq. (28) with some modifications by the gauge transformation. Instead of taking into account the path-ordered integral, we put and complex parameters a and b are determined by the method of least squares so that Eq. (51) reproduces the BE formula. Then, using the determined parameters, we numerically calculate the transition probability that the state remains the ground state. The zeroth-order approximation of the transition probability from the initial ground state to the final highest state is plotted in Fig. 2. We see that the approximation gives a reasonable description of the BE formula at small δ = J/ √ γ. We find a small deviation at large δ. In Fig. 3, we plot the transition probability to the ground state. The numerical calculations are compared with the results of the zeroth-order approximation and of the estimation by Eq. (51). We see that the numerical result for Eq. (50) is well described by theoretical estimations. A large deviation for Eq. (46) is considered to be the presence of the diagonal elements ofẐ 0 . Since the diagonal components do not contribute to the BE formula, it is impossible to estimate the undetermined parameters in the present method. In any case, the theoretical analysis works well in the nonadiabatic regime.

Three-body interacting spin model
For comparison, we study the three-body interacting Hamiltonian We note thatẐ 0 = 0 in this case. The BE formula is evaluated as This probability goes to unity at N → ∞. Then, compared to the two-body case, the groundstate probability takes a considerably smaller value. The comparison of the BE formula and the zeroth-order approximation is shown in Fig. 4. The ground-state probability is plotted in Fig. 5. We again find deviations at large δ, basically the same conclusion as the two-body case.

Size dependence
The BE formula shows that the size dependence of the probability to the highest state is totally different between the two-and three-body systems. We find that, at least for nonadiabatic regime, the probability to the ground state has the scaling Prob ∼ exp(−N c), as we show in Fig. 6. This scaling is expected to be changed in the superadiabatic regime δ ≫ 1 where the static picture holds. In the adiabatic approximation, the transition probability is determined by the energy gap between the ground state and the first-excited state [18]. The energy gap behaves differently between the two-and three-body systems as their models have different kinds of quantum phase transitions. In our numerical calculation, δ takes smaller values and we find that the scaling is unchanged in the present systems. The same scaling holds for probabilities to the other states except the highest state, which shows that the BE formula is specific to the highest state and does not describe a typical behavior.

Probability distribution
Finally, we examine the transition-probability distribution among the energy levels. The numerical results are shown in Fig. 7. For nonadiabatic regime δ ∼ 1, we find that, independently of the models, the transitions to higher levels are dominant. For adiabatic regime δ ≫ 1, the behavior is dependent of the models as we expect from the BE formula. We also see that our approximations work well in all the energy levels.

Conclusion
In conclusion, we have studied the multiple LZ transitions by using the method of the auxiliary Hamiltonians. We find that the auxiliary Hamiltonians introduced in Ref. [19] are equivalent to the counterdiabatic terms in the method of STA and are directly related to the original Hamiltonian. For real symmetric Hamiltonians, the time-dependent parameters are easily specified and we derive the differential equations for the counterdiabatic Hamiltonians. Using the counterdiabatic Hamiltonians and the deformation of the integration contour, we evaluated the time-evolution operator and calculated the transition probabilities. It is instructive to see that the general Hamiltonian contains several counterdiabatic terms. Even though the present study focused only on the LZ Hamiltonian, our method is general and is applicable to any Hamiltonians. The spectral representation of the counterdiabatic term in Eq. (11) is not so useful for complicated systems and we can use the differential equation such as Eq. (24) instead. To solve the equation, we can use some familiar methods such as perturbative expansions as we demonstrated in the present work.
As for the problem of the LZ transitions, our method becomes a generalization of the BE formula. The BE formula is applicable to a specific state only and it is hard to generalize the formula for different states. We treat the time-evolution operator, which allows us to estimate transition probabilities between arbitrary eigenstates. Although it is generally difficult to solve the differential equation in our method, a simple approximation gives a reasonable description in the nonadiabatic regime, which becomes a complementary method of the adiabatic approximation.
It is an interesting problem to apply the present method to more general systems. In the present study, we treated the following conditions for simplicity: (i). real symmetric Hamiltonians, (ii). nondegenerate eigenvalues ofX, and (iii). vanishing of diagonal elements ofẐ inX basis. We note the conditions (i) and (iii) are satisfied for the standard transversefield Ising models. Although the constraint (ii) does not hold in general, this problem can be circumvented by changing the basis. We can also apply our method to other dynamical systems with nonlinear time-dependent terms. Thus, our method is general enough and various applications can be found, which will be useful to understand general dynamical systems. Figure 6: Size dependence of the ground-state probability forĤ 2−body andĤ 2−body,0 (top) and forĤ 3−body (bottom). We take δ = 2.0. Legend is the same as in Fig. 3.

A Two-level system
We study the two-level system. By calculating the asymptotic form ofξ(λ, g) at λ → ∞, we demonstrate the decomposition of the time-evolution operator in Eq. (2).
In the two-dimensional case, we only have three independent operators except the identity operator. They are represented by the Pauli operators: We can also putξ (λ, g) = α(λ, g)X + β(λ, g)Ŷ + γ(λ, g)Ẑ, to write the differential equation (24) as ∂ λ β(λ, g) = sin λ 2 2 + gα(λ, g) cos ∂ λ γ(λ, g) = cos λ 2 2 − gα(λ, g) sin Using these equations, we can derive the differential equation for α as This is a linear nonhomogeneous differential equation and the solution is obtained by combining the particular solution and the general solutions of the corresponding homogeneous equation. We are only interested in the asymptotic form and the solution is basically expressed in power series of 1/λ. The particular solution is expanded in terms of g and we have α(λ, g) = g ln λ + 3 4λ 4 + · · · + g 3 1 − 2 ln λ 4λ 2 + · · · + · · · .
The general solution of the homogeneous equation has nonoscillating terms and oscillating ones: c(g), h(g), and θ(g) represent undetermined functions. Combining these results and neglecting irrelevant terms, we obtain α(λ, g) ∼ g ln λ + gc(g) − gh(g) λ cos Since the differential equation is in terms of λ, not of g, and only the asymptotic form is obtained in the present analysis, c 0 (g), h(g), and θ(g) are left undetermined. As we can understand from the differential equation, α(λ, g) is an odd function in g, which shows that c 0 (g), h(g), and θ(g) are even functions of g. We also see that this result shows that Eq. (35) holds in the two-level case. Since the integral in Eq. (35) is independent of the operatorẐ, the equation holds in the general case as well.
Although the last term in Eq. (63) is negligible at the asymptotic limit, this term is important to find β and γ. By using Eq. (57), we conclude that These are even functions of g. As we mentioned in the main body of the present paper, the divergence ofξ is logarithmic in λ and is contained only in the diagonal part. The log term in the trigonometric functions of the offdiagonal part can be removed by the gauge transformation as we show below. Now that we have obtained the asymptotic form ofξ, we make a gauge transformation represented byφ (λ, g) = g 2 2 ln λ + dg gc(g) X .
Then, the gauge-transformed operator is given bŷ This operator is independent of λ, which means that the divergent contribution can be removed by the gauge transformation. The time-evolution operator is decomposed aŝ where δ = J/ √ γ. Furthermore, the second line can be rewritten as This relation can be shown by using thatβ(g) andγ(g) are even functions of g. Now we have the decomposition in Eq. (2), as it should be. It is not necessary to show this decomposition in the present system because we know the exact solution. However, we see below that this demonstration is useful to apply the method to general systems.

B General properties of the auxiliary Hamiltonians
In Appendix A, we showed that the logarithmic divergence appears only in the diagonal part ofξ at first order in g. We can show that this property holds in general systems. For simplicity, here we only consider the case that the diagonal part ofẐ is zero:Ẑ 0 = 0. The general form of the diagonal part ofξ at first order in g is written in Eq. (35) and has a logarithmic contribution. The divergence is suppressed in the offdiagonal part of the same term because the integrand is more oscillating. We note that, to find this property, it is required that the spectrum ofX is nondegenerate. The same consideration is applied to the higher-order terms in g. At second order and third order, the diagonal part is respectively written as where X jk = X j − X k . Each term is oscillating rapidly and the integration basically gives a finite contribution. A possibly-divergent contribution arises for a specific choice of the subscript indices so that some cancellation occurs in the phase. For example, by choosing (j, k, ℓ, m) = (j, k, j, k) with j = k in Eq. (73), we have The integral is the same as that appears in two-level systems. We have already shown that this integral does not give any divergent contribution. We can generalize this consideration to the higher-order terms. Thus, we conclude that, in general systems, the divergence appears only in the diagonal part at first order in g.