Accelerated adiabatic passage in cavity magnomechanics

The TQD approach was proposed in the first decade of this century, depending on the full knowledge about the instantaneous eigenstructure of the original Hamiltonian. Conventionally, if the original time-dependent Hamiltonian $H(t)$ could be formally expressed in the spectral representation as $H(t)=\sum_{n}E_{n}(t)|n(t)\rangle\langle n(t)|$, then assisted by an ancillary Hamiltonian, or called the counterdiabatic (CD) Hamiltonian Berry (2009); Guéry-Odelin et al. (2019)

$\displaystyle H_{\rm CD}(t)=i\sum_{n}\left[1-|n(t)\rangle\langle n(t)|\right]|\dot{n}(t)\rangle\langle n(t)|,$

(10)

the system could keep track of the instantaneous eigenstates of $H(t)$ in a much faster speed. This approach is also believed to highly robust against the control-parameter variations Guéry-Odelin et al. (2019). One can understand that it applies usually to the discrete systems, since it is explicitly represented by the eigenstates and their time-derivatives.

In this work, we derive the CD term using operators in the continuous-variable systems. Using the Bogolyubov transformation Tikochinsky (1978), the Hamiltonian in Eq. (7) is rewritten as

$$H(t)=\omega_{A}A^{\dagger}A+\omega_{B}B^{\dagger}B,$$

(11)

where $\omega_{A,B}=\pm\sqrt{\Delta^{2}+4g^{2}}/2$ and

		$\displaystyle A\equiv\cos\theta m+\sin\theta b,$		(12)
		$\displaystyle B\equiv\sin\theta m-\cos\theta b$

with $\tan(2\theta)=2g(t)/\Delta(t)$. To simplify the formation of the TQD protocol, the coupling strength $g(t)$ in this protocol is set to be real, i.e., $g(t)=g^{*}(t)$.

In the subspace with a fixed and arbitrary excitation number $N$, the CD Hamiltonian for the system Hamiltonian in Eq. (11) can be written as $H_{\rm CD}=i\sum_{n=0}^{N}|\dot{\epsilon}_{n}\rangle\langle\epsilon_{n}|$, where the orthonormal eigenstates read

$$|\epsilon_{N-n}\rangle=\frac{1}{\sqrt{(N-n)!n!}}(A^{\dagger})^{N-n}(B^{\dagger})^{n}|0\rangle,$$

(13)

with $|0\rangle$ the vacuum state for both modes. Due to the fact that

		$\displaystyle\dot{A}=-\dot{\theta}\sin\theta m+\dot{\theta}\cos\theta b=-\dot{\theta}B,$		(14)
		$\displaystyle\dot{B}=\dot{\theta}\cos\theta m+\dot{\theta}\sin\theta b=\dot{\theta}A,$

we have

	$\displaystyle|\dot{\epsilon}_{N-n}\rangle$	$\displaystyle=-\dot{\theta}\sqrt{(n+1)(N-n)}|\epsilon_{N-n-1}\rangle$		(15)
		$\displaystyle+\dot{\theta}\sqrt{n(N-n+1)}|\epsilon_{N-n+1}\rangle.$

Then the CD Hamiltonian by Eq. (10) becomes

	$\displaystyle H_{\rm CD}$	$\displaystyle=i\sum_{n=0}^{N}|\dot{\epsilon}_{n}\rangle\langle\epsilon_{n}|$		(16)
		$\displaystyle=-\dot{\theta}\sum_{n=0}^{N-1}\sqrt{(N-n)(n+1)}|\epsilon_{N-n-1}\rangle\langle\epsilon_{N-n}|$
		$\displaystyle+\dot{\theta}\sum_{n=1}^{N}\sqrt{n(N-n+1)}|\epsilon_{N-n+1}\rangle\langle\epsilon_{N-n}|.$

According to the definition in Eq. (13), the first term in Eq. (16) expands

		$\displaystyle\sum_{n=0}^{N-1}\sqrt{(N-n)(n+1)}|\epsilon_{N-n-1}\rangle\langle\epsilon_{N-n}|$		(17)
	$\displaystyle=$	$\displaystyle\sum_{n=0}^{N-1}\sqrt{(N-n)(n+1)}\frac{(A^{\dagger})^{N-n-1}(B^{\dagger})^{n+1}A^{N-n}B^{n}}{(N-n-1)!n!\sqrt{(N-n)(n+1)}}$
	$\displaystyle=$	$\displaystyle\sum_{n=0}^{N-1}\frac{(A^{\dagger}A)^{N-n-1}(B^{\dagger}B)^{n}B^{\dagger}A}{(N-n-1)!n!},$
	$\displaystyle=$	$\displaystyle\frac{(A^{\dagger}A+B^{\dagger}B)^{N-1}}{(N-1)!}B^{\dagger}A=IB^{\dagger}A=B^{\dagger}A.$

where $I$ is the identify operator in the subspace with $N-1$ excitations. In the last line, we have applied the binomial theorem $(x+y)^{N}=\sum^{N}_{n=0}C^{n}_{N}x^{n}y^{N-n},C_{N}^{n}=N!/[n!(N-n)!]$. Similarly, the second term in Eq. (16) turns out to be

$$\sum_{n=1}^{N}\sqrt{n(N-n+1)}|\epsilon_{N-n+1}\rangle\langle\epsilon_{N-n}|=A^{\dagger}B.$$

(18)

Note $N$ is arbitrary, then in the whole Hilbert space, the CD Hamiltonian could be expressed as

$$H_{\rm CD}=i\dot{\theta}(A^{\dagger}B-B^{\dagger}A)=i\dot{\theta}(b^{\dagger}m-m^{\dagger}b).$$

(19)

And regarding the original Hamiltonian (7), then the total Hamiltonian for transitionless quantum driving reads

		$\displaystyle H_{\rm tot}=H(t)+H_{\rm CD}$		(20)
		$\displaystyle=\frac{\Delta(t)}{2}(m^{\dagger}m-b^{\dagger}b)+[g(t)-i\dot{\theta}]m^{\dagger}b+[g(t)+i\dot{\theta})]mb^{\dagger},$

where the time-dependence and boundary-condition of

$$\dot{\theta}=\frac{\dot{g}\Delta-\dot{\Delta}g}{\Delta^{2}+4g^{2}}$$

(21)

determines the speed of the accelerated adiabatic passage.

By virtue of the definitions in Eqs. (12) and (13), the adiabatic path from $|k\rangle_{m}|0\rangle_{b}$ to $|0\rangle_{m}|k\rangle_{b}$ under the total Hamiltonian (20) is constructed by

$$|k(t)\rangle=|\epsilon_{k}(t)\rangle=\frac{1}{\sqrt{k!}}(A^{\dagger})^{k}|0\rangle=|\epsilon_{k}\rangle$$

(22)

under the boundary condition $\theta(t=0)=0$ and $\theta(t=T)=\pi/2$. And the system wavefunction can be written as

$$|\psi_{k}(t)\rangle=e^{-i\chi_{k}(t)}|k(t)\rangle,$$

(23)

where the quantum phase is

	$\displaystyle\chi_{k}(t)$	$\displaystyle=\int^{t}_{0}dt^{\prime}E_{k}(t^{\prime})-i\int^{t}_{0}dt^{\prime}\langle\epsilon_{k}(t^{\prime})|\partial_{t^{\prime}}\epsilon_{k}(t^{\prime})\rangle$		(24)
		$\displaystyle=\int^{t}_{0}dt^{\prime}E_{k}(t^{\prime})=k\chi(t).$

Note $E_{k}(t)$ is the instantaneous eigenvalue of $|\epsilon_{k}\rangle$, $\langle\epsilon_{k}(t^{\prime})|\partial_{t^{\prime}}\epsilon_{k}(t^{\prime})\rangle=0$, and $\chi(t)\equiv\int^{t}_{0}dt^{\prime}\omega_{A}(t^{\prime})$. Alternatively, by setting $\theta(0)=\pi/2$ and $\theta(T)=0$, the adiabatic path followed by the system can be constructed by $|k(t)\rangle=|\epsilon_{0}(t)\rangle=\frac{1}{\sqrt{k!}}(B^{\dagger})^{k}|0\rangle$. To be self-consistent, we stick to the boundary condition $\theta(0)=0$ and $\theta(T)=\pi/2$ in this work.

In general situations, a superposed state $|\psi(0)\rangle=\sum_{k}C_{k}|k0\rangle$ with normalized coefficients $C_{k}$’s will adiabatically evolve to

$$|\psi(T)\rangle=\sum_{k}C_{k}e^{-ik\chi(T)}|0k\rangle,$$

(25)

at the desired moment $T$ and $C_{k}$ is invariant with time. The quantum phase $\chi_{k}(T)$ for the Fock state with $k$ excitations is proportional to $k$. By the preceding analyse, the phase difference can be periodically cancelled by the bare Hamiltonian of the mechanical mode. The state transfer thus has been indeed completed by Eq. (25).

Now we can verify the TQD approach in the state transfer from mode-$m$ to mode-$b$ by presenting the practical dynamics of the target-state population $P$ given by Eq. (9). We choose the time-dependence of the effective frequency $\Delta$ and the driving-enhanced coupling strength $g$ to be in a sinusoid shape, and by Eq. (21),

$$\Delta=2\Omega\cos(2\theta),\quad g=\Omega\sin(2\theta),\quad\theta=\frac{\pi}{2}\frac{t}{T},$$

(26)

where $\Omega$ is the coupling strength determined by the desired transfer time $T$ as $\Omega T=\pi$. The shape functions of $\Delta$, $g$, and $\dot{\theta}$ are plotted in Fig. 2(a). With various initial states of mode-$m$, the blue solid lines and the red dashed lines in Figs. 2(b), (c), and (d) describe the dynamics of the target-state population under $H_{\rm tot}$ in Eq. (20) with the counterdiabatic term and that under $H(t)$ in Eq. (7) without the counterdiabatic term, respectively.

One can find that practically the accelerated state-transfer could be perfectly completed by the TQD approach for either Fock-state and superposed states. The latter of the hybrid mode is prepared as an even cat-state $(|\zeta\rangle+|-\zeta\rangle)/\sqrt{2+2e^{-2\zeta^{2}}}$, where $|\zeta\rangle$ is the Glauber coherent state. For the Fock-state transfer $|1\rangle_{m}|0\rangle_{b}\rightarrow|0\rangle_{m}|1\rangle_{b}$ in Fig. 2(b), $P$ approaches $0.97$ by the original Hamiltonian $H(t)$. In Figs. 2(c) and (d) for the cat state with $\zeta=1$ and $4$, respectively, it is found that the TQD approach manifests its power for larger cat-state by achieving unit transfer population. In contrast, under the original Hamiltonian $H(t)$, $P$ approaches respectively $0.98$ and $0.65$ for $\zeta=1$ and $4$.

Another main-stream accelerated adiabatic-passage is the invariant-based inverse engineering  Guéry-Odelin et al. (2019); Lewis and Riesenfeld (1969), where the parametrical adiabatic-path of the system is designed through a Hermitian operator $I(t)$ termed the Levis-Riesenfeld invariant. For an arbitrary original Hamiltonian $H(t)$, the invariant satisfies

$$\frac{\partial I(t)}{\partial t}=-i[H(t),I(t)].$$

(27)

In the framework of invariant-based inverse engineering, the wavefunction of a time-dependent Schrödinger equation $i\partial_{t}|\psi(t)\rangle=H(t)|\psi(t)\rangle$ can be expressed as $|\psi(t)\rangle=\sum_{n}C_{n}e^{-i\kappa_{n}(t)}|\epsilon_{n}(t)\rangle$, where $C_{n}$ is a time-independent amplitude, $|\epsilon_{n}(t)\rangle$ is the eigenstates of the invariant $I(t)$, and $\kappa_{n}$ is the Lewis-Riesenfeld phase defined by

$$\dot{\kappa}_{n}(t)=\langle\epsilon_{n}(t)|[-i\partial_{t}+H(t)]|\epsilon_{n}(t)\rangle.$$

(28)

To ensure the desired state transfer rather than the full-time adiabatic passage, $I(t)$ and $H(t)$ have to share the same eigenstates at both initial and final moments.

To carry out the derivation of a general LR invariant for our two-coupled-resonator system with arbitrary target state, we rewrite the system Hamiltonian $H(t)$ (7) into

	$\displaystyle H(t)$	$\displaystyle=\frac{\Delta(t)}{2}m^{\dagger}m+[g_{R}(t)-ig_{I}(t)]m^{\dagger}b$		(29)
		$\displaystyle+[g_{R}(t)+ig_{I}(t)]mb^{\dagger}-\frac{\Delta(t)}{2}b^{\dagger}b.$

where $g_{R}(t)$ and $g_{I}(t)$ represent the real and imaginary parts of the complex coupling strength $g(t)$.

Its corresponding Lewis-Riesenfeld invariant can be formulated by

	$\displaystyle I(t)$	$\displaystyle=\cos\beta(m^{\dagger}m-b^{\dagger}b)+\sin\beta(e^{-i\alpha}m^{\dagger}b+e^{i\alpha}mb^{\dagger}),$		(30)
		$\displaystyle=A^{\dagger}A-B^{\dagger}B,$

where

	$\displaystyle A=\cos\left(\frac{\beta}{2}\right)e^{\frac{i\alpha}{2}}m+\sin\left(\frac{\beta}{2}\right)e^{-\frac{i\alpha}{2}}b,$		(31)
	$\displaystyle B=\sin\left(\frac{\beta}{2}\right)e^{\frac{i\alpha}{2}}m-\cos\left(\frac{\beta}{2}\right)e^{-\frac{i\alpha}{2}}b,$

are the normalized annihilation operators for $I(t)$ and both $\beta\equiv\beta(t)$ and $\alpha\equiv\alpha(t)$ are time-dependent functions to be determined. Substituting Eq. (30) into Eq. (27), we have

	$\displaystyle\dot{\beta}$	$\displaystyle=2g_{I}\cos\alpha-2g_{R}\sin\alpha,$		(32)
	$\displaystyle\dot{\alpha}$	$\displaystyle=\Delta-\cot\beta(2g_{R}\cos\alpha+2g_{I}\sin\alpha).$

With the annihilation and creation operators of $I(t)$, the original Hamiltonian (29) can be rewritten as

$$H=\omega(A^{\dagger}A-B^{\dagger}B)+g_{AB}A^{\dagger}B+g^{*}_{AB}B^{\dagger}A,$$

(33)

where

	$\displaystyle\omega$	$\displaystyle\equiv\frac{\Delta\cos\beta}{2}+g_{R}\sin\beta\cos\alpha+g_{I}\sin\beta\sin\alpha,$		(34)
	$\displaystyle g_{AB}$	$\displaystyle\equiv\frac{\Delta\sin\beta}{2}-g_{R}(\cos\beta\cos\alpha+i\sin\alpha)$
		$\displaystyle-g_{I}(\cos\beta\sin\alpha-i\cos\alpha).$

And then in the subspace with a fixed excitation number $N$, the general solution of the Schrödinger equation can be expressed as a superposition of the eigenstates of the invariant,

$$|\psi_{N}(t)\rangle=\sum^{N}_{n=0}p_{n}|\epsilon_{N-n}\rangle e^{-i\kappa_{N-n}(t)},$$

(35)

where $p_{n}$ is a normalized coefficient, $|\epsilon_{N-n}\rangle$ is the normalized eigentstates of the LR invariant taking the same form as in Eq. (13), and $\kappa_{N-n}(t)$ is the quantum phase defined in Eq. (28).

We consider the state transfer from $|N\rangle_{m}|0\rangle_{b}$ to $|0\rangle_{m}|N\rangle_{b}$. With the definitions in Eqs. (31) and (13), one can construct a particular solution via

$$|\psi_{N}(t)\rangle=|\epsilon_{N}\rangle e^{-i\kappa_{N}(t)}=\frac{e^{-i\kappa_{N}(t)}}{\sqrt{N!}}(A^{\dagger})^{N}|0\rangle.$$

(36)

under the boundary conditions $\beta(0)=0$ and $\beta(T)=\pi$. As for the Lewis-Riesenfeld phase, we have

	$\displaystyle\dot{\kappa}_{N}$	$\displaystyle=-i\langle\epsilon_{N}|\partial_{t}|\epsilon_{N}\rangle+\langle\epsilon_{N}|H|\epsilon_{N}\rangle$		(37)
		$\displaystyle=-i\left(-iN\frac{\dot{\alpha}\cos\beta}{2}\right)+N\frac{\Delta\cos\beta}{2}$
		$\displaystyle+N(g_{R}\sin\beta\cos\alpha+g_{I}\sin\beta\cos\alpha)$
		$\displaystyle=N\frac{g_{R}\cos\alpha+g_{I}\sin\alpha}{\sin\beta}=N\dot{\kappa},$

with $\dot{\kappa}\equiv(g_{R}\cos\alpha+g_{I}\sin\alpha)/\sin\beta$, in which we have applied the time-derivative of the operators $A$ and $B$

	$\displaystyle\frac{\partial}{\partial t}A^{\dagger}$	$\displaystyle=-\frac{\dot{\beta}}{2}B^{\dagger}-\frac{i\dot{\alpha}}{2}\left(\cos\beta A^{\dagger}+\sin\beta B^{\dagger}\right),$		(38)
	$\displaystyle\frac{\partial}{\partial t}B^{\dagger}$	$\displaystyle=\frac{\dot{\beta}}{2}B^{\dagger}-\frac{i\dot{\alpha}}{2}\left(\sin\beta A^{\dagger}-\cos\beta B^{\dagger}\right),$

and Eqs. (32) and (33).

Given the time-dependent parameters $\beta(t)$, $\alpha(t)$, and $\kappa(t)$ in Eqs. (32) and (37), the coupling strengths and the effective frequency of the hybrid mode that are directly relevant in quantum control can be expressed by

	$\displaystyle g_{R}$	$\displaystyle=\dot{\kappa}\cos\alpha\sin\beta-\frac{\dot{\beta}}{2}\sin\alpha,$		(39)
	$\displaystyle g_{I}$	$\displaystyle=\dot{\kappa}\sin\alpha\sin\beta+\frac{\dot{\beta}}{2}\cos\alpha,$
	$\displaystyle\Delta$	$\displaystyle=\dot{\alpha}+2\dot{\kappa}\cos\beta.$

By virtue of their excitation-number-independence, the transfer of an arbitrary superposed state from mode-$m$ to mode-$b$ can be achieved under the boundary conditions $\beta(0)=0$ and $\beta(T)=\pi$. Thus in a general situation, an initial state

$$|\psi(0)\rangle=\sum_{k}C_{k}|k0\rangle=\sum_{k}C_{k}|\psi_{k}(0)\rangle,$$

(40)

where $C_{k}$ is the time-independent normalized coefficients, will evolve to

$$|\psi(T)\rangle=\sum_{k}C_{k}e^{-ik(\kappa+\alpha)}|0k\rangle,$$

(41)

at the final time $T$, according to Eqs. (36) and (31).

The state transfer assisted by the Lewis-Riesenfeld invariant can be verified in Fig. 3 by the state population $P$ of the phonon mode-$b$ given by Eq. (9). With the selected control functions of $\beta$, $\alpha$, and $\kappa$, one can directly find the time-dependence of the real and imaginary parts of the coupling strength $g_{R}$ and $g_{I}$ and the frequency of the hybrid mode $\Delta$ through Eq. (39). They are plotted in Fig. 3(a). The initial states in Fig. 3(b) are various cat states with $\zeta=1,2,4$. It is found that a perfect transfer can always be achieved via the LR-invariant-based inverse engineering.

Before working on the two accelerated adiabatic passages, we first consider a straightforward $\pi$-pulse for the state transfer. Then in Eq. (7), it is found that $\Delta(t)=0$ and $\int^{T}_{0}dtg(t)=\pi/2$. Accordingly, Eq. (42) becomes $H_{\rm exp}=(1+\gamma)H_{g}$, which can be diagonalized with the new operators $A=(m+b)/\sqrt{2}$ and $B=(m-b)/\sqrt{2}$, similarly to the transformation in Eqs. (11) and (12). The special initial state $|\Psi(0)\rangle=|N\rangle_{m}|0\rangle_{b}=|N0\rangle$ could then expand by the eigenstates in Eq. (13)

$$|N0\rangle=\frac{1}{\sqrt{2^{N}}}\sum^{N}_{n=0}\sqrt{C^{n}_{N}}|\epsilon_{N-n}\rangle.$$

(46)

The eigenvalue of $|\epsilon_{N-n}\rangle$ is now $(N-2n)(1+\gamma)g(t)$ with respect to $H_{\rm exp}$. So that we have

$\displaystyle|\Psi(T)\rangle=\frac{1}{\sqrt{2^{N}}}\sum^{N}_{n=0}\sqrt{C^{n}_{N}}e^{-i(N-2n)(1+\gamma)\frac{\pi}{2}}|\epsilon_{N-n}\rangle.$

(47)

by the Schrödinger equation in Eq. (44). According to Eqs. (11) and (12), the target state is

$$|0N\rangle=\frac{1}{\sqrt{2^{N}}}\sum^{N}_{n=0}(-1)^{n}\sqrt{C^{n}_{N}}|\epsilon_{N-n}\rangle.$$

(48)

Then the state-transfer fidelity measured by the target-state population $P$ in mode-$b$ is

$\displaystyle P=|\langle 0N|\Psi(T)\rangle|^{2}=\cos^{2N}\left(\frac{\pi}{2}\gamma\right).$

(49)

by virtue of Eq. (9). And by Eq. (45), the systematic-error sensitivity to coupling strength for the initial Fock state $|N0\rangle$ is

$$q_{g}=\frac{N\pi^{2}}{4}.$$

(50)

The preceding derivation can be straightforwardly extended to arbitrary superposed state $|\psi(0)\rangle=\sum_{k}C_{k}|k0\rangle$ by virtue of its independence of the excitation number. In general situations, we have

		$\displaystyle P=\sum_{C_{k}\neq 0}|C_{k}|^{2}\cos^{2k}\left(\frac{\pi}{2}\gamma\right),$		(51)
		$\displaystyle q_{g}=\frac{\pi^{2}}{4}\sum_{C_{k}\neq 0}k|C_{k}|^{2}=\frac{\pi^{2}}{4}\bar{n}_{m}.$

It is interesting to find that the systematic-error sensitivity in the $\pi$-pulse protocol is proportional to the average excitation number $\bar{n}_{m}$ of the initial state.

Now we analysis the stability for state transfer and error sensitivity of the TQD protocol provided in Sec. III.1. Note the unperturbed Hamiltonian $H(t)$ in the total Hamiltonian in Eq. (20) is replaced with $H_{\rm exp}$ in Eq. (42). In Figs. 4(a) and (c) $\eta=0$, and in Figs. 4(b) and (d) $\gamma=0$. In Fig. 4, the TQD protocols are performed by the time-dependence of the effective frequency $\Delta$ and the driving-enhanced coupling strength $g(t)$ (assumed to be real for TQD), that are determined by the shape functions of $\theta=\pi/2(t/T)$ in Eq. (26) or $\theta=\pi/2(t/T)^{2}$. One can observe that the protocol stability is not sensitive to the choice of the target states. It is found that the impact of the coupling-strength fluctuation of the interaction Hamiltonian $H_{g}$ is asymmetrical to the parameter $\gamma$ in the negative and positive axis. With the same magnitude, the decrement in the state population $P$ induced by a positive $\gamma$ is clearly smaller than that by a negative $\gamma$. In particular, $P=0.99$ for $\gamma/\Omega=0.2$ and $P=0.96$ for $\gamma/\Omega=-0.2$. In contrast, the state-transfer population is roughly symmetrical to the energy fluctuation $\eta$ of the free Hamiltonian $H_{\Delta}$. Another difference between Figs. 4(a), (c) and Figs. 4(b), (d) manifests in the error sensitivity to the shape of the control parameter $\theta(t)$. For a nonvanishing $\gamma$ ($\eta$), the protocol is more robust with the linear function $\theta=\pi/2(t/T)$ [the quadratic function $\theta=\pi/2(t/T)^{2}$] than that with the quadratic function $\theta=\pi/2(t/T)^{2}$ [the linear function $\theta=\pi/2(t/T)$].

In Fig. 5, we switch simultaneously on the systematic error in both interaction Hamiltonian and free Hamiltonian and fix the target state as a cat state with $\zeta=1$. Much to our anticipation, a nearly unit state-transfer is found in a remarkable regime with $\gamma\approx\eta$. It could be readily understood that when $\gamma\approx\eta$, the experimental Hamiltonian $H_{\rm exp}$ in Eq. (42) is approximated by $(1+\gamma)H(t)$, equivalent to a rescaling over the whole original Hamiltonian that renders the same counter-diabatic Hamiltonian $H_{CD}$ in Eq. (19).

From both Figs. 4 and 5, with even $20\%$ fluctuations in system parameters, the state-transfer fidelity could be maintained above $0.95$. One can generally find that the TQD approach is robust against the systematic errors.

In this subsection, we construct an optimal protocol against the systematic errors by using the Lewis-Riesenfeld invariant in Eq. (30). Now we employ the ideal Hamiltonian $H(t)$ in Eq. (29) with a complex coupling strength $g(t)$.

From the Schrödinger equation in Eq. (44), one can obtain the time-evolved state at the final time $T$ up to the second-order of $O(\gamma^{2})$ and $O{(\eta^{2})}$:

		$\displaystyle|\Psi(T)\rangle=|\psi(T)\rangle-i\gamma\int^{T}_{0}dtU_{0}(T,t)H_{g}|\psi(t)\rangle$		(52)
		$\displaystyle-\gamma^{2}\int^{T}_{0}dt\int^{t}_{0}dt^{\prime}U_{0}(T,t)H_{g}(t)U_{0}(t,t^{\prime})H_{g}(t^{\prime})|\psi(t^{\prime})\rangle$
		$\displaystyle-i\eta\int^{T}_{0}dtU_{0}(T,t)H_{\Delta}|\psi(t)\rangle$
		$\displaystyle-\eta^{2}\int^{T}_{0}dt\int^{t}_{0}dt^{\prime}U_{0}(T,t)H_{\Delta}(t)U_{0}(t,t^{\prime})H_{\Delta}(t^{\prime})|\psi(t^{\prime})\rangle$
		$\displaystyle+\cdots,$

where $|\psi(t)\rangle$ is the unperturbed solution determined by the ideal evolution operator $U_{0}(t,0)$. Under the assumption that the initial state is a Fock state $|N\rangle_{m}|0\rangle_{b}$ and by virtue of Eqs. (35) and (36), the unperturbed solution reads

$$|\psi(t)\rangle=|\epsilon_{N}(t)\rangle e^{-iN\kappa(t)},$$

(53)

where $\kappa$ is defined in the last line of Eq. (37), and we have

$$U_{0}(s,t)=\sum^{N}_{n=0}e^{-i(N-2n)[\kappa(s)-\kappa(t)]}|\epsilon_{N-n}(s)\rangle\langle\epsilon_{N-n}(t)|.$$

(54)

Substituting Eqs. (53) and (54) to Eq. (52), the final state population is

	$\displaystyle P$	$\displaystyle=|\langle 0N|\Psi(T)\rangle|^{2}=|\langle\psi(T)|\Psi(T)\rangle|^{2}$		(55)
		$\displaystyle\approx 1-\gamma^{2}\sum^{N}_{n=1}\left|\int^{T}_{0}dte^{-2ni\kappa(t)}\langle\epsilon_{N-n}(t)|H_{g}|\epsilon_{N}(t)\rangle\right|^{2}$
		$\displaystyle-\eta^{2}\sum^{N}_{n=1}\left|\int^{T}_{0}dte^{-2ni\kappa(t)}\langle\epsilon_{N-n}(t)|H_{\Delta}|\epsilon_{N}(t)\rangle\right|^{2}.$

By virtue of

	$\displaystyle\langle\epsilon_{N-1}(t)|H_{g}|\epsilon_{N}(t)\rangle$	$\displaystyle=-\sqrt{N}\dot{\kappa}\sin\beta\cos\beta+\frac{i\sqrt{N}\dot{\beta}}{2},$		(56)
	$\displaystyle\langle\epsilon_{N-n}(t)|H_{g}|\epsilon_{N}(t)\rangle$	$\displaystyle=0,\quad n\neq 1$
	$\displaystyle\langle\epsilon_{N-1}(t)|H_{\Delta}|\epsilon_{N}(t)\rangle$	$\displaystyle=\frac{\sqrt{N}\dot{\alpha}}{2}\sin\beta+\sqrt{N}\dot{\kappa}\sin\beta\cos\beta,$
	$\displaystyle\langle\epsilon_{N-n}(t)|H_{\Delta}|\epsilon_{N}(t)\rangle$	$\displaystyle=0,\quad n\neq 1$

and the boundary condition $\beta(0)=0$ and $\beta(0)=\pi$, the systematic error sensitivities for $|\psi(0)\rangle=|N0\rangle$ are found to be

	$\displaystyle q_{g}$	$\displaystyle=N\left|\int^{T}_{0}dt\dot{\beta}\sin^{2}\beta e^{-2i\kappa}\right|^{2},$		(57)
	$\displaystyle q_{\Delta}$	$\displaystyle=N\left|\int^{T}_{0}dt\sin\beta\left(\frac{\dot{\alpha}}{2}+\dot{\kappa}\cos\beta\right)e^{-2i\kappa}\right|^{2},$

according to their definitions in Eq. (45). Specially when $\kappa$ is constant, we restore the result in the $\pi$-pulse case $q_{g}=N\pi^{2}/4$, which is independent of $\beta(t)$.