November 17, 2023
We present the interesting resonance of two kinds of geometric quantities, namely the Aharonov-Anandan (AA) phase and the time-energy uncertainty, and reveal the relation between resonance and the hidden symmetry in the asymmetric Rabi model by numerical and analytical methods. By combining the counter-rotating hybridized rotating-wave method with time-dependent perturbation theory, we solve systematically the time evolution operator and then obtain the geometric phase of the Rabi model. In comparison with the numerically exact solutions, we find that the analytical results accurately describe the geometric quantities in a wide parameter space. We unveil the effect of the bias on the resonance of geometric quantities, (1) the positions of all harmonic resonances stemming from the shift of the Rabi frequency at the presence of the bias; (2) the occurrence of even order harmonic resonance due to the bias. When the driving frequency is equal to the subharmonics of the bias, the odd higher-order harmonic resonances disappear. Finally, the hidden symmetry has a resemblance to that of the quantum Rabi model with bias, which indicates the quasienergy spectra are similar to the energy spectra of the latter.
Since geometric phase of Berry’s seminal work was introduced in the cyclic evolution of a system under adiabatic condition[1], its effects have been discovered and experimentally measured in various fields of physics[2]–[11] as well as chemistry[12]. The geometric property of Berry’s phase visually lies in that it merely depends on the solid angle subtended by the closed path that the parameters traverse. Aharonov and Anandan (AA) extended Berry’s phase to non-adiabatic cases[13], [14], removing the adiabatic restriction and only assuming that the initial and final states of the quantum system differ by a total phase factor. Besides theoretical significance, AA phase has important applications to non-adiabatic geometric gates in quantum computation, with internal resiliency to certain noises and control errors[15]–[17]. Furthermore, Samuel and Bhandari developed the geometric phase in a general setting without the assumption of cyclic evolution[18]. All these phases are determined by the curve \(\mathcal{C}\) that the quantum state traverses in the projective Hilbert space \(\mathcal{P}\), and are thus defined as geometric quantities[19].
Time-energy uncertainty, as one of the most fundamental quantities in quantum mechanics, is also proved to be a geometric quantity[19], [20]. Although the proposal of the uncertainty principle for time and energy even dates back to the origin of quantum mechanics in 1920s, it is understood as setting a fundamental limit on the rate of quantum dynamics in its modern formulation[21]–[24], and thus widely employed in studies on the speed of gate operations in quantum computation[25]–[28], the precision of measurement in quantum metrology[29]–[32] and so on. With the rapid development of technologies in manipulating quantum systems, it has received increasing attention as the basis for optimal control theory[33]–[36].
The semiclassical Rabi model [37], [38], which describes a two-level system coupled with a classical monochromatic periodic field, is a typical quantum model exhibiting the properties of geometric quantities. Its studies have a rich history for both experimental and theoretical investigations[39]–[47], and nowadays this model is widely applied in quantum information technology[48]–[50]. It is found that there are a wide variety of interesting dynamical features in this model, such as well-known Rabi oscillations, Bloch-Siegert shifts, coherent destruction of tunneling (CDT)[51], driving-induced tunneling oscillations (DITO)[52], and plateau dynamics[53]. In this work, we consider the asymmetric semiclassical Rabi model, where a bias term is taken into account. The bias breaks the symmetry of the semiclassical Rabi model[54], which contributes to the complexity and phenomena different from those of the symmetric Rabi model. This allows for the representation of energy biases between the two states of a flux qubit[55].
While previous studies have primarily concentrated on the dynamics of the asymmetric semiclassical Rabi model[56], [57], including phenomena such as CDT and DITO, relatively little attention has been paid to the geometric quantities of this model. The geometric quantities, which are robust against certain types of noise and errors, have the potential to be powerful tools for quantum computation and quantum information processing. Therefore, our work not only investigates the properties of geometric quantities of the asymmetric semiclassical Rabi model but also provides insights into the behavior in a wide range of parameter space.
In this work, we investigate the harmonic resonance of AA phase and the time-energy uncertainty of the asymmetric semiclassical Rabi model, which is beyond the work on the unbiased case[58]. First, we generalize the harmonic resonance features of both geometric quantities in parameter space, and further explore the feature of the harmonic resonance. To derive the analytical expression for AA phase, we employ counterrotating-hybridized rotating-wave (CHRW) method, which takes account of the influences of counterrotating terms and bias together and is more reliable than rotating wave approximation. Further, considering the second harmonic terms of the Hamiltonian in the rotating frame, we combine perturbation theory with the CHRW method to calculate the AA phase and reveal its subtle resonant picture. Interestingly, at certain points in the parameter space, the odd harmonic resonance disappears in the numerical results. Using Floquet theory[39], [59], we explain this interesting phenomenon, and reveal the relation between the hidden symmetry of the asymmetric semiclassical Rabi model and odd harmonic resonance. While the hidden symmetry of the asymmetric quantum Rabi model has been extensively studied, our results provide the first instructive evidence for the similar phenomenon in the semiclassical counterpart.
The structure of this work is as follows. In Sec. 2, after a brief review of both geometric quantities, we perform numerical calculations to demonstrate the resonance phenomenon, and show the features of harmonic resonance. In Sec. 3, we employ the CHRW method to analytically calculate the AA phase and the positions of the higher-order harmonic resonances. In Sec. 4, we apply perturbation theory based on the CHRW method to take into account the effects of higher-order harmonic terms in calculating the AA phase in the higher-order harmonic resonance regime. In Sec. 5, we shed light on the absence of harmonic resonance in the results of numerical calculation as hidden symmetry appears. This is further discussed through comparison between the asymmetric semiclassical Rabi model and the asymmetric quantum Rabi model. Finally, we give the conclusion of this paper in Sec. 6.
The Hamiltonian of the asymmetric semiclassical Rabi model reads \[\label{Eq32classical32Rabi32H} H(t)= -\frac{\Delta}{2} \sigma_x - \frac{\epsilon+A\cos(\omega t)}{2} \sigma_z,\tag{1}\] where \(\sigma_x\) and \(\sigma_z\) are the Pauli matrices, \(\Delta\) is the tunneling strength, \(\epsilon\) is the static bias. \(A, \omega\) and \(T=2\pi/\omega\) are the amplitude, frequency and period of the driving field, respectively. We set \(\hbar = 1\) throughout this paper and also use \(\varepsilon(t)=\epsilon+A\cos(\omega t)\) to denote the bias-modulated driving field[44], [56], [60]. This Hamiltonian can represent the systems moving in an effective double well potential modulated by an ac field, under the condition that only the ground state in either well can be occupied[56], [61]. A system initially localized in one well will oscillate between the eigenstates in the left and right well due to quantum mechanical tunneling, which is exemplified by the two equivalent configurations of \(\mathrm{NH}_3\), as well as the two current states of a flux qubit[55], [60].
By implementing a rotation about the y-axis on this Hamiltonian, we obtain a transformed representation \(\exp (i\pi \sigma_y/4) H(t) \exp (-i\pi \sigma_y/4) =-\frac{\Delta}{2} \sigma_z + \frac{\epsilon+A\cos(\omega t)}{2} \sigma_x\). This particular form of the Hamiltonian is prevalently employed in the fields of quantum optics and nuclear magnetic resonance, where \(\Delta\) signifies the energy difference between the two levels, and the driving term causes transitions between these levels[62], [63].
In the following, we first present a succinct overview of the concepts and geometric properties of the AA phase and the time-energy uncertainty in Sec. 2.1 and Sec. 2.2. Then we display the \(\Delta\)-dependence of both geometric quantities using numerical methods in Sec. 2.3, and further discuss the features of resonance in Sec. 2.4.
Since the AA phase is a generalization of Berry phase, it is pertinent to briefly review the fundamental concepts and geometric properties of Berry’s phase.
Consider a time-dependent Hamiltonian \(H(\boldsymbol{R}(t))\) where \(\boldsymbol{R}=(R_1, R_2, ..., R_n)\) and \(\boldsymbol{R}(T)=\boldsymbol{R}(0)\), with the initial state being the \(m\)th eigenstate \(|m(0)\rangle\) of \(H(\boldsymbol{R}(0))\). If the system is subjected to adiabatic processes, i.e., \(|\langle j | \partial H/\partial t| k\rangle | \ll (E_j-E_k)^2, (j\ne k)\), where \(E_j\) is the \(j\)th eigenenergy of \(H(\boldsymbol{R}(t))\), the state of the system will remain in the \(m\)th eigenstate \(|m(t)\rangle\) while also obtaining a phase factor \(\exp \left[ i (\alpha_m(t) +\gamma_m(t) \right]\), where \(\alpha_m(t)=- \int_0^t E_m(t')dt'\) is called the dynamical phase, and \(\gamma_m(t)=i\int_0^t \langle m(t')|\dot{m}(t')\rangle dt'\) is called the adiabatic phase. Berry’s phase is defined as the adiabatic phase acquired over a cycle, i.e.,\(\gamma_m(T)\), which has been proved to be real and measurable.
The geometric property of Berry’s phase is corresponding to its relation to the solid angle subtended by the closed path in the parameter space spanned by \(\boldsymbol{R}\)[1]. Consequently, it is referred to as the geometric phase.
Based on Berry’s phase, AA removed the adiabatic restriction, while assuming that the initial and final states of the quantum system differ by a total phase factor \(e^{i\theta}\)[13]. Since the Hamiltonian in Eq. (1 ) is periodic, the cyclic state \(|\psi(t)\rangle\) satisfies \[|\psi(T)\rangle =e^{i\theta}|\psi(0)\rangle,\] where \(\theta\) and \(|\psi(0)\rangle\) are defined as the total phase and the cyclic initial state, respectively. Then the AA phase is defined by subtracting the dynamical phase \(\alpha\) from \(\theta\)[13]: \[\label{Eq32definition32for32AA32phase} \gamma=\theta-\alpha,\tag{2}\] where the dynamical phase is calculated by \[\label{Eq32definition32for32dynamical32phase} \alpha=-\int_{0}^{T}\langle\psi|H|\psi\rangle d\tau.\tag{3}\]
The geometric property of AA phase is embodied in the projective Hilbert space \(\mathcal{P}\) (see Appendix 7), in which the curve \(\mathcal{C}\) traversed by the system decides the value of AA phase. Quantities with such property are referred to as geometric quantities by AA.
Another geometric quantity under discussion is the time-energy uncertainty, defined as the time integral of the standard of energy: \[\label{Eq32uncertainty} s=2\int \frac{\Delta E(t)}{\hbar}dt,\tag{4}\] where \[\label{Eq32energy32uncertainty} \Delta E(t)=\left[\langle \psi | H^2 |\psi \rangle-\langle \psi | H |\psi \rangle^2\right]^{1/2}.\tag{5}\] The coefficients in Eq. (4 ) enables \(s\) to be equivalent to the distance of curve \(\mathcal{C}\) along which the state evolves in projective Hilbert space \(\mathcal{P}\), measured by the Fubini-Study metric[19]. Owing to the periodicity of the Hamiltonian in Eq. (1 ), the curve \(\mathcal{C}\) will overlap beyond a single period. As such, it is adequate to compute the time-energy uncertainty within one period.
Suppose that \(|\psi\rangle\) and \(|\psi+d\psi\rangle\) are separated by an infinitesimal distance. Then the following expression for an infinitesimal length of path \(ds\) traversed by the state vector can be derived: \[\label{Eq32ds} ds^2=4(1-|\langle \psi |\psi+d\psi \rangle^2|)\tag{6}\]
The projective Hilbert space for a two-level system can be identified as the Bloch sphere. In the present case, the Fubini-Study metric is the usual metric on the Bloch sphere with unit radius. Therefore, \(ds\) can also be obtained as: \[\label{Eq32d32theta} ds=d\phi\tag{7}\] where \(\phi\) is the angle between \(|\psi\rangle\) and \(|\psi+d\psi\rangle\). The equivalence of Eq. (6 ) and Eq. (7 ) can also be easily verified by means of numerical calculation.
We present the \(\Delta\)-dependence of the AA phase, the time-energy uncertainty and the quasienergy (see Appendix 8.1) for different bias using numerical methods in Fig. 1. With the Hamiltonian in Eq. (1 ), we first numerically solve the Schrödinger equation to obtain the evolution operator \(U(t)\): \[\label{Eq32SchrodingerU} i\partial_t U(t) = H(t) U(t).\tag{8}\] Then both the total phase \(\theta\) and the cyclic initial state \(|\psi(0)\rangle\) can be obtained from the eigenvalues and eigenvectors of \(U(t=T)\). Finally, through numerical integration, we can calculate the AA phase using Eq. (2 ) and Eq. (3 ), and the time-energy uncertainty using Eq. (4 ) and Eq. (5 ). On the other hand, the quasienergy is obtained by diagonalizing the Floquet Hamiltonian numerically, which we will discuss in detail in Sec. 5.
For the semiclassical Rabi model, we obtain AA phases \(\gamma_{\pm}\) which correspond to different cyclic initial states and satisfy \(\gamma_+ + \gamma_- = 2l \pi,l\in \mathbb{Z}\) (see Appendix 8). Nevertheless, the values of their time-energy uncertainty are equal, which means their trajectories in the projective Hilbert space are identical. Note that we have selected the interval \([0,2\pi]\) as the principal value range for the AA phase and set \(\gamma_+ + \gamma_- = 2 \pi\).




Figure 1: AA phases \(\gamma_{\pm}\), time-energy uncertainty \(s\) and quasienergy \(q\) (see Appendix 8.1) as a function of \(\Delta/\omega\) for different values of \(\epsilon/\omega\) with a fixed driving strength \(A/\omega =1\). (a) \(\epsilon/\omega=0\). (b) \(\epsilon/\omega=0.5\). (c) \(\epsilon/\omega=1\). (d) \(\epsilon/\omega=1.5\). \(\gamma_+\) and \(\gamma_-\) are AA phases corresponding to different cyclic initial states, and we set here \(\gamma_-\equiv 2\pi -\gamma_+\) according to the complementary relation. \(\gamma_{\pm}\) and \(s\) are all divided by \(\pi\), and \(q\) is divided by \(\omega\) in the figures. We have selected the interval \([0,2\pi]\) as the principal value range for the AA phase..
We show the characters of the resonance from the three quantities: AA phases \(\gamma_{\pm}\), time-energy uncertainty \(s\) and quasienergy \(q\). The resonance feature of the time-energy uncertainty is quite pronounced. At the harmonic resonance, the uncertainty reaches its local maxima, which are close to \(2\pi, 4\pi\) and \(6\pi\) for the main, second and third resonance respectively. Meanwhile, in the vicinity of the harmonic resonance corresponding to the uncertainty, both branches of the AA phase \(\gamma_{\pm}\) pass through \(\gamma=\pi\), and the number of intersections increases with the order of resonance. At the same time, the branches of quasienergy exhibits avoided crossing. Through numerical calculations, we have found that the harmonic resonance of the three quantities tend to coincide.
In the top and middle panels of Fig. 1 (a) we displays the geometric quantities for the Rabi model without bias, i.e., \(\epsilon=0\). In this case, only the main and third harmonic resonance occur at \(\Delta/\omega \approx 1\) and \(\Delta/\omega \approx 3\), respectively. Near these resonances, \(\gamma_+\) (or \(\gamma_-\)) passes through \(\pi\) once and more than once, respectively. In the bottom panel of Fig. 1 (a), quasienergy spectrum near the main and third harmonic resonance regime exhibits anti-crossing but that near \(\Delta/\omega=2\) exhibits crossing. It indicate that the corresponding relation between geometric quantities and quasienergy spectra.
When \(\epsilon /\omega >0\), there happens the second harmonic resonance contributed by the bias term at \(\Delta/\omega\approx 2\) besides the main and third harmonic resonances, as observed in Fig. 1 (b). With the increase of \(\epsilon /\omega\), all the resonance shift toward smaller values of \(\Delta/\omega\), and the main harmonic resonance eventually vanishes as \(\epsilon/\omega \ge 1\), as shown in Fig. 1 (c) and Fig. 1 (d). When \(\epsilon/\omega =1\), the third harmonic resonance also disappears, because both geometric quantities are relevant to the hidden symmetry of the asymmetric semiclassical Rabi model, which we will come back to discuss in Sec. 5.
To explore the condition of resonance and visualize its features, in Fig. 2, we demonstrate the population dynamics of the Rabi model and trajectories on the Bloch sphere for different tunneling strength \(\Delta/\omega\) and initial states. For \(A/\omega=1\) and \(\epsilon/\omega= 0.8\), the third harmonic resonance happens at \(\Delta/\omega=2.7993\), and we use \(|\psi (0)\rangle\equiv |\psi_{\rm{res}}(0)\rangle\) to denote one of its cyclic initial states, while at the near-resonance position \(\Delta/\omega=2.7\), one of its cyclic initial states is denoted as \(|\psi (0)\rangle|_{\Delta/\omega=2.7}\). Figures 2(a)-(d) show the population of the spin-up state \([1\quad0]^{T}\) in the \(\sigma_z\) basis as a function of \(t\) for different cases, calculated by \[{P_{\rm{up}}}( t)=\overline{ \frac{\sigma_z( t)+1}{2}}=\left\langle \psi(0)\left|U(t)^\dagger\frac{\sigma_z+1}{2}U(t)\right|\psi(0)\right\rangle,\] which represents the dynamics of the two-level systems. Figures 2(e)-(h) show the corresponding trajectories of \(|\psi(t)\rangle =U(t)| \psi(0)\rangle\) on the Bloch sphere, with red dots labeling the initial states and black dots labeling the final states.
In Fig. 2 (a) and Fig. 2 (e), \(\Delta/\omega=2.7993\), which is the resonance position, and the initial state is \(|\psi_{\rm{res}}(0)\rangle\). The features of resonance are manifestly shown by the rapid oscillation of \(P_{\rm{up}}\) and the trajectory on the Bloch sphere whose length is very close to \(6\pi\). As seen in Fig. 2 (e), the initial and final states converge because the initial state is the cyclic initial state for \(\Delta/\omega=2.7993\). In Fig. 2 (b) and Fig. 2 (f), \(\Delta/\omega=2.7\), which is corresponding to the nonresonance position, and the initial state is \(|\psi (0)\rangle|_{\Delta/\omega=2.7}\). In this case, resonance is not observed as \(P_{\rm{up}}\) varies slowly and the trajectory greatly contracts.
In Figs. 2 (c) and 2 (g), we show the population dynamics and the trajectory on the Bloch sphere for \(\Delta/\omega=2.7993\) initially prepared from the initial cyclic state with \(|\psi (0)\rangle|_{\Delta/\omega=2.7}\), respectively. We find that the dynamics and trajectory on the Bloch sphere in this case are similar to those results in Fig. 2 (b) and Fig. 2 (f), indicating that resonance disappears. In Fig. 2 (d) and Fig. 2 (h), if the initial state is initially prepared as \(|\psi_{\rm{res}}(0)\rangle\) and a nonresonant parameter \(\Delta/\omega=2.7\) is set, it is found that like-resonance phenomena recur which seems to those in Fig. 2 (a) and Fig. 2 (e). Note that the initial and final states diverge in Fig. 2 (g) and Fig. 2 (h) because in both cases the initial states and the parameters \(\Delta/\omega\) do not match.








Figure 2: (a-d) Time evolution of \(P_{\rm{up}}\) for different tunneling strength \(\Delta/\omega\) and initial states, with \(A/\omega=1\) and \(\epsilon/\omega= 0.8\) fixed: (a) \(\Delta/\omega=2.7993\) and the initial state is \(|\psi_{\rm{res}}(0)\rangle\); (b) \(\Delta/\omega=2.7\) and the initial state is \(|\psi (0)\rangle|_{\Delta/\omega=2.7}\); (c) \(\Delta/\omega=2.7993\) and the initial state is \(|\psi (0)\rangle|_{\Delta/\omega=2.7}\) ; (d) \(\Delta/\omega=2.7\) and the initial state is \(|\psi_{\rm{res}}(0)\rangle\) . (e-h) Trajectories on the Bloch sphere with the parameters and the initial states corresponding to (a-d), respectively. The red arrows and dots on the Bloch spheres represent the initial states, while the black dots label the final states. The arrows on the trajectories represent the directions of quantum state evolution. In (f) and (g), we zoom in on the divergence and the convergence of the initial and the final states, respectively..
From Fig. 2, we conclude that harmonic resonance mainly comes from the initial state \(|\psi (0)\rangle\). As long as one cyclic initial state corresponding to its resonance position is selected, like-Resonance phenomena can happen even if the parameter \(\Delta/\omega\) does not match that of the resonance. Meanwhile, \(|\psi (0)\rangle\) is highly sensitive to \(\Delta/\omega\) in the vicinity of harmonic resonance, as is observed from Fig. 2 (e) and Fig. 2 (g).
In this section, we apply the CHRW method to analytically calculate AA phase of the asymmetric semiclassical Rabi model. We first introduce the CHRW methodology in Sec. 3.1, and then present the analytical expression for AA phase in Sec. 3.2. Finally, we calculate the positions of the higher-order harmonic resonances by the CHRW method in Sec. 3.3.
We perform the unitary transformation with a generator \(S(t)=-i\frac{A}{2\omega}\sin(\omega t)\left(\xi \sigma_z+\zeta \sigma_x\right)\), to the Hamiltonian in Eq. (1 ). The two parameters \(\xi\) and \(\zeta\) are determined later. In the rotating frame, we obtain the evolution operator \[i\frac{d {U'(t)}}{dt} =H'(t) U'(t)\] in which \(U'(t)=e^{S(t)}U(t)\) and \(H'(t)=e^{S(t)}H(t)e^{-S(t)}-ie^{S(t)}\partial_t e^{-S(t)}\). Then, the wave function is obtained \(|\psi'(t)\rangle=e^{S(t)} |\psi\rangle\), and the transformed Hamiltonian is written \[\begin{align} H'=&-\frac{\Delta}{2}\left[\sigma_x-\frac{1-\cos \Theta}{x^2} \xi\left(\xi \sigma_x-\zeta \sigma_z\right)+\frac{\sin \Theta}{x} \xi \sigma_y\right] \\ &-\frac{\varepsilon(t)}{2}\left[\sigma_z+\frac{1-\cos \Theta}{x^2} \zeta\left(\xi \sigma_x-\zeta \sigma_z\right)-\frac{\sin \Theta}{x} \zeta \sigma_y\right] \\ &+\frac{A}{2}(\xi \sigma_z+\zeta \sigma_x)\cos (\omega t), \end{align}\] where \(\Theta=z \sin (\omega t)\), \(z=\frac{A}{\omega}x\), and \(x=\sqrt{\xi^2+\zeta^2}\). After using the identity \[\begin{align} \exp(iy \sin \alpha)=\sum_{-\infty}^{\infty}J_n(y)e^{in\alpha}, \end{align}\] where \(J_n(y)\) are the nth-order Bessel functions of the first kind, we divide the Hamiltonian into four parts according to the order of the harmonics, \[H'=H_0'+H_1'(t)+H_2'(t)+V(t),\] where \[H_0'=-\frac{\tilde{\Delta}}{2}\sigma_x-\frac{\tilde{\epsilon}}{2}\sigma_z,\] \[\begin{align} H_1^{\prime}=&-\frac{(\Delta \xi-\epsilon \zeta)}{x} J_1(z) \sin (\omega t) \sigma_y-\frac{A}{2}\left[1-\xi-\zeta^2 j_c\right]\\ & \cos (\omega t) \sigma_z +\frac{A}{2} \zeta\left[1-\xi j_c\right] \cos (\omega t) \sigma_x, \end{align}\] \[\begin{align} H_2^{\prime}=&\frac{A}{2} \frac{\zeta}{x} J_1(z) \sin (2 \omega t) \sigma_y \\ &-\frac{(\Delta \xi-\epsilon \zeta)}{x^2} J_2(z) \cos (2 \omega t)\left(\xi \sigma_x-\zeta \sigma_z\right), \end{align}\] \[\begin{align} V=&\frac{A}{2} \frac{\zeta}{x^2} J_2(z) \cos (3 \omega t)\left(\xi \sigma_x-\zeta \sigma_z\right)\\ &-\frac{[\Delta \xi-\varepsilon(t) \zeta]}{x^2} \sum_{n=2}^{\infty}\bigl\{x J_{2 n-1}(z) \sin [(2 n-1) \omega t] \sigma_y \\ &+J_{2 n}(z) \cos (2 n \omega t)\left(\xi \sigma_x-\zeta \sigma_z\right) \bigr\}. \end{align}\] The parameters \(\tilde{\Delta},\tilde{\epsilon}\) and \(j_c\) are defined as \[\tilde{\epsilon}=\epsilon+\frac{\zeta}{x^2}\left[1-J_0(z)\right](\Delta\xi-\epsilon\zeta),\] \[\tilde{\Delta}=\Delta-\frac{\xi}{x^2}\left[1-J_0(z)\right](\Delta\xi-\epsilon\zeta),\] \[j_c=\frac{1-J_0(z)-J_2(z)}{x^2},\] respectively. Note that the zero-\(\omega\) Hamiltonian \(H_0'\) consists of the renormalized tunneling term and renormalized bias one, single-\(\omega\) Hamiltonian \(H_1'\) corresponds to single-harmonic processes, and double-\(\omega\) Hamiltonian \(H_2'\) specifies second-order harmonic processes. Finally, \(V(t)\) includes all higher-order harmonic terms \(n\omega,n\geq 3\). Until now, no approximations are made, and the effects of the counterrotating terms and bias are all taken into account. Next, we keep all the zeroth and first harmonics of the transformed Hamiltonian \((n\omega, n = 0,1)\) i.e., \(H'\approx H_0'+H_1'\) and neglect the higher-order harmonic terms that involves all multi-\(\omega\) or multi-photon assisted transitions \((n\omega, n = 2,3,4,...)\). The validity of the omission of \(H_2'+V\) depends on the effects of the higher-frequency driving terms whose contribution to the dynamics is not prominent except for the ultra-strong driving-strength case[57].
The Hamiltonian \(H'_0\) is diagonalized by a unitary matrix \(D=u\sigma_z-v\sigma_x\) with \(u=\sqrt{\frac{1}{2}-\frac{\tilde{\epsilon}}{2\tilde{\Xi}}}\) and \(v=\sqrt{\frac{1}{2}+\frac{\tilde{\epsilon}}{2\tilde{\Xi}}}\) to the form \[\tilde{H}_0=\frac{\tilde{\Xi}}{2}\tau_z,\] where \(\tilde{\Xi}=\sqrt{\tilde{\Delta}^2+\tilde{\epsilon}^2}\) is the renormalized energy splitting,\(\tau_z\) is the z-component pseudo-spin operator in the energy eigenbasis.
The Hamiltonian \(H'_1\) is then changed to : \[\begin{align} \tilde{H}_1&=D^{\dagger}H_1'D\\ &= \frac{(\Delta \xi-\epsilon \zeta)}{x} J_1(z) \sin (\omega t) \tau_y \\ &+\frac{A}{2}\left[1-\xi-\zeta^2 j_c\right] \cos (\omega t) \left(\frac{\tilde{\epsilon}}{\widetilde{\Xi}} \tau_z+\frac{\tilde{\Delta}}{\widetilde{\Xi}} \tau_x\right)\\ &+\frac{A}{2} \zeta\left[1-\xi j_c\right] \cos (\omega t)\left(\frac{\tilde{\epsilon}}{\widetilde{\Xi}} \tau_x-\frac{\tilde{\Delta}}{\widetilde{\Xi}} \tau_z\right) \end{align}\] \(\tau_x\) and \(\tau_y\) are, respectively, the x-component and the y-component spin operators in the energy eigenbasis. In order to make the driving interaction term \(D^{\dagger}H_1'D\) hold the RWA-like form, we choose the two proper parameters \(\xi\) and \(\zeta\) to satisfy the following two self-consistent equations: \[\label{xizeta1} 0=\frac{A}{2}\left[\frac{\tilde{\Delta}}{\tilde{\Xi}}\left(1-\xi-\zeta^2j_c\right)+\frac{\tilde{\epsilon}}{\tilde{\Xi}}\zeta(1-\xi j_c)\right]-\frac{\Delta \xi-\epsilon\zeta}{x}J_1(z),\tag{9}\] \[\label{xizeta2} 0=\tilde{\epsilon}(1-\xi-\zeta^2j_c)-\tilde{\Delta}\zeta(1-\xi j_c).\tag{10}\] Therefore, we obtain: \[\begin{align} \tilde{H}=&D^{\dagger}(H_0'+H_1')D\\ &=\frac{\tilde{\Xi}}{2}\tau_z + \frac{\tilde{A}}{2} \left[ \tau_+ \exp (-i\omega t) + \tau_- \exp(i\omega t) \right], \end{align}\] where \(\tau_\pm=(\tau_x\pm i\tau_y)/2\) and \(\tilde{A}\) is the renormalized amplitude of the driving field resulting from the combination of the counterrotating coupling and static bias, \[\begin{align} \tilde{A}=\frac{\Delta \xi-\epsilon \zeta}{x}2J_1 \left(z\right). \end{align}\]
The solutions to the Schrödinger equation corresponding to \(\tilde{H}\) is denoted as \(|\tilde{\psi}\rangle\), and satisfy \(|\tilde{\psi}\rangle=D^{\dagger}|\psi'\rangle=D^{\dagger}e^S |\psi\rangle\). Meanwhile, the evolution operator \(\tilde{U}(t)\) corresponding to \(\tilde{H}\) is solved analytically \[\label{Eq32U95tilde} \begin{align} &\tilde{U}(t)=\\ &\left[ \begin{array}{cc} e^{-i \frac{\omega t}{2}}\left[\cos \left(\frac{\tilde{\Omega} t}{2}\right)-\frac{i \tilde{\delta}}{\tilde{\Omega}} \sin \left(\frac{\tilde{\Omega} t}{2}\right)\right] & -e^{-i \frac{\omega t}{2}} \frac{i \tilde{A}}{\tilde{\Omega}} \sin \left(\frac{\tilde{\Omega} t}{2}\right) \\ -e^{i \frac{\omega t}{2}} \frac{i \tilde{A}}{\tilde{\Omega}} \sin \left(\frac{\tilde{\Omega} t}{2}\right) & e^{i \frac{\omega t}{2}}\left[\cos \left(\frac{\tilde{\Omega} t}{2}\right)+\frac{i \tilde{\delta}}{\tilde{\Omega}} \sin \left(\frac{\tilde{\Omega} t}{2}\right)\right] \end{array}\right], \end{align}\tag{11}\] where \(\tilde{\Omega} = \sqrt{\tilde{\delta}^2 + \tilde{A}^2 }\) and \(\tilde{\delta} = \tilde{\Xi} -\omega\) are, respectively, the modulated Rabi frequency and the renormalized detuning parameter given by the CHRW method.
By the analytical expression of transformed evolution operator \(\tilde{U}(t)\), we calculate the cyclic states \(|\psi_{\pm}(t)\rangle\) by \[|\psi_{\pm}(t)\rangle=e^{-S(t)}D |\tilde{\psi}_{\pm}(t)\rangle=e^{-S(t)}D\tilde{U}(t)|\tilde{\psi}_{\pm}(0)\rangle.\] and total phases \(\theta_{\pm}\) by their eigenvalues \(e^{i\theta_{\pm}}\). The dynamical phases are obtained \[\label{Eq32alpha32formula} \begin{align} \alpha_{\pm} &=- \int_{0}^{T} \langle \psi_{\pm}|H|\psi_{\pm} \rangle d\tau\\ &=- \int_{0}^{T} \langle \tilde{\psi}_{\pm}(0)|\tilde{U}^{\dagger}D^{\dagger}e^SHe^{-S}D\tilde{U}|\tilde{\psi}_{\pm}(0) \rangle d\tau. \end{align}\tag{12}\] By \(\tilde{U}(T)\) in Eq. (11 ), \[\tilde{U}(T) = -\cos \left(\frac{\tilde{\Omega}T}{2}\right)I +i\frac{\tilde{A}}{\tilde{\Omega}}\sin \left(\frac{\tilde{\Omega}T}{2}\right)\tau_x + i\frac{\tilde{\delta}}{\tilde{\Omega}}\sin \left(\frac{\tilde{\Omega}T}{2}\right)\tau_z,\] we get its eigenvectors and phases of eigenvalues \[|\tilde{\psi}_{\pm}(0)\rangle=\sqrt{\frac{2 \tilde{\Omega}}{\tilde{\Omega} \mp \tilde{\delta}}}\left(\begin{array}{c} \frac{1}{2} \mp \frac{\tilde{\delta}}{2 \tilde{\Omega}} \\ \mp \frac{\tilde{A}}{2 \tilde{\Omega}} \end{array}\right),\] \[\theta_{\pm}=\pm \frac{\tilde{\Omega}-\omega}{2}T,\] respectively. Dynamic phases read \[\label{Eq32alpha95CHRW} \begin{align} \alpha_{\pm}&=\pm \frac{T}{4\tilde{\Omega}} \Bigg\{ \epsilon\Bigg[ \frac{\tilde{\epsilon}}{\tilde{\Xi}}\tilde{\delta}\left(1+J_0(z)\right) +\frac{2\xi \zeta}{x^2} \frac{\tilde{\Delta}}{\tilde{\Xi}} \tilde{\delta} \left(1-J_0(z)\right)\\ & - \frac{2\zeta}{x}\tilde{A} J_1(z)+\frac{\xi^2-\zeta^2}{x^2}\frac{\tilde{\epsilon}}{\tilde{\Xi}}\tilde{\delta}\left(1-J_0(z)\right) \Bigg]\\ &+A\tilde{A} \Bigg[\frac{\xi^2}{x^2} \frac{\tilde{\Delta}}{\tilde{\Xi}} +\frac{\xi \zeta}{x^2} \frac{\tilde{\epsilon}}{\tilde{\Xi}}\left(\frac{2J_1(z)}{z}-1\right) +\frac{2\zeta^2}{x^2}\frac{\tilde{\Delta}}{\tilde{\Xi}}\frac{J_1(z)}{z}\Bigg]\\ &+2 \Delta \Bigg[\frac{\zeta^2}{x^2} \frac{\tilde{\Delta}}{\tilde{\Xi}}\tilde{\delta} + \frac{\xi^2}{x^2} \frac{\tilde{\Delta}}{\tilde{\Xi}} J_0(z)\tilde{\delta} +\frac{\xi\zeta}{x^2}\frac{\tilde{\epsilon}}{\tilde{\Xi}}\tilde{\delta}(1-J_0(z)) \\ &+\frac{\xi}{x}\tilde{A}J_1(z) \Bigg] \Bigg\} \end{align}\tag{13}\] Finally we obtain the analytical expression of AA phases
\[\label{Eq32theta95CHRW} \begin{align} &\gamma_{\pm} =\theta_{\pm} -\alpha_{\pm} \\ &=\pm\frac{T}{4\tilde{\Omega}} \Bigg\{ 2\tilde{\Omega}\left(\tilde{\Omega}-\omega\right) -\epsilon\Bigg[ \frac{\tilde{\epsilon}}{\tilde{\Xi}}\tilde{\delta}\left(1+J_0(z)\right) +\frac{2\xi \zeta}{x^2} \frac{\tilde{\Delta}}{\tilde{\Xi}} \tilde{\delta} \left(1-J_0(z)\right) - \frac{2\zeta}{x}\tilde{A} J_1(z)+\frac{\xi^2-\zeta^2}{x^2}\frac{\tilde{\epsilon}}{\tilde{\Xi}}\tilde{\delta}\left(1-J_0(z)\right) \Bigg] -A\tilde{A} \\ &\Bigg[\frac{\xi^2}{x^2} \frac{\tilde{\Delta}}{\tilde{\Xi}} +\frac{\xi \zeta}{x^2} \frac{\tilde{\epsilon}}{\tilde{\Xi}}\left(\frac{2J_1(z)}{z}-1\right) + \frac{2\zeta^2}{x^2}\frac{\tilde{\Delta}}{\tilde{\Xi}}\frac{J_1(z)}{z} \Bigg]-2 \Delta \left[ \frac{\zeta^2}{x^2} \frac{\tilde{\Delta}}{\tilde{\Xi}}\tilde{\delta} + \frac{\xi^2}{x^2} \frac{\tilde{\Delta}}{\tilde{\Xi}} J_0(z)\tilde{\delta} +\frac{\xi\zeta}{x^2}\frac{\tilde{\epsilon}}{\tilde{\Xi}}\tilde{\delta}(1-J_0(z)) \;+\frac{\xi}{x}\tilde{A}J_1(z) \right] \Bigg\}. \end{align}\tag{14}\]
For the symmetric case, i.e., \(\epsilon=0\), the result of the parameter \(\zeta=0\) is self-consistently solved by Eq. (9 ) and Eq. (10 ). This simplifies Eq. (14 ) to the previous result derived in [58] albeit with slight differences in the coefficients due to the present definition of \(\tilde{A}\): \[\gamma_\pm = \pm \left[ \frac{\tilde{\Omega}-\omega}{2}-\frac{\tilde{A}}{2\tilde{\Omega}}\left(\frac{\tilde{\Delta}\tilde{\delta}}{\tilde{A}}+\frac{A}{2}+\frac{\tilde{A}}{2}\right)\right]T.\]
By the Rabi frequency of the CHRW method, we calculate the positions of the higher-order harmonic resonance[58]. The second harmonic resonance occurs when the modulated effective Rabi frequency equals the frequency of external driving field, i.e., \[\label{Eq32condition95for95res2} \tilde{\Omega} = \omega.\tag{15}\] Likewise, the condition for the third harmonic resonance is \[\label{Eq32condition95for95res3} \tilde{\Omega} = 2\omega.\tag{16}\] An analytical expression for the Rabi frequency up to the second order in the driving strength \(A\) is derived in Ref. [57]: \[\label{Eq32Rabi95freq} \tilde{\Omega}_{\rm 2nd}=\left(\omega-\Xi_0\right)^2+\frac{A^2\Delta^2}{2\Xi_0\left(\omega+\Xi_0\right)},\tag{17}\] where \(\Xi_0=\sqrt{\Delta^2+\epsilon^2}\). We use the condition \(\tilde{\Omega}_{\rm 2nd}= \omega\) to estimate the position of the second harmonic resonance.


Figure 3: (a) Position of the second harmonic resonance. (b) Position of the third harmonic resonance. \(0<\epsilon/\omega<0.5, A/\omega=1\) for both figures. The red lines are calculated by numerically exact method, the dashed dotted lines are obtained by Eq. (15 ) and Eq. (16 ) with the CHRW method, and the blue line in (a) is given by \(\tilde{\Omega}_{\rm 2nd}=\omega\) as an approximation..
It is obvious to see that, in Fig. 3, the positions of the harmonic resonance of numerically exact results agrees well with those of the CHRW method. Both the positions of the second harmonic resonance and third harmonic resonance show the tendency of decrease of the tunneling strength as the bias increases, indicating a compensating effect between \(\Delta\) and \(\epsilon\). The result obtained by using \(\tilde{\Omega}=\omega\), which is shown by the blue line in Fig. 3 (a), is in agreement with the numerically exact solution when \(\epsilon < \Delta\), despite a tiny deviation, which originates from neglecting higher order terms of \(A\) in deriving Eq. (17 ). However, the relation \(\tilde{\Omega}_{\rm 2nd}=2\omega\) could not correctly predict the third harmonic resonance because \(\tilde{\Omega}_{\rm 2nd}\) does not take into account the competition effect of bias and counter-rotating terms in higher-order resonance parameter regime, such as \(\Delta \sim 3 \omega\).
Although the CHRW method can predict the positions of the higher-order harmonic resonance of the AA phase and depict its asymptotic tendency as a function of \(\Delta/\omega\), it is necessary to combine it with perturbation theory in order to reveal the harmonic resonance features. Consider the time-dependent Hamiltonian \(\tilde{H}'\) split into a unperturbed Hamiltonian \(\tilde{H}\) and a perturbation with a dimensionless parameter \(\lambda\) representing the order of perturbation, \[\tilde{H}'=\tilde{H}+\lambda\tilde{H}_{pt},\] and its corresponding time evolution operator \[\tilde{U}'=\tilde{U}+\lambda\tilde{U}_{pt}+...,\] where \(i\frac{d\tilde{U}}{dt} =\tilde{H}\tilde{U}\). Since \(\tilde{U}'(t)\) satisfies the time-dependent Schrödinger equation \(i\frac{d\tilde{U}'}{dt} =\tilde{H}'\tilde{U}'\), we can solve \(\tilde{U}_{pt}\) to the first order \[i\frac{d\tilde{U}_{pt}}{dt} = \tilde{H}\tilde{U}_{pt}+\tilde{H}_{pt}\tilde{U}.\] Then we analytically solve \[\label{Eq32U95pt} \tilde{U}_{pt}(t) = -i\tilde{U}\int_0^t \tilde{U}^{-1}\tilde{H}_{pt}\tilde{U} d\tau.\tag{18}\] In this work we choose the parameter, \(A/\omega <2\), and therefore the influences of second harmonic terms in \(H_2'(t)\) are dominant compared with the other higher harmonic terms in \(V(t)\). In the following, by considering the second harmonic terms in \(H_2'\), we calculate evolution operator under perturbation and then derive modified cyclic initial states and more accurate analytical results for the AA phase. In the energy eigenbasis representation of diagonalization, \(H_2'\) becomes \[\tilde{H}_2=D^{\dagger}H_{2}'D=\tilde{H}_{2y}+\tilde{H}_{2x}+\tilde{H}_{2z},\] where \[\tilde{H}_{2y}=-\frac{A\zeta}{2x}J_1(z)\sin (2\omega t)\tau_y,\] \[\tilde{H}_{2x}=\frac{\epsilon \zeta-\Delta \xi}{x^2}J_2(z) q \cos (2\omega t) \tau_x,\] \[\tilde{H}_{2z}=\frac{\epsilon \zeta-\Delta \xi}{x^2}J_2(z) p \cos (2\omega t) \tau_z,\] and \(p = (v^2-u^2)\zeta-2\xi uv\) and \(q =(v^2-u^2)\xi +2uv\zeta\). Here we take into consideration the effects of the terms with \(\tau_x,\tau_y\) and \(\tau_z\) separately. First, we explore the effects of \(\tilde{H}_{2y}\), calculating the integral: \[\begin{align} &\int_{0}^{T} \tilde{U}^{-1}\tilde{H}_{2y}\tilde{U} d\tau = \frac{2A\zeta J_1(z) \sin \left(\frac{\tilde{\Omega}}{2}T\right)}{x\tilde{\Omega}^2 (9\omega^4 -10\omega^2\tilde{\Omega}^2 +\tilde{\Omega}^4)} \\ &\Bigg[ \left( 3\tilde{\delta}\tilde{\Omega}\omega^3+ 2\tilde{\delta}^2\tilde{\Omega}\omega^2-\tilde{\delta}\tilde{\Omega}^3\omega \right) \cos\left(\frac{\tilde{\Omega}}{2}T\right) \tau_x \\ &+ \left( \tilde{\Omega}^4\omega-2\tilde{\Omega}^2\tilde{\delta}\omega^2-3\tilde{\Omega}^2\omega^3 \right)\sin \left(\frac{\tilde{\Omega}}{2}T\right) \tau_y \\ &+ \tilde{A} \left(\tilde{\Omega}^3 \omega -3\tilde{\Omega}\omega^3-2\tilde{\Omega}\tilde{\delta}\omega^2 \right)\cos\left(\frac{\tilde{\Omega}}{2}T\right) \tau_z \Bigg]. \end{align}\] Substituting Eq. (11 ) into Eq. (18 ), we obtain its first order correction to \(\tilde{U}'\): \[\tilde{U}_{2y}(T)=i \frac{2A\zeta \omega J_1(z) \sin (\frac{ \tilde{\Omega}}{2}T) (3\omega+2\tilde{\delta}\omega-\tilde{\Omega}^2) }{x\tilde{\Omega}(9\omega^4 -10\omega^2\tilde{\Omega}^2 +\tilde{\Omega}^4)}(\tilde{\delta}\tau_x -\tilde{A}\tau_z).\] The whole time evolution operator \(\tilde{U}'_{y}(T) = \tilde{U}(T) + \tilde{U}_{2y}(T)\) yields \[\begin{align} \tilde{U}'_{y}(T)&=-\cos \left(\frac{\tilde{\Omega}T}{2}\right)I +\frac{i}{\tilde{\Omega}}\sin \left(\frac{\tilde{\Omega}T}{2}\right)\\ &\left[ (\tilde{A}+{k_y}\tilde{\delta}) \tau_x+(\tilde{\delta}-{k_y}\tilde{A} )\tau_z \right], \end{align}\] where \[\label{Eq32ky} {k_y}=\frac{2A\zeta \omega J_1(z) (3\omega+2\tilde{\delta}\omega-\tilde{\Omega}^2) }{x(9\omega^4 -10\omega^2\tilde{\Omega}^2 +\tilde{\Omega}^4)}.\tag{19}\] We immediately obtain one of its eigenvectors and its corresponding eigenvalue \[\label{Eq32tilde32p32y39} |\tilde{\psi}_{+y}'(0) \rangle= \frac{1}{{L_y}}\left[ \begin{array}{c} \tilde{A}+{k_y}\tilde{\delta} \\ {k_y}\tilde{A}-\tilde{\delta}-\sqrt{1+{k_y}^2}\tilde{\Omega} \end{array} \right],\tag{20}\] \[\label{Eq32eigenvalue95pt} \lambda_+ = -\cos \left(\frac{\tilde{\Omega}T}{2}\right) -i\sqrt{1+{k_y}^2} \sin \left(\frac{\tilde{\Omega}T}{2}\right),\tag{21}\] where the normalization factor \[\label{Eq32Ly} {L_y}=\left[2(1+{k_y}^2)\tilde{\Omega}^2-2(\tilde{\delta}+{k_y}\tilde{A})\sqrt{1+{k_y}^2}\tilde{\Omega}\right]^{1/2}.\tag{22}\] Because the eigenvalue \(|\lambda_+|\simeq 1\), we can regard \[\label{Eq32total95phase32pt} \theta_+' = {\rm{arg}}(\lambda_+)=\arctan \left[\sqrt{1+k_y^2} \tan \left(\frac{\tilde{\Omega}T}{2}\right)\right]-\frac{\omega T}{2}\tag{23}\] as the modified total phase under perturbation, which includes the effects of counterrotating terms. Also, the perturbed Rabi frequency in the vicinity of harmonic resonance can be derived accordingly (see Appendix 9). In fact, \(\theta_+'\) is very close to \(\theta_+\) given by the CHRW method, as the eigenvalues of \(\tilde{U}'_{y}(T)\) are very close to those of \(\tilde{U}(T)\). Thus we can use the latter to approximate \(\theta_+'\), and the less \(|{k_y}|\) is, the better the approximation is.
Owing to the great contribution \(\tilde{U}_{2y}(T)\) to \(\tilde{U}'_{y}(T)\), both dynamical phase and AA phase change sharply near higher-order harmonic resonance. Using Eq. (2 ) and Eq. (12 ), we obtain the AA phase corresponding to Eq. (20 ) \[\gamma_+=\theta_+ +\int_{0}^{T} \langle \tilde{\psi}_{+y}'(0)|\tilde{U}^{\dagger}D^{\dagger}e^SHe^{-S}D\tilde{U}|\tilde{\psi}_{+y}'(0)\rangle d\tau\] Here we use \(\tilde{U}\) instead of \(\tilde{U}'\) in the integrand, which causes neglectable difference in the calculation. In light of the complementary relation \(\gamma_+ +\gamma_- = 2l\pi,l\in \mathbb{Z}\), we obtain both AA phases approximately
\[\label{Eq32gamma95pt952y} \begin{align} \gamma_\pm\approx &\pm\frac{\tilde{\Omega}-\omega}{2}T\mp\frac{\pi}{2{L_y}^2} \Bigg\{ (\tilde{\delta}+\sqrt{{k_y}^2+1}\tilde{\Omega}-{k_y}\tilde{A}) \Big\{\frac{A^2 \tilde{\delta}}{2\omega^3}\left[\epsilon(\nu^2-\mu^2)-2\Delta \mu \nu\right] -\frac{\tilde{\delta}}{128\omega}(3z^4-64z^2+512)\\ &\left(u^2\epsilon-2uv\Delta-v^2\epsilon\right)-\frac{2A\tilde{A}}{\omega^2}\left[(u\epsilon-v\Delta)\nu-(u\Delta+v\epsilon)\mu\right]+\frac{A\tilde{A}}{\omega}\left[\frac{2}{z}J_1(z)+1\right]\left[{k_y}(u^2-v^2)+2uv\right]\Big\}\\ &+\frac{4{k_y}}{x^2}J_2(z)\left(u\nu-v\mu\right)\left[{k_y}\tilde{A}(\tilde{\Omega}+\tilde{\delta})-\left(1+\sqrt{{k_y}^2+1}\right)\tilde{\Omega}\tilde{\delta}-\sqrt{{k_y}^2+1}\tilde{\Omega}^2-\tilde{\delta}^2\right]\Bigg\}, \end{align}\tag{24}\]
where \(\mu = \xi u-\zeta v, \nu = \xi v+\zeta u\). The details of derivation are given in Appendix 9.
Likewise, when we take account of the influences of \(\tilde{H}_{2x}\) and \(\tilde{H}_{2z}\), we repeat the processes of derivation above and obtain \[\begin{align} \tilde{U}_{2x}(T)=\frac{ik_x\sin \left(\frac{\tilde{\Omega}}{2}T\right)}{\tilde{\Omega}}\left(\tilde{\delta}\tau_x -\tilde{A}\tau_z\right), \\ \tilde{U}_{2z}(T)=\frac{ik_z\sin \left(\frac{\tilde{\Omega}}{2}T\right)}{\tilde{\Omega}}\left(\tilde{\delta}\tau_x -\tilde{A}\tau_z\right). \end{align}\] where \[\label{Eq32kx} k_x=\frac{2(\Delta\xi -\epsilon \zeta)J_2(z)}{x^2} \frac{\left[(v^2-u^2)\xi +2uv\zeta\right](3\omega^3+5\omega^2\tilde{\delta}+\omega\tilde{\Omega}^2-\tilde{\delta}\tilde{\Omega}^2)}{(9\omega^4-10\omega^2\tilde{\Omega}^2+\tilde{\Omega}^4)},\tag{25}\] \[\label{Eq32kz} k_z=-\frac{2(\Delta\xi -\epsilon \zeta)J_2(z)}{x^2}\frac{\tilde{A}\left[2\xi uv+\zeta (u^2-v^2)\right]}{\tilde{\Omega}^2-4\omega^2}.\tag{26}\]
It is noticeable that both \(\tilde{U}_{2x}(T)\) and \(\tilde{U}_{2z}(T)\) possess a similar mathematical structure to \(\tilde{U}_{2y}(T)\) expect their distinct coefficients. If combining them together, we obtain the first order correction to the time evolution operator \(\tilde{U}\). Therefore, it turns out that the same forms of cyclic initial state and dynamical phase as those in Eq. (20 ) and Eq. (24 ) remains unchanged. Substituting \({k_y}\) in Eq. (22 ) and Eq. (24 ) with \(k={k_y}+{k_x}+{k_z}\) results in the perturbed result. Moreover, they are more accurate results for the AA phase in comparison with the unperturbed and numerical results which is shown in Fig. 4. Meanwhile, we can see clearly from Eq. (19 ),Eq. (25 ) and Eq. (26 ) that \(\tilde{H}_{2y}\) and \(\tilde{H}_{2x}\) contributes to the higher-order harmonic resonance for \(\tilde{\Omega}\approx \omega\) and \(\tilde{\Omega}\approx 3\omega\), while \(\tilde{H}_{2z}\) contributes to the higher-order harmonic resonance for \(\tilde{\Omega}\approx 2\omega\).




Figure 4: (a-b) AA phase \(\gamma_+/\pi\) as a function of \(\Delta/\omega\) for different values of \(\epsilon/\omega\) with \(A/\omega=1\) fixed: (a) \(\epsilon/\omega=0.5\); (b) \(\epsilon/\omega=1.5\). (c-d) Modulated effective Rabi frequency \(\tilde{\Omega}/\omega\) as a function of \(\Delta/\omega\) for different values of \(\epsilon/\omega\) with \(A/\omega=1\) fixed: (c) \(\epsilon/\omega=0.5\); (d) \(\epsilon/\omega=1.5\). The results of perturbation theory based on the CHRW method, denoted as CHRW+PT, agree well with the numerically exact results, while those obtained by CHRW show asymptotic tendency. In (c) and (d), we include magnified views to highlight the sudden changes at the harmonic resonance at \(\tilde{\Omega}/\omega\approx2\)..
To evaluate the effectiveness of perturbation theory based on the CHRW method, we depict the AA phase \(\gamma_+\) and the modulated effective Rabi frequency \(\tilde{\Omega}\) as a function of \(\Delta/\omega\) in Fig. 4, obtained by the numerical method, the CHRW method, and perturbation theory based on the CHRW method. In Fig. 4 (a) and Fig. 4 (b), we remove the limitation of AA phase to \([0,2\pi]\) to better exhibit the harmonic resonance peaks of AA phase. It is obvious to see that the results obtained by the CHRW method exhibit an asymptotic tendency for both the AA phase and the modulated effective Rabi frequency. In contrast, the combination of perturbation theory and the CHRW method (denoted as CHRW+PT) accurately depicts the harmonic resonance peaks of \(\gamma_+\) and sudden changes of \(\tilde{\Omega}\) in the second and the third harmonic resonance regimes, defined as the full width at half maximum (FWHM) of AA phase[58]. In Fig. 4 (a) and Fig. 4 (b), we observe not only the third harmonic resonance driven by the external field, but also the second harmonic resonance arising from the combination of bias and driving field. It is known that the odd higher-order harmonic resonance happens in the symmetric Rabi model. However, the asymmetry in the asymmetric Rabi model induces even higher-order harmonic resonance which is a significant character in the asymmetric Rabi model, compared with the results of the symmetric Rabi model[58]. In Fig. 4 (a) where \(\epsilon/\omega=0.5\), the second harmonic resonance occurs at \(\Delta/\omega \approx 1.75\), and the third harmonic resonance occurs at \(\Delta/\omega \approx 2.865\). As the bias increases to \(\epsilon/\omega=1.5\) in Fig. 4 (b), the positions of the second and the third harmonic resonances shift to \(\Delta/\omega \approx 1.25\) and \(\Delta/\omega \approx 2.515\) respectively, which reflects the compensating effect between \(\Delta\) and \(\epsilon\) as demonstrated in Sec. 3.3. Moreover, as shown in Fig. 4 (c) and Fig. 4 (d), the second and the third harmonic resonances happen when \(\tilde{\Omega}/\omega \approx 1\) and \(\tilde{\Omega}/\omega \approx 2\) respectively, in agreement with the prediction of the CHRW method.
Now we come back to discuss the hidden symmetry of the asymmetric semiclassical Rabi model that appears as \(\epsilon=m\omega, m\in \mathbb{Z}\). After careful calculation, we confirm that as \(\epsilon=m\omega, m\in \mathbb{Z}\), the harmonic resonance of geometric quantities is absent for the numerical results at the third harmonic resonance \(\Delta_{\rm{res}}\). Furthermore, we find that in this case, the geometric quantities are actually determined by the arbitrary choices of cyclic initial states. In this section, we first illustrate this discovery using Floquet theory (about which we give a short overview in Appendix 8.1). We also explain this phenomenon in terms of hidden symmetry by comparing the asymmetric semiclassical Rabi model with the asymmetric quantum Rabi model.
First of all, one can numerically verify the fact that when \(\epsilon=\omega\) and \(\Delta= \Delta_{\rm{res}}\), the quasi-energies satisfy: \[q_- - q_+ = m\omega, m\in \mathbb{Z},\] which is equivalent to \[\label{Eq32equivalent32total32phase} \theta_+ - \theta_-= 2m\pi, m\in\mathbb{Z}.\tag{27}\] According to Floquet theory, both quasi-energies are physically identical, which means the quasi-energies are degenerate, and any superposition of the cyclic initial states \(|\Psi\rangle = c_1 |\psi_{-}(0)\rangle + c_2 |\psi_{+}(0)\rangle\) is also a cyclic initial state with the equivalent quasienergy[59], i.e., \[U(T)|\Psi\rangle=e^{-iq_+T}|\Psi\rangle.\] Because \(|\psi_{-}(0)\rangle\) and \(|\psi_{+}(0)\rangle\) are linear independent, \(|\Psi\rangle\) can be any state vector in the Hilbert space, and therefore the invariant space of \(U(T)\) is the whole space. Thus, we can naturally infer that \(U(T)=e^{-iq_+T}I\). Furthermore, combining Eq. (27 ) with the complementary relation \(\theta_+ + \theta_-= 2l\pi,l\in \mathbb{Z}\), we can set \(\theta_+=\theta_-=0\), without loss of generality, and in this case, \(U(T)=I\).
Since the cyclic initial state can be, in principle, any vector in the Hilbert space, and geometric quantities which depend on it cannot be uniquely defined, including AA phase and time-energy uncertainty as our focuses in this work. In this sense, the absent resonance obtained numerically is only one of the infinity solutions. It can be verified that the degeneracy of quasi-energies as \(\Delta= \Delta_{\rm{res}}\) happens not only when \(\epsilon=\omega\), but when \(\epsilon=m\omega, m\in \mathbb{Z}\) as well.
This intriguing phenomenon stems from the hidden symmetry of the asymmetric Rabi model, which has been reported by a number of researches for the asymmetric quantum Rabi model[54], [64]–[66]. The Hamiltonian for the asymmetric quantum Rabi model can be written as: \[\label{Eq32Hq} H_q=\omega a^\dagger a- \frac{g}{4} \sigma_z (a^\dagger +a)-\frac{\Delta}{2}\sigma_x-\frac{\epsilon}{2}\sigma_z,\tag{28}\] where \(\omega\) is the frequency of a quantized resonator field with annihilation and creation operators \(a\) and \(a^\dagger\), and \(g\) is the coupling strength between the qubit and the field. For this Hamiltonian, the hidden symmetry appears as long as \(\epsilon=m\omega, m\in \mathbb{Z}\), and conserved quantity which commutes with the Hamiltonian, and corresponding symmetric unitary operator can be found in the case [64]–[66].
It is necessary to mention that the coefficients and signs chosen in Eq. (28 ) are corresponding to those in Eq. (1 ), and the properties of the asymmetric quantum Rabi model are not reliant on those particular choices. The quantum field in Eq. (28 ) can be reduced to the classical field in Eq. (1 ) with \(g\left|\langle a \rangle\right|=A\) fixed and \(g\rightarrow 0, \left|\langle a \rangle\right| \rightarrow \infty\), where \(\left|\langle a \rangle\right|\) is expectation value of the annihilation operator \(a\) [67].
According to Floquet theory, the Floquet Hamiltonian for the asymmetric semiclassical Rabi model, is \[\label{Eq32HF} H_F = -\frac{\Delta}{2} \sigma_x - \frac{\epsilon+A\cos(\omega t)}{2} \sigma_z -i\partial_t.\tag{29}\] A comparison of Eq. (29 ) and Eq. (28 ) in the matrix form demonstrates the similarity between the asymmetric semiclassical Rabi model and the asymmetric quantum Rabi model. For the asymmetric semiclassical Rabi model, we introduce a set of basis: \[\label{Eq32basis32for32asymmetric32semiclassical32Rabi32model} \left| \uparrow \rm{or} \downarrow, n\right\rangle = \left| \uparrow \rm{or} \downarrow \right\rangle e^{-in\omega t},\tag{30}\] where \(\left| \uparrow \right\rangle\) and \(\left| \downarrow \right\rangle\) are the eigenstates of \(\sigma_z\). In the extended Hilbert space \(\mathcal{H} \otimes \mathcal{T}\), those vectors become \(\left| \uparrow \rm{or} \downarrow, n\rangle \right\rangle\) (see Appendix 8.1 for the details). In this basis, the matrix of \(H_F\) reads:
\[\label{Eq32HF32matrix} H_F=\left[ \begin{array}{c|cccccc} \ddots & |\uparrow,-n\rangle\rangle & |\downarrow,-n\rangle\rangle & |\uparrow,-(n+1)\rangle\rangle & |\downarrow,-(n+1)\rangle\rangle & |\uparrow,-(n+2)\rangle\rangle & |\downarrow,-(n+2)\rangle\rangle \\ \hline|\uparrow,-n\rangle\rangle & n \omega-\frac{\epsilon}{2} & -\frac{\Delta}{2} & -\frac{A}{4} & 0 & 0 & 0 \\ |\downarrow,-n\rangle\rangle & -\frac{\Delta}{2} & n \omega+\frac{\epsilon}{2} & 0 & \frac{A}{4} & 0 & 0 \\ |\uparrow,-(n+1)\rangle\rangle & -\frac{A}{4} & 0 & (n+1) \omega-\frac{\epsilon}{2} & -\frac{\Delta}{2} & -\frac{A}{4} & 0 \\ |\downarrow,-(n+1)\rangle\rangle & 0 & \frac{A}{4} & -\frac{\Delta}{2} & (n+1) \omega+\frac{\epsilon}{2} & 0 & \frac{A}{4} \\ |\uparrow,-(n+2)\rangle\rangle & 0 & 0 & -\frac{A}{4} & 0 & (n+2) \omega-\frac{\epsilon}{2} & -\frac{\Delta}{2} \\ |\downarrow,-(n+2)\rangle\rangle & 0 & 0 & 0 & \frac{A}{4} & -\frac{\Delta}{2} & (n+2) \omega+\frac{\epsilon}{2} \end{array} \right].\tag{31}\]
For the asymmetric quantum Rabi model, we introduce another similar set of basis \[\label{Eq32basis32for32asymmetric32quantum32Rabi32model} \left| \uparrow \rm{or} \downarrow, n \right\rangle,\tag{32}\] where \(\left| \uparrow \right\rangle\) and \(\left| \downarrow \right\rangle\) are the same as the symbols in Eq. (30 ) and \(| n \rangle\) is the fork state for the periodic driving field. In this basis, \(H_q\) takes the form \[\label{Eq32Hq32matrix} H_q=\left[ \begin{array}{c|cccccc} \ddots & |\uparrow,n\rangle & |\downarrow,n\rangle & |\uparrow,n+1\rangle & |\downarrow,n+1\rangle & |\uparrow,n+2\rangle & |\downarrow,n+2\rangle \\ \hline|\uparrow,n\rangle & n \omega-\frac{\epsilon}{2} & -\frac{\Delta}{2} & -\frac{g}{4}\sqrt{n+1} & 0 & 0 & 0 \\ |\downarrow,n\rangle & -\frac{\Delta}{2} & n \omega+\frac{\epsilon}{2} & 0 & \frac{g}{4}\sqrt{n+1} & 0 & 0 \\ |\uparrow,n+1\rangle & -\frac{g}{4}\sqrt{n+1} & 0 & (n+1) \omega-\frac{\epsilon}{2} & -\frac{\Delta}{2} & -\frac{g}{4}\sqrt{n+2} & 0 \\ |\downarrow,n+1\rangle & 0 & \frac{g}{4}\sqrt{n+1} & -\frac{\Delta}{2} & (n+1) \omega+\frac{\epsilon}{2} & 0 & \frac{g}{4}\sqrt{n+2} \\ |\uparrow,n+2\rangle & 0 & 0 & -\frac{g}{4}\sqrt{n+2} & 0 & (n+2) \omega-\frac{\epsilon}{2} & -\frac{\Delta}{2} \\ |\downarrow,n+2\rangle & 0 & 0 & 0 & \frac{g}{4}\sqrt{n+2} & -\frac{\Delta}{2} & (n+2) \omega+\frac{\epsilon}{2} \end{array} \right].\tag{33}\]
The structural similarity between \(H_F\) and \(H_q\) is apparent in the matrix representation: the exactly equal tridiagonal elements for both Hamiltonian (note that \(n\in \mathbb{Z}^{+}\) for both matrices and \(-n\) is used in Eq. (31 )). Meanwhile, the off-tridiagonal elements of \(H_q\) depend on \(n\) and those of \(H_F\) depend on \(A\), which is related to the average number of photons in the semiclassical limit[67], [68]. However, in \(H_F\), \(n\) runs from negative infinity to positive infinity, while in \(H_q\), \(n\) cannot be negative because it labels the number of photons [39].




Figure 5: Comparison of (quasi-)energy spectrum as a function of \(\Delta/\omega\) between the asymmetric semiclassical Rabi model and the asymmetric quantum Rabi model under different parameters. (a) and (c) Quasi-energy spectrum of the asymmetric semiclassical Rabi model with \(A/\omega=1, \epsilon/\omega=1,0.8\). (b) and (d) Eigen-energy spectrum of the asymmetric quantum Rabi model with \(g/\omega=1, \epsilon/\omega=1,0.8\). For \(\epsilon/\omega=1\), the quasi-energies or eigen-energies are degenerate at certain points, while no degeneracy occurs for \(\epsilon/\omega=0.8\)..
With the matrices of both Hamiltonian, we calculate the convergence values of the eigenvalues to obtain accurate energy spectra by incrementally increasing their dimensions. Figure 5 (a) and Figure 5 (c) display the quasienergy of the asymmetric semiclassical Rabi model against \(\Delta/\omega\) with \(A/\omega=1\) for \(\epsilon/\omega=1\) and \(\epsilon/\omega=0.8\), respectively, while Figure 5 (b) and Figure 5 (d) the eigenenergy of the asymmetric quantum Rabi model with \(g/\omega=1\), for \(\epsilon/\omega=1\) and \(\epsilon/\omega=0.8\). From these figures, we can see that for \(\epsilon/\omega=1\), the quasi-energies are degenerate at \(\Delta_{\rm{res}}/\omega\approx 2.74\), and degeneracies can also be observed for the eigenenergy spectrum in the vicinity of \(\Delta/\omega\approx 2.5\). However, no level crossing occurs when \(\epsilon/\omega\ne 1\). Though not shown here, we also verified the general existence of level crossings for \(\epsilon/\omega= m,m\in \mathbb{Z}\), and anti-crossings for \(\epsilon/\omega \ne m\), and we found that these characters were not influenced by the value of \(A/\omega\) or \(g/\omega\), which mainly shifted the spectra along the \(\Delta/\omega\) axis. This fact provides strong evidence for the hidden symmetry of both the asymmetric semiclassical Rabi model and the asymmetric quantum Rabi model as structural properties of the Hamiltonian, which probably have impacts in all coupling regimes of light-matter interaction or driven two-level system. In consequence, the hidden symmetry contributes to the non-uniqueness of the geometric quantities, AA phase and time-energy uncertainty.
In this work, we investigate the harmonic resonance of AA phase and time-energy uncertainty as geometric quantities and hidden symmetry in the asymmetric Rabi model by numerical method. At the same time, we apply perturbation theory based on the CHRW method to analytically give accurate time evolution operator and then calculate accurate geometric phase which is in good agreement with numerically exact results. The asymmetry in the Rabi model have important effects on the geometric phase. In comparison with the existence of only odd harmonic resonances at the absence of bias, various resonance phenomena take place at the presence of bias, including the emergence of the even higher-order harmonic resonances induced by the asymmetry, and the disappearance of odd higher-order resonance caused by the hidden symmetry present when \(\epsilon=m\omega,m\in \mathbb{Z}\). Besides, with the increase of the bias, all the harmonic resonance shift toward smaller values of \(\Delta/\omega\) owing to the compensating effect between \(\epsilon\) and \(\Delta\), confirmed by both the numerical method and the CHRW method using the condition \(\tilde{\Omega}=m\omega, m\in \mathbb{Z}\).
We illustrate the features of harmonic resonance encapsulated by the repeated intersection of AA phase at \(\gamma=\pi\), attainment of a local maximum in time-energy uncertainty, and rapid oscillations in the population of the spin-up state, \(\rm{P_{up}}\). The harmonic resonance was mainly steered by the initial state \(|\psi (0)\rangle\), which is highly sensitive to \(\Delta/\omega\) in the vicinity of harmonic resonance. The results obtained by the CHRW method, in which the influences of counterrotating terms and bias were all taken into account in the rotating frame, exhibit asymptotic tendencies for both AA phase and modulated effective Rabi frequency. With the combination of perturbation theory, we further effectively depicting the harmonic resonance peaks of \(\gamma_+\) and sudden changes of \(\tilde{\Omega}\) in the second and the third harmonic resonance regimes.
Upon investigation of the absence of harmonic resonance in the results of numerical calculation for \(\epsilon=m\omega,m \in \mathbb{Z},\Delta= \Delta_{\rm{res}}\), we shed light on the dependence of geometric quantities on the choice of cyclic initial states due to the degeneracy of quasienergy of the asymmetric semiclassical Rabi model. Using Floquet theory, we found the structural similarity between \(H_F\) and \(H_q\) in the matrix representation and spectra, which reveals the hidden symmetry in the asymmetric semiclassical Rabi model as an analogy of that in the asymmetric quantum Rabi model reported already in [54], [64]–[66]. Our results suggest that the hidden symmetry reflects the intrinsic property of asymmetric Rabi model regardless of the amplitude of driving field or the coupling strength between the qubit and oscillating field, which plays a potential role in all coupling regimes of light-matter interaction or driven two-level system. Meanwhile, since the the time-energy uncertainty equals the length of path traversed by the state on the Bloch sphere, it can be measured in superconducting qubit systems via quantum process tomography.
The work is supported by National Natural Science Foundation of China (Grants No. 11774226 and No. 61927822).
In quantum mechanics, states are described by vectors in the Hilbert space \(\mathcal{H}\). However, any two vectors \(\psi,\phi \in \mathcal{H}\), satisfying \[\psi=c\phi,c\in \mathbb{C},\] are physically equivalent \((\psi\sim\phi)\). Therefore, the projective Hilbert space \(\mathcal{P}\), defined as: \[\mathcal{P} :=\mathcal{H}/\sim\] is the proper space for representing quantum states that can be physically distinguished.
For a time-dependent Hamiltonian which satisfies \(H(t+T)=H(t)\), the solution to the Schrödinger equation \[\label{Schrödinger32equation} i\partial_t |\psi(t)\rangle = H(t) |\psi(t)\rangle\tag{34}\] takes the form: \[\label{Eq32floquet} |\psi_\alpha(t)\rangle= |u_\alpha(t)\rangle e^{-iq_\alpha t},\tag{35}\] where \(q_\alpha\) is called quasienergy, and \(|u_\alpha(t+T)\rangle = |u_\alpha(t)\rangle\). This is often referred to as Floquet theorem. Since \(|\psi_\alpha(0)\rangle\) merely acquires a phase factor after \(t=T\), i.e., the total phase \(\theta = -q_\alpha T\), it is called a cyclic initial state. Inserting Eq. (35 ) into Eq. (34 ), the original problem can be transformed into an eigenvalue problem with Floquet Hamiltonian \(H_F\): \[H_F \left|u_\alpha(t)\right\rangle =\left[H(t)-i\partial_t\right]|u_\alpha(t)\rangle = q_\alpha |u_\alpha(t)\rangle.\] Note that the subscript \(\alpha\) is to label the order of the eigenvectors or eigenvalues, but there are actually infinite set of solutions, because \(\left|u_{\alpha}(t),n\right\rangle = e^{-in\omega t} \left|u_{\alpha}(t),0\right\rangle \equiv e^{-in\omega t} \left|u_{\alpha}(t)\right\rangle\) is physically identical to Eq. (35 ) only with the shifted quasienergy \(q_{\alpha,n}\equiv q_\alpha-n\omega\). Thus, it is sufficient to consider the set of quasi-energies within an interval of width \(\omega\) and corresponding total phases within an interval of width \(2\pi\).
More importantly, \(H_F\) can be written as a time-independent infinite matrix through the expansion of the original Hilbert space \(\mathcal{H}\) to \(\mathcal{H} \otimes \mathcal{T}\)[39], where \(\mathcal{T}\) contains all T-periodic functions. A particular set of orthonormal complete basis in \(\mathcal{T}\) is written as \(\{|n\rangle\}\) with the plane wave function \[\langle t|n \rangle = e^{-in\omega t}.\]
Considering a T-periodic state vector \(\left|u_\alpha(t),n\right\rangle \in \mathcal{H}\), we can expand it in a Fourier series, plane wave terms are derived: \[\left|u_{\alpha}(t), n\right\rangle=e^{-i n \omega t}\left|u_\alpha(t)\right\rangle=\sum_l e^{-i l \omega t}|u_\alpha^{(n-l)}\rangle,\] where \(|u_\alpha^{(k)}\rangle\) are time-independent Fourier coefficients. Now in \(\mathcal{H} \otimes \mathcal{T}\), the corresponding state vector is defined as \[\left| u_\alpha, n\rangle \right\rangle \equiv \sum_l\left|u_\alpha^{(n-l)}\right\rangle \otimes \left|l\right\rangle,\] and finally the inner product is defined as \[\langle \langle u_{\alpha},n|u_{\beta},m\rangle\rangle \equiv \frac{1}{T}\int_{0}^{T} dt \left\langle u_{\alpha}(t),n|u_{\beta}(t),m\right\rangle.\]
To derive the complementary relation of AA phases, we start from the Schrödinger equation in Eq. (8 ). For this purpose, here we denote the Hamiltonian as \(H(t)=[h_1,h_2;h_2,-h_1]\), where \(h_1=-\left(A \cos \omega t+\epsilon\right)/2, h_2=-\Delta/2\), and \(U(t)=[u_1,u_2;u_3,u_4]\). Then we can list the following ordinary differential equations (ODE’s): \[\begin{align} i\partial_t u_1=h_1u_1+h_2u_3, \tag{36}\\ i\partial_t u_2=h_1u_2+h_2u_4, \tag{37}\\ i\partial_t u_3=h_2u_1-h_1u_3, \tag{38}\\ i\partial_t u_4=h_2u_2-h_1u_4, \tag{39} \end{align}\] with boundary conditions: \[\begin{align} u_1(0)=u_4(0)=1, \nonumber\\ u_2(0)=u_3(0)=0. \nonumber \end{align}\] Note that, if we substitute the variables \(u_1\) and \(u_2\) in Eq. (38 ) and Eq. (39 ) with \(u_4^*\) and \(-u_3^*\) respectively and take the conjugate of both equations, we will have the same equations as Eq. (36 ) and Eq. (37 ), with the boundary conditions still satisfied. Due to the uniqueness theorem for ODE’s, the relations \[\label{Eq32elementU} u_1= u_4^*, \quad u_2= -u_3^*\tag{40}\] must be kept in the solution for Eq. (8 ). Since the evolution operator is unitary, we can get \[\begin{align} U(t)U(t)^{\dagger}&= \left[\begin{array}{cc} u_1 & u_2 \\ -u_2^* & u_1^* \end{array}\right] \left[\begin{array}{cc} u_1^* & -u_2 \\ u_2^* & u_1 \end{array}\right]\\ &=\left[\begin{array}{cc} \left|u_1\right|^2 +\left|u_2\right|^2 & 0 \\ 0 & \left|u_1\right|^2 +\left|u_2\right|^2 \end{array}\right]\\ &=I, \end{align}\] which implies that \[\label{Eq32detU} det\left[U(t)\right] =\left|u_1\right|^2 +\left|u_2\right|^2=1, \forall t\ge 0.\tag{41}\] Because \(e^{i\theta_+}\) and \(e^{i\theta-}\) are the two eigenvalues of \(U(T)\), \(e^{i\theta_+} e^{i\theta-} = e^{i \theta_+ + \theta_-}=det\left[U(T)\right]=1\). Therefore, we get the relation between total phases \[\label{Eq32complementary32relation32for32total32phase} \theta_+ + \theta_-= 2l\pi,l\in \mathbb{Z}.\tag{42}\]
Moreover, with Eq. (40 ) and Eq. (41 ) it is easy to prove through some calculations that the two eigenvectors for \(U(T)\), i.e., the cyclic initial states \(|\psi_{\pm}(0)\rangle\) are orthogonal to each other, and thus \[\langle \psi_{-}(0)|U^{\dagger}HU|\psi_{-}(0)\rangle=-\langle \psi_{+}(0)|U^{\dagger}HU|\psi_{+}(0)\rangle,\] which leads to \[\label{Eq32complementary32relation32for32dynamical32phase} \alpha_-=-\alpha_+.\tag{43}\] Finally we obtain from Eq. (2 ),Eq. (3 ),Eq. (42 ) and Eq. (43 ) that \(\gamma_+ + \gamma_- = 2l\pi,l\in \mathbb{Z}\).
Near the third harmonic resonance, \(\tilde{\Omega}\approx 2\omega\). Thus, the total phase given in Eq. (23 ) can be simplified as \[\theta_+'\approx \frac{\sqrt{1+k^2}\left(\tilde{\Omega}-2\omega\right)T-\omega T}{2}.\] By \(\theta_{+}'= \frac{\tilde{\Omega}'}{2}T+\frac{2n+1}{2}\omega T, n\in \mathbb{Z}\), we get the modified Rabi frequency \[\tilde{\Omega}'\approx \sqrt{1+k^2}\left(\tilde{\Omega}-2\omega\right)+2\omega.\] To calculate the dynamical phases, our target is to solve the integral in Eq. (12 ), since the total phases given by the CHRW method are reserved due to their tiny variation. In the first step, we write the middle part of the integrand
\[\label{Eq32DeHeD} \begin{align} &D^{\dagger}e^SHe^{-S}D=-\frac{\tau_y}{2} \Bigg\{ \frac{ \sin\Theta}{x}\left[ \nu(u\varepsilon-v\Delta)-\mu(u\Delta +v\varepsilon) \right] \Bigg\} -\frac{\tau_x}{2}\Bigg\{ \cos^2\left(\frac{\Theta}{2}\right)\left( v^2\Delta-2uv\varepsilon-u^2\Delta \right)+\frac{1}{x^2}\sin^2\left(\frac{\Theta}{2}\right)\\ &\left( \Delta\left( \mu^2-\nu^2 \right) -2\varepsilon \mu\nu \right)\Bigg\}-\frac{\tau_z}{2}\Bigg\{ \cos^2\left(\frac{\Theta}{2}\right)(u^2 \varepsilon-2uv\Delta-v^2\varepsilon)+\frac{1}{x^2}\sin^2\left(\frac{\Theta}{2}\right)\left[ \varepsilon\left( \mu^2-\nu^2 \right)+2\Delta \mu \nu\right] \Bigg\}. \end{align}\tag{44}\]
Next, we calculate the state vector at time \(t\), using the evolution operator \(\tilde{U}\) without perturbation
\[\label{Eq32U95tilde32p95tilde39} \tilde{U} |\tilde{\psi}_{+}'(0) \rangle = \frac{1}{{L}} \left(\begin{array}{c} e^{-i\omega t/2}\left[ (\tilde{A}+{k}\tilde{\delta})\cos \left(\frac{\tilde{\Omega}t}{2}\right)+i(\tilde{A}\sqrt{1+{k}^2}-{k}\tilde{\Omega})\sin \left(\frac{\tilde{\Omega}t}{2}\right) \right] \\ e^{i\omega t/2}\left[ ({k}\tilde{A}-\tilde{\delta}-\sqrt{1+{k}^2}\tilde{\Omega})\cos \left(\frac{\tilde{\Omega}t}{2}\right)-i(\tilde{\Omega}+\sqrt{1+{k}^2}\tilde{\delta})\sin \left(\frac{\tilde{\Omega}t}{2}\right) \right] \end{array}\right),\tag{45}\] where the initial state under perturbation \(|\tilde{+}' \rangle\) takes the form of Eq. (20 ) with all components \(\tau_x,\tau_y\) and \(\tau_z\) taken into consideration.
To calculate the integral in Eq. (12 ) analytically, we also use \(\omega\) to approximate \(\tilde{\Omega}\), which leads to negligible effects. Last but not least, in the process, Taylor expansion is adopted \[\cos\left(\frac{\Theta}{2}\right)\approx 1-\frac{1}{2}\left(\frac{\Theta}{2}\right)^2,\sin \left(\frac{\Theta}{2}\right)\approx \frac{\Theta}{2}.\] Combining Eq. (44 ) and Eq. (45 ) with the approximations above, we then calculate the dynamical phase given by perturbation theory \[\label{Eq32alpha95p21} \begin{align} &\alpha_+' \approx- \int_{0}^{T} \langle \tilde{\psi}_{+}'(0)|\tilde{U}^{\dagger}D^{\dagger}e^SHe^{-S}D\tilde{U}|\tilde{\psi}_{+}'(0) \rangle d\tau \approx \frac{\pi}{2{L}^2} \Bigg\{ (\tilde{\delta}+\sqrt{{k}^2+1}\tilde{\Omega}-{k}\tilde{A}) \bigg\{\frac{A^2 \tilde{\delta}}{2\omega^3}\left[\epsilon(\nu^2-\mu^2)-2\Delta \mu \nu\right]\\ &-\frac{\tilde{\delta}}{128\omega}(3z^4-64z^2+512)(u^2\epsilon-2uv\Delta-v^2\epsilon)-\frac{2A\tilde{A}}{\omega^2}\left[(u\epsilon-v\Delta)\nu-(u\Delta+v\epsilon)\mu\right]+\frac{A\tilde{A}}{\omega}\left[\frac{2}{z}J_1(z)+1\right][k(u^2-v^2)+2uv]\bigg\} \\ &+\frac{4{k}}{x^2}J_2(z)\left(u\nu-v\mu\right)\left[{k}\tilde{A}(\tilde{\Omega}+\tilde{\delta})-(1+\sqrt{{k}^2+1})\tilde{\Omega}\tilde{\delta}-\sqrt{{k}^2+1}\tilde{\Omega}^2-\tilde{\delta}^2\right]\Bigg\}. \end{align}\tag{46}\]