††thanks<EMAIL_ADDRESS><EMAIL_ADDRESS>
# Robust population inversion in three-level systems by composite pulses

Cheng Zhang  Fujian Key Laboratory of Quantum Information and Quantum Optics
(Fuzhou University), Fuzhou 350108, China  Department of Physics, Fuzhou
University, Fuzhou 350108, China  Yang Liu  Fujian Key Laboratory of Quantum
Information and Quantum Optics (Fuzhou University), Fuzhou 350108, China
Department of Physics, Fuzhou University, Fuzhou 350108, China  Zhi-Cheng Shi
Fujian Key Laboratory of Quantum Information and Quantum Optics (Fuzhou
University), Fuzhou 350108, China  Department of Physics, Fuzhou University,
Fuzhou 350108, China  Jie Song  Department of Physics, Harbin Institute of
Technology, Harbin 150001, China  Yan Xia  Fujian Key Laboratory of Quantum
Information and Quantum Optics (Fuzhou University), Fuzhou 350108, China
Department of Physics, Fuzhou University, Fuzhou 350108, China  Shi-Biao Zheng
Fujian Key Laboratory of Quantum Information and Quantum Optics (Fuzhou
University), Fuzhou 350108, China  Department of Physics, Fuzhou University,
Fuzhou 350108, China

###### Abstract

In this work, we exploit the idea of composite pulses to achieve robust
population inversion in a three-level quantum system. The scheme is based on
the modulation of the coupling strength, while the other physical parameters
remain unchanged. The composite pulses sequence is designed by vanishing high-
order error terms, and can compensate the systematic errors to any desired
order. In particular, this scheme keeps a good performance under the
disturbance of waveform deformations. This trait ensures that population
inversion can be nearly obtained even when the pulse sequence has a short jump
delay. As an example, we employ the designed composite pulse sequence to
prepare the $W$ state in a robust manner in the superconducting circuits. The
numerical results show that the fidelity can still maintain a high level in a
decoherence environment.

###### pacs:

02.30.Yy, 03.65.Yz, 31.15.-p, 42.25.Kb

##  I Introduction

Robust control of quantum states is one of the crucial topics in quantum
information processing (QIP) Krause  _et al._ (1998); Ramakrishna  _et al._
(2000); Vilela Mendes and Man’ko (2003); Brańczyk  _et al._ (2007); Wang and
Niu (2020); Zhang  _et al._ (2020a); Dalgaard  _et al._ (2020); Chen and Xue
(2020). Through controlling the external fields, the initial state of a
quantum system can be driven into the desired target state. During this
process, many requirements need to be satisfied. The first one is to quickly
accomplish quantum operations, because the coherence time of quantum systems
is generally very short. Once the operation time is sufficiently long, the
advantages of quantum computation would gradually be lost. The second one is
to perform quantum operations in a robust way. As is well known, the
infidelity of quantum logic gates needs to be below the fault-tolerant
threshold ($10^{-4}$) to protect the reliability of quantum computation
Nielsen and Chuang (2000). However, quantum systems always suffer from the
harmful effect of external noises, which may be caused by the imperfect
knowledge of quantum systems or the fluctuations of external fields, etc.
These external noises finally lead to a sharp decline in the fidelity of
quantum gates. In response, quantum control theory is developed to tackle
these problems.

Quantum control theory has been successfully applied in single Avinadav  _et
al._ (2014); Børkje  _et al._ (2013) and multibody systems Doria  _et al._
(2011); Caneva  _et al._ (2011). At the same time, this theory has made new
progress in quantum state preparations Baksic  _et al._ (2016); Wang  _et al._
(2015); Bao  _et al._ (2018); Chen  _et al._ (2017); Leung  _et al._ (2017);
Wu  _et al._ (2018); Zhang and Yang (2021), quantum gate operations Lin and Li
(2009); Wu  _et al._ (2008); Kang  _et al._ (2018); Kosut  _et al._ (2013);
Hegde  _et al._ (2020); Wu  _et al._ (2020); Guo  _et al._ (2020); Dridi  _et
al._ (2020a); Kang  _et al._ (2020); Chen  _et al._ (2020); Zhang  _et al._
(2020b); Xue  _et al._ (2016), and quantum error corrections Ge  _et al._
(2020); Zhong and Jing (2020); Ma  _et al._ (2020a); Takeda  _et al._ (2020);
Goldman  _et al._ (2020); Hayes  _et al._ (2020). During the quantum control
process, the primary mission is to design an appropriate pulse shape of the
control fields to drive the system evolution as we expect. The traditional
method is to employ the adiabatic passage (AP) Allen and Eberly (1987);
Melinger  _et al._ (1992, 1994); Solá  _et al._ (1999); Band and Magnes
(1994); Vitanov  _et al._ (2001); Goswami (2003), which matches the initial
(target) state to one instantaneous eigenstate as an evolution path of the
system. The AP requests that the physical parameters must vary slowly enough
so as to meet the adiabatic condition. It is widely known that there are two
relevant well-established techniques: the chirped rapid adiabatic passage
(CHIRAP) Melinger  _et al._ (1992, 1994); Solá  _et al._ (1999); Band and
Magnes (1994) and the stimulated Raman adiabatic passage (STIRAP) Gaubatz  _et
al._ (1990); Marte  _et al._ (1991); Shore  _et al._ (1992). The former
proposes to apply a chirped pulse to modulate the frequency detuning of the
system. The slow chirp rate ensures that the adiabatic condition is strictly
satisfied. Therefore, one can drive the system evolution along a given
eigenstate and suppress transitions to other eigenstates Solá  _et al._
(1999). The latter is one of the most popular methods for quantum control in
three-level systems and has since been widely expanded to other fields Vitanov
_et al._ (2017). In the STIRAP, the Stokes pulse and the pump pulse are
applied in a counterintuitive order, and one can achieve a successful
population transfer between two states under the two-photon resonance
condition in three-level systems Vitanov  _et al._ (2017). The main feature of
AP is that the robustness against parameter fluctuations could be achieved at
the expense of operation speed and accuracy. However, once the adiabatic
condition cannot be satisfied well, the accuracy would be sharply dropped due
to the imperfect adiabatic path.

To optimize the speed and accuracy of AP, an improved and acclaimed technique
named shortcut to adiabaticity (STA) Guéry-Odelin  _et al._ (2019); Chen  _et
al._ (2011, 2010) has been developed. There are two common methods in the STA
technique: Lewis-Riesenfeld invariants Lewis and Riesenfeld (1969) and
transitionless quantum driving Demirplak and Rice (2003, 2005); Berry (2009).
The main idea of STA is to design the pulse shape to canceling nonadiabatic
transitions through an additional Hamiltonian. Nevertheless, it depends
heavily on the precisely known physical parameters, and strong fluctuations of
parameters might make this method inefficient. Recently, optimal control (OP)
has become another popular way to achieve high-fidelity quantum operations
Doria  _et al._ (2011); Lin  _et al._ (2020a); Holland  _et al._ (2020); Lin
_et al._ (2020b); Larocca  _et al._ (2020); Tian  _et al._ (2020); An  _et
al._ (2021); Plantenberg  _et al._ (2007); Glaser  _et al._ (2015). One
disadvantage of OP is that the pulse shape is always intricate and nonanalytic
because it is achieved by using advanced mathematical algorithms. As a result,
the design process only yields certain numerical solutions. In order to take
account of the robustness and accuracy of quantum control, one can turn to the
composite pulses (CPs) technique, where the physical parameters are flexible
to design and the shapes are readily implemented in experiments due to its
commonality.

The CPs technique was born in the field of nuclear magnetic resonance Wimperis
(1994); Levitt (1986); Wimperis (1991); Husain  _et al._ (2013). Over the past
decade, CPs have provided a lot of solutions in quantum computation Genov  _et
al._ (2014); Ivanov and Vitanov (2011); Torosov and Vitanov (2019, 2011);
Jones (2013); Genov  _et al._ (2020); Dridi  _et al._ (2020b) due to their
extremely high accuracy, robustness against errors, and extraordinary
flexibility Torosov  _et al._ (2021). Generally speaking, the CPs technique is
composed of a series of precise constant pulses with different relative
phases. These pulses are usually imposed on external laser fields, electric
fields, magnetic fields, or radio-frequency fields, etc. It is shown that the
phase-modulated composite pulses (PMCPs) can not only improve the robustness
against errors (broadband pulse), but also enhance the sensitivity and
selectivity of excitation (narrowband pulse) Torosov  _et al._ (2020a).
Recently, the PMCPs have also been designed to produce a rotation in the Bloch
sphere Torosov and Vitanov (2019), where the errors resulting from the pulse
area are gradually compensated with the increase of the pulse number. Another
application of PMCPs is to implement the NOT gate, and this scheme is robust
against both offset uncertainties and control field variations by a very small
number of modulation parameters Dridi  _et al._ (2020b).

However, when the system phase cannot be modulated or the system does not
carry phase information, PMCPs would be invalid. In order to supplement the
deficiency of composite pulses, researchers turn their attention to some other
adjustable parameters. Therefore, the detuning-modulated composite pulses
(DMCPs) and the strength-modulated composite pulses (SMCPs) have become two
additional promising methods. Recently, based on photonic integrated circuits,
Greener _et al._ applied DMCPs in the coupled optical waveguide model to
achieve complete light transfer Kyoseva  _et al._ (2019). SMCPs are used to
implement dynamically corrected single-qubit gates on singlet-triplet qubits
Wang  _et al._ (2012, 2014), where qubit manipulations are handled by
adjusting the electrically controlled exchange coupling strength. Note that
most works Genov  _et al._ (2014); Torosov and Vitanov (2019); Ivanov and
Vitanov (2011); Torosov and Vitanov (2011); Jones (2013); Genov  _et al._
(2020); Dridi  _et al._ (2020b); Torosov  _et al._ (2021, 2020a); Wang  _et
al._ (2014); Kyoseva  _et al._ (2019); Wang  _et al._ (2012) of CPs are based
on two-level systems. Recently, the CPs have been extended to the three-level
systems and have provided many creative works Torosov and Vitanov (2020); Shi
_et al._ (2021); Torosov  _et al._ (2020b); Mansourzadeh-Ashkani  _et al._
(2021); Ghosh  _et al._ (2017); Friesen  _et al._ (2017); Shi  _et al._
(2022), such as the implementation of high-fidelity composite quantum gates
Torosov and Vitanov (2020); Shi  _et al._ (2021) and the efficient detection
of chiral molecules Torosov  _et al._ (2020b). It is worth mentioning that
those works Torosov and Vitanov (2020); Shi  _et al._ (2021); Torosov  _et
al._ (2020b) still focus on the phase modulation in three-level systems.

In this paper, based on the composite pulses technique, we develop a general
SMCPs scheme for robust population inversion in three-level systems. According
to the Taylor expansion, the final transition probability expression is
rearranged into a series of error terms. With predetermined phases and pulse
area, these error terms could be effectively eliminated from low order to high
order by properly adjusting coupling strengths. In addition, we investigate
the influence of waveform deformation on the transition probability. The
results show that slight deformation still allows the SMCPs to work well,
except for the heavily deformed case. In the end, we further apply the current
SMCPs scheme to prepare the $W$ state in the superconducting quantum
interference device (SQUID) model, and the numerical simulations demonstrate
that the SMCPs can still maintain a robust performance in a decoherence
environment.

The structure of this paper is organized as follows. In Sec. II, we present
the design procedure of the SMCPs in the three-level system, which are
directly derived from the total composite propagator. In Sec. III, we
illustrate how to eliminate the effect of pulse area error by SMCPs with the
specific number of pulses. Here, the sequence of up to seven pulses has been
studied, and the longer sequence could also be obtained in a similar way.
Then, we analyze the impact of the waveform deformation of CPs on the
performance of the transition probability and compare our scheme with other
composite pulses. In Sec. IV, we give the application of the SMCPs scheme in
robustly preparing the $W$ state in the superconducting circuits. Finally, we
give the conclusion in Sec. V.

##  II Theoretical model

Figure 1: A $\Lambda$-type three-level system driven by two control fields in
the two-photon resonance regime.

In this section, we elaborate on the design procedure of the general pulse
waveform in a three-level system with $\Lambda$-type structure. This system
has been regarded as a paradigmatic model in many branches of physics,
including atomic and molecular Pillet  _et al._ (1993); Featonby  _et al._
(1996); Král and Shapiro (2001), quantum information science Lacour  _et al._
(2006), and some other fields Klein  _et al._ (2007, 2008); Longhi (2009);
Wang and Clerk (2012). The basic structure is shown in Fig. 1, where
$|g\rangle$ and $|f\rangle$ are two ground states and $|e\rangle$ is the
excited state. There are two control fields that drive two transitions:
$|g\rangle\leftrightarrow|e\rangle$ and $|f\rangle\leftrightarrow|e\rangle$.
Meanwhile, the transition frequency of the three-level system and the carrier
frequency of the control fields satisfy the two-photon resonance condition
Scully and Zubairy (1997). Note that the direct transition between two ground
states is forbidden. In the interaction picture, the dynamics of the system
can be described by the following Hamiltonian ($\hbar=1$ hereafter):

$\displaystyle H=\Delta|e\rangle\langle e|+\frac{1}{2}\Big{(}\Omega
e^{i\phi}|g\rangle\langle e|+\lambda e^{i\varphi}|f\rangle\langle
e|+\mathrm{H.c.}\Big{)},$ (1)

where $\Omega$ and $\phi$ ($\lambda$ and $\varphi$) are the coupling strength
and the phase of the transition $|g\rangle\leftrightarrow|e\rangle$
($|f\rangle\leftrightarrow|e\rangle$), and $\Delta$ is the single-photon
detuning. When the detuning is much larger than two coupling strengths (i.e.,
$|\Delta|\gg\Omega,\lambda$), the excited state $|e\rangle$ could be
adiabatically eliminated Brion  _et al._ (2007). Then, the three-level system
would be reduced as an effective two-level system. In the following, we study
the system dynamics including the excited state, and thus do not focus on this
situation.

The propagator $U$ of this three-level system satisfies the following
$\mathrm{Schr\ddot{o}dinger}$ equation:

$\displaystyle i\dot{U}=HU.$ (2)

When the Hamiltonian given by Eq. (1) is time independent, the solution of
this equation becomes $U=\mathrm{exp}\big{(}{-iHT}\big{)}$, where $T$ is the
time duration. It is instructive to adopt its matrix form in the basis
$\\{|g\rangle,|f\rangle,|e\rangle\\}$, which reads:

$\displaystyle U(\Theta)=\left[\begin{array}[]{ccc}U_{11}&U_{12}&U_{13}\\\
U_{21}&U_{22}&U_{23}\\\ U_{31}&U_{32}&U_{33}\\\ \end{array}\right],$ (6)

where

$\displaystyle U_{11}$ $\displaystyle=$
$\displaystyle\cos^{2}{\Theta}+e^{-i\delta}(\cos{\frac{A}{2}}+\frac{i\Delta}{\sqrt{1+\Delta^{2}}}\sin{\frac{A}{2}})\sin^{2}{\Theta},$
$\displaystyle U_{22}$ $\displaystyle=$
$\displaystyle\sin^{2}{\Theta}+e^{-i\delta}(\cos{\frac{A}{2}}+\frac{i\Delta}{\sqrt{1+\Delta^{2}}}\sin{\frac{A}{2}})\cos^{2}{\Theta},$
$\displaystyle U_{33}$ $\displaystyle=$ $\displaystyle
e^{-i\delta}(\cos{\frac{A}{2}}-\frac{i\Delta}{\sqrt{1+\Delta^{2}}}\sin{\frac{A}{2}}),$
$\displaystyle U_{12}$ $\displaystyle=$
$\displaystyle\frac{1}{2}e^{-i\delta}e^{i\Psi}(\cos{\frac{A}{2}}+\frac{i\Delta}{\sqrt{1+\Delta^{2}}}\sin{\frac{A}{2}}-e^{i\delta})\sin{2\Theta},$
$\displaystyle U_{21}$ $\displaystyle=$
$\displaystyle\frac{1}{2}e^{-i\delta}e^{-i\Psi}(\cos{\frac{A}{2}}+\frac{i\Delta}{\sqrt{1+\Delta^{2}}}\sin{\frac{A}{2}}-e^{i\delta})\sin{2\Theta},$
$\displaystyle U_{13}$ $\displaystyle=$
$\displaystyle-e^{-i2\delta}U_{31}^{*}=-ie^{-i\delta}e^{i\phi}\frac{\sin{\Theta}}{\sqrt{1+\Delta^{2}}}\sin{\frac{A}{2}},$
$\displaystyle U_{23}$ $\displaystyle=$
$\displaystyle-e^{-i2\delta}U_{32}^{*}=-ie^{-i\delta}e^{i\varphi}\frac{\cos{\Theta}}{\sqrt{1+\Delta^{2}}}\sin{\frac{A}{2}}.$

Here, $\Theta=\arctan{(\Omega/\lambda)}$, $\delta=\int_{0}^{T}\Delta/2\ dt$,
the phase difference $\Psi=\phi-\varphi$, and the expression of the total
pulse area $A$ is

$\displaystyle A=\int_{0}^{T}\sqrt{\Omega^{2}+\lambda^{2}+\Delta^{2}}\ dt.$
(7)

Obviously, when the pulse area $A=2\pi$, two ground states are decoupled from
the excited state $|e\rangle$ so that we can only concentrate on the subspace
$\\{|g\rangle,|f\rangle\\}$. Then, when the system is initially in the ground
state $|g\rangle$, the transition probability $P$ of the target state
$|f\rangle$ becomes

$\displaystyle
P=|U_{21}|^{2}=\cos^{2}{\left(\frac{\pi\Delta}{2\sqrt{1+\Delta^{2}}}\right)}\sin^{2}{2\Theta}.$
(8)

It is easily found from Eq. (8) that the detuning $\Delta$ is required to be
zero in order to achieve complete population inversion for the single pulse.

Then, one can see from Eq. (7) that two kinds of uncertainties would cause the
pulse area error when $\Delta=0$. The first one is the coupling strength
error, which may be caused by the inhomogeneity of the control fields. Another
one is the inaccurate pulse duration, which may be limited to the experiment
condition or manual operation accuracy. Note that tiny errors in the pulse
area would have a significant effect on population inversion. This can be
demonstrated as follows. Suppose that the pulse area error is represented by
$\epsilon$, and a deviation in the original pulse area reads
$A=A(1+\epsilon)$. As a result, the transition probability given by Eq. (8) is
rewritten as

$P(\epsilon)=\sin^{4}{\frac{A(1+\epsilon)}{4}}\sin^{2}{2\Theta}.$ (9)

In order to study how the pulse area error makes an impact on population
inversion, the transition probability $P(\epsilon)$ given by Eq. (9) is
expanded as the Taylor series,

$\begin{split}&P(\epsilon)=\alpha_{0}+\alpha_{1}\epsilon+\alpha_{2}\epsilon^{2}+O(\epsilon^{3}),\\\
\end{split}$ (10)

where $\alpha_{j}$ is the $j\mathrm{th}$-order coefficient, $j=0,1,2,\cdots$.
Actually, $\alpha_{0}$ is the transition probability in the absence of the
pulse area error. By setting the pulse area $A=2\pi(1+\epsilon)$, the first
three order coefficients are

$\begin{split}\alpha_{0}=\sin^{2}{2\Theta},~{}\alpha_{1}&=0,~{}\alpha_{2}=-\frac{1}{2}\pi^{2}\sin^{2}{2\Theta}.\\\
\end{split}$ (11)

Notice that the first-order coefficient $\alpha_{1}$ is automatically
eliminated because all odd-order coefficients contain the term $\sin{(A/2)}$.
Here, $\Theta=\pi/4+k\pi/2$ is the indispensable condition to achieve
population inversion, i.e., $\alpha_{0}=1$. We can observe from Eq. (11) that
the second-order coefficient $\alpha_{2}$ has the same monotony with
$\alpha_{0}$ because they have the same term $\sin^{2}{2\Theta}$. This means
that the error would also reach maximum if we achieve complete population
inversion.

Figure 2: (a) Transition probability $P_{N}(\epsilon)$ vs the pulse area error
$\epsilon$ under different pulses sequences. (b) Inversion error
$1-P_{N}(\epsilon)$ vs the pulse area error $\epsilon$. All $\Theta_{n}$ come
from Table 1.

In Figs. 2 and 2, the red solid curve represents $P(\epsilon)$ as a function
of the pulse area error $\epsilon$, which shows that the transition
probability sharply drops near $\epsilon=0$. Hence, the single-pulse scheme
does not possess a robust manner against the pulse area error.

In order to solve this problem, next we put forward the SMCPs scheme. The
sequence is composed of multiple single pulses, where only coupling strengths
for each pulse are different. For the compact expression, we label the
Hamiltonian and the propagator of the $n\mathrm{th}$ pulse as $H_{n}$ and
$U_{n}(\Theta_{n})$, respectively. Then, the general form of the total
propagator for the $N$-pulses sequence can be written as

$\displaystyle U^{(N)}\\!=\\!U_{N}(\Theta_{N})U_{N-1}(\Theta_{N-1})\cdots
U_{2}(\Theta_{2})U_{1}(\Theta_{1}).$ (12)

The detailed derivation of the total propagator $U^{(N)}$ without the pulse
area error is presented in Appendix A. Note that the transition probability of
every state is associated with both the detuning and the coupling strengths of
each pulse. For simplicity, we assume each pulse is in the resonant regime
($\Delta_{n}=0$) in the following.

The target of this work is to achieve robust population inversion by SMCPs.
Without loss of generality, we assume the initial state is
$|\psi_{i}\rangle=|g\rangle$; then the target state becomes
$|\psi_{f}\rangle=|f\rangle$. Each pulse area is represented as
$2\pi(1+\epsilon)$, and we label $P_{N}(\epsilon)$ as the transition
probability of the state $|f\rangle$ in the $N$ pulses sequence. Similar to
the derivations in Eq. (10), $P_{N}(\epsilon)$ can be written in the following
form by the Taylor expansion:

$P_{N}(\epsilon)=\alpha_{N,0}+\alpha_{N,1}\epsilon+\alpha_{N,2}\epsilon^{2}+\cdots+O(\epsilon^{n}),$
(13)

where $\alpha_{N,j}$ is the $j\mathrm{th}$-order coefficient in the $N$ pulses
sequence, and its expression can be deduced from the total propagator given by
Eq. (12).

Here, we do not intend to change the pulse area and phases of this system, and
thus the pulse area and phases are predetermined by constants during the
evolution process. Note that all odd-order coefficients $\alpha_{N,j}(j=2n-1)$
would vanish if the pulse area of each pulse is equal to $2\pi$, i.e.,
$A_{1}=\cdots=A_{N}=2\pi$. As a result, only even-order coefficients are left
in Eq. (13). By modulating the parameters $\Theta_{n}$ ($n=1,\dots,N$), we
demand that the zeroth-order coefficient be unity, which ensures complete
population inversion, and other even-order coefficients vanish as many as
possible. Note that according to the expression
$\Theta_{n}=\arctan{(\Omega_{n}/\lambda_{n})}$, we can keep the coupling
strength $\lambda_{n}$ unchanged and alter the coupling strength $\Omega_{n}$
to implement the modulation of $\Theta_{n}$. This means that we only adjust
one of the coupling strengths to achieve robust population inversion in the
three-level system.

To be specific, if the sequence is composed of $N$ pulses, we first solve $N$
equations to obtain the solutions of $\Theta_{n}$, i.e.,

$\displaystyle\alpha_{N,0}=1,$ $\displaystyle\alpha_{N,2}=0,$ (14)
$\displaystyle~{}~{}~{}\cdots$ $\displaystyle\alpha_{N,2N-2}=0.$

Then, we calculate different coupling strengths $\Omega_{n}$ of each pulse
according to the expression $\Omega_{n}=\lambda_{n}\tan\Theta_{n}$.
Consequently, we have the composite pulses sequence to achieve population
inversion. In this situation, the transition probability is accurate to the
order of $O(\epsilon^{2N-2})$. Compared with the single-pulse scheme, the
transition probability in the SMCPs scheme has a better robust manner against
the pulse area error $\epsilon$ since the high-order coefficients disappear.
In other words, eliminating more high-order coefficients in Eq. (13) could
further suppress the detrimental effect of the pulse area error. When the
composite pulses sequence is sufficiently short, the analytical results are
readily obtained by solving the corresponding equations. However, the long
pulse sequence would increase the number of equations. As a result, it is
difficult to derive the analytical results. Instead, we provide the numerical
solutions for the long pulse sequence. In order to demonstrate this issue in
more detail, we exemplify the SMCPs scheme with a specific number of pulses
below.

##  III Strength-modulated composite pulses scheme

In this section, we study how to design the SMCPs to achieve robust population
inversion when the pulse area and phases remain unchanged. For simplicity, we
uniformly set $A_{n}=2\pi$, $\phi_{n}=\pi/2$, $\varphi_{n}=0$, and assume the
system always works in the resonance regime (i.e., $\Delta_{n}=0$). Note that
$\Theta_{n}$ are only the adjustable parameters.

###  III.1 Two pulses

The propagator in the two-pulses sequence can be expressed by

$\displaystyle U^{(2)}=U_{2}(\Theta_{2})U_{1}(\Theta_{1}).$ (15)

According to the propagator $U^{(2)}$ given by Eq. (15), the transition
probability $P_{2}(\epsilon)$ of the state $|f\rangle$ is

$\displaystyle P_{2}(\epsilon)$ $\displaystyle=$
$\displaystyle\frac{1}{4}\Big{\\{}\sin^{2}{\\!\pi\epsilon}\big{[}\cos{\\!\Theta_{1}}\\!\sin{\\!\Theta_{2}}\\!+\\!\cos{\\!\Theta_{2}}(2\sin{\\!\Theta_{1}}\\!+\\!\sin{\\!\Theta_{2}})\big{]}$
(18)
$\displaystyle\\!-\\!2\cos^{4}{(\frac{\pi\epsilon}{2})}\sin{2(\Theta_{1}\\!-\\!\Theta_{2})}\Big{\\}}^{2}.$

The detailed derivation process is given in Appendix B. By the Taylor
expansion, $P_{2}(\epsilon)$ can be regrouped as

$P_{2}(\epsilon)=\alpha_{2,0}+\alpha_{2,2}\epsilon^{2}+O(\epsilon^{4}),$ (19)

where the expressions of the first two coefficients read

$\displaystyle\alpha_{2,0}$ $\displaystyle=$
$\displaystyle\sin^{2}\big{[}2(\Theta_{1}-\Theta_{2})\big{]},$ (20)
$\displaystyle\alpha_{2,2}$ $\displaystyle=$ $\displaystyle\pi^{2}\sin
2(\Theta_{2}-\Theta_{1})\Big{[}\sin 2(\Theta_{1}-\Theta_{2})$
$\displaystyle+2\cos\Theta_{2}\sin\Theta_{1}+\frac{1}{2}\sin
2\Theta_{1}+\frac{1}{2}\sin\Theta_{2}\Big{]}.$

In the two-pulses sequence, we need to simultaneously solve the following
equations:

$\displaystyle\alpha_{2,0}=1,$ (24) $\displaystyle\alpha_{2,2}=0.$ (25)

After some calculations, one solution of Eqs. (24) and (25) can be written as

$\begin{split}\Theta_{1}=m\pi+\frac{\pi}{8},~{}\Theta_{2}=l\pi+\frac{3\pi}{8},\end{split}$
(26)

where $m$ and $l$ are arbitrary integers. By this choice of $\Theta_{1}$ and
$\Theta_{2}$, the transition probability is accurate to the fourth order in
the pulse area error $\epsilon$,

$P_{2}(\epsilon)=1+O(\epsilon^{4}).$ (27)

The gray solid curve in Fig. 2 represents the relation between the transition
probability and the pulse area error when $\Theta_{1}=\pi/8$ and
$\Theta_{2}=3\pi/8$. The result shows that the transition probability appears
as a small flat top profile against the pulse area error in the two-pulses
sequence. Correspondingly, the gray solid curve in Fig. 2 demonstrates that
the inversion error in the two-pulses sequence is much lower than that in the
single pulse under the same condition. Obviously, the two-pulses sequence has
a better fault tolerance than the single pulse in the error-prone environment.

Note that the above designed sequence can also used for implementing a NOT
gate with a well-defined phase. Although the phases are inessential for
population inversion, they play a very important role when achieving a quantum
gate because different phases would determine different types of quantum
gates. For our objective, by the group of solution $\Theta_{1}=\pi/8$ and
$\Theta_{2}=3\pi/8$ and the phase difference $\Psi_{n}=\pi/2$, the final
propagator in the basis $\\{|g\rangle,|f\rangle\\}$ becomes (up to a global
phase)

$\displaystyle U^{(2)}=\left[\begin{array}[]{ccc}0&1\\\ 1&0\\\
\end{array}\right],$

which is actually the $X$ gate. However, when the phase difference is chosen
as $\Psi_{n}=\pi$, the final propagator would become

$\displaystyle U^{(2)}=\left[\begin{array}[]{ccc}0&-i\\\ i&0\\\
\end{array}\right],$

which is the $Y$ gate; see Appendix B for details.

Furthermore, this current sequence is readily extended to yield an arbitrary
rotation gate. The general form of an arbitrary rotation gate in the basis
$\\{|g\rangle,|f\rangle\\}$ is Nielsen and Chuang (2000)

$\displaystyle\hat{R}=\left[\begin{array}[]{cc}\cos{(\theta/2)}e^{i\beta^{\prime}}&\sin{(\theta/2)e^{i\beta}}\\\
-\sin{(\theta/2)}e^{-i\beta}&\cos{(\theta/2)}e^{-i\beta^{\prime}}\\\
\end{array}\right],$

where $\theta$ is the rotation angle ($0\leq\theta\leq\pi$), and $\beta$ and
$\beta^{\prime}$ are relative phases. To implement this rotation gate, Eq.
(24) needs to be modified as $\alpha_{2,0}=\sin^{2}{(\theta/2)}$, and then we
properly adjust the phases of the coupling strengths. As an example, to
implement the Hadamard gate by the two-pulses sequence, we first need to
obtain a maximum superposition state of $|g\rangle$ and $|f\rangle$, and the
equations read

$\displaystyle\alpha_{2,0}=1/2,$ $\displaystyle\alpha_{2,2}=0.$

One group of the solutions is $\Theta_{1}=0.1047$ and $\Theta_{2}=0.4974$.
Then, we need to choose two proper phases $\phi_{n}=\varphi_{n}$, and obtain
the Hadamard gate

$\displaystyle H=\frac{1}{\sqrt{2}}\left[\begin{array}[]{cc}1&1\\\ -1&1\\\
\end{array}\right].$

Note that different from the previous work Torosov and Vitanov (2020) for the
Hadamard gate in the three-level system, the total pulse area of the current
sequence can be $4\pi$, while it is $6\pi$ for the sequence $\mathcal{H}_{6}$
in Ref. Torosov and Vitanov (2020). This flexibility of the sequence length
benefits from the current sequence starting straight from the three-level
system rather than nesting the existing sequence.

###  III.2 Three pulses

The three-pulses sequence produce the following propagator:

$\displaystyle U^{(3)}=U_{3}(\Theta_{3})U_{2}(\Theta_{2})U_{1}(\Theta_{1}).$
(32)

Thus, the exact expression of the transition probability $P_{3}(\epsilon)$ can
be written as

$\displaystyle P_{3}(\epsilon)$ $\displaystyle=$
$\displaystyle\Big{\\{}\\!\sin\\!2\Theta_{3}(\sin\\!2\Theta_{1}\sin\\!2\Theta_{2}\\!-\\!2\sin\\!\Theta_{1}\sin\\!\Theta_{2}\sin\\!\pi\epsilon\sin\\!2\pi\epsilon)\\!+\\!8\cos^{2}\\!\Theta_{1}\\!\big{[}\\!\cos\\!\Theta_{2}\sin\\!\Theta_{3}\cos(\Theta_{2}\\!-\\!\Theta_{3})\\!+\\!\sin\\!\Theta_{2}\cos\\!\Theta_{3}\sin^{2}\\!\pi\epsilon\big{]}\\!$
$\displaystyle+4\sin^{2}\\!\Theta_{3}\\!\big{[}\sin\\!2\Theta_{1}\sin^{2}\\!\Theta_{2}\\!+\\!2\sin\\!\Theta_{1}\cos\\!\Theta_{2}\sin^{2}\\!\pi\epsilon\big{]}\\!-\\!4\sin^{2}\\!\pi\epsilon(\sin\\!2\Theta_{1}\\!\cos\\!\Theta_{2}\\!\cos\\!\Theta_{3}\\!-\\!4\sin\\!\Theta_{1}\sin\\!\Theta_{2}\sin\\!2\Theta_{3})\\!$
$\displaystyle-\cos\pi\epsilon\big{[}\\!\sin 2\Theta_{1}\\!+\\!\sin
2(\Theta_{1}\\!-\\!\Theta_{2})\\!+\\!\sin
2(\Theta_{1}\\!-\\!\Theta_{3})\\!-\\!3\sin\\!2(\Theta_{1}\\!-\\!\Theta_{2}\\!+\\!\Theta_{3})\\!+\\!4\sin\Theta_{2}\cos\Theta_{3}\cos(\Theta_{2}\\!-\\!\Theta_{3})$
$\displaystyle+8\sin\Theta_{1}\cos\Theta_{3}\sin^{2}\pi\epsilon(\cos(\Theta_{1}\\!-\\!\Theta_{2})\\!+\\!\cos\Theta_{2}\\!\cos\Theta_{3}\\!+\\!1)\\!\big{]}\\!-\\!\cos^{2}\\!\pi\epsilon\big{[}\\!\sin
2\Theta_{1}\\!+\\!\sin 2\Theta_{2}\\!+\\!\sin 2\Theta_{3}\\!-\\!\sin
2(\Theta_{1}\\!-\\!\Theta_{2})\\!$ $\displaystyle-\\!\sin
2(\Theta_{1}\\!-\\!\Theta_{3})\\!-\\!\sin
2(\Theta_{2}\\!-\\!\Theta_{3})\\!-\\!3\sin
2(\Theta_{1}\\!-\\!\Theta_{2}\\!+\\!\Theta_{3})\\!\big{]}\\!+\\!8\cos^{3}\\!\pi\epsilon\sin\\!\Theta_{1}\cos\Theta_{3}\cos(\Theta_{1}\\!-\\!\Theta_{2})\cos(\Theta_{2}\\!-\\!\Theta_{3})\Big{\\}}^{2}.$

Similarly, through the Taylor expansion, the transition probability
$P_{3}(\epsilon)$ can be expanded as

$P_{3}(\epsilon)=\alpha_{3,0}+\alpha_{3,2}\epsilon^{2}+\alpha_{3,4}\epsilon^{4}+O(\epsilon^{6}).$
(33)

Since only three parameters $\\{\Theta_{1},\Theta_{2},\Theta_{3}\\}$ are
contained in the three-pulses sequence, it is sufficient to extend Eq. (13) to
the fourth-order term. The expressions of the first three coefficients are

$\displaystyle\alpha_{3,0}$ $\displaystyle\\!=\\!$
$\displaystyle\sin^{2}{2(\Theta_{1}-\Theta_{2}+\Theta_{3})},$
$\displaystyle\alpha_{3,2}$ $\displaystyle\\!=\\!$
$\displaystyle-\frac{\pi^{2}}{2}\sin{2(\Theta_{1}\\!-\\!\Theta_{2}\\!+\\!\Theta_{3})}\Big{\\{}\sin{2(\Theta_{1}\\!-\\!\Theta_{2})}\\!+\\!2\cos{(\Theta_{1}\\!-\\!\Theta_{2})}\sin{(\Theta_{1}\\!+\\!\Theta_{2}\\!-\\!2\Theta_{3})}\\!+\\!3\sin{2(\Theta_{1}\\!-\\!\Theta_{2}\\!+\\!\Theta_{3})}\\!$
(35)
$\displaystyle+4\sin{\Theta_{1}}\cos{(\Theta_{2}\\!-\\!2\Theta_{3})}\\!+\\!4\cos{\Theta_{3}}\big{[}\sin{(2\Theta_{1}\\!-\\!\Theta_{2})}\\!+\\!\sin{\Theta_{1}}\big{]}\Big{\\}},$
$\displaystyle\alpha_{3,4}$ $\displaystyle\\!=\\!$
$\displaystyle\frac{\pi^{4}}{48}\Big{\\{}\\!48\big{[}\sin\\!\Theta_{1}\cos(\Theta_{2}\\!-\\!2\Theta_{3})\\!+\\!\sin(2\Theta_{1}\\!-\\!\Theta_{2})\cos\Theta_{3}+\sin\\!\Theta_{1}\cos\Theta_{3}\\!+\\!\cos(\Theta_{1}\\!-\\!\Theta_{2})\sin(\Theta_{1}\\!-\\!\Theta_{3})\cos(\Theta_{2}\\!-\\!\Theta_{3})$
$\displaystyle+51\sin^{2}\\!2(\Theta_{1}\\!-\\!\Theta_{2}\\!+\\!\Theta_{3})\big{]}^{2}\\!\\!+\\!2\sin\\!2(\Theta_{1}\\!-\\!\Theta_{2}\\!+\\!\Theta_{3})\big{[}64\sin\\!\Theta_{1}\cos(\Theta_{2}\\!-\\!2\Theta_{3})\\!+\\!4\cos\\!\Theta_{3}\big{[}16\sin(2\Theta_{1}\\!-\\!\Theta_{2})\\!+\\!19\sin\\!\Theta_{1}$
$\displaystyle+3\sin\\!\Theta_{2}\big{]}\\!+\\!64\cos(\Theta_{1}\\!-\\!\Theta_{2})\sin(\Theta_{1}\\!-\\!\Theta_{3})\cos(\Theta_{2}\\!-\\!\Theta_{3})\\!+\\!12\sin\\!\Theta_{1}\\!\cos\\!\Theta_{2}\\!+\\!3(\sin\\!2\Theta_{1}\\!+\\!\sin\\!2\Theta_{2}\\!+\\!\sin\\!2\Theta_{3})\big{]}\Big{\\}}.$

Therefore, in the three-pulses sequence, the solutions need to simultaneously
satisfy the following equations:

$\displaystyle\alpha_{3,0}=1,$ (38) $\displaystyle\alpha_{3,2}=0,$ (39)
$\displaystyle\alpha_{3,4}=0.$ (40)

Here, the solution for Eq. (38) can be written as
$\Theta_{1}\\!-\\!\Theta_{2}\\!+\\!\Theta_{3}\\!=\\!\pi/4+m\pi/2$, where $m$
is an arbitrary integer. However, it is difficult to derive analytical
formulas for Eqs. (39) and (40). Instead, we employ the numerical method to
obtain the solutions of Eqs. (38)–(40). One numerical solution is presented in
Table 1. We figure out the profile of the transition probability
$P_{3}(\epsilon)$ by the blue dashed curves in Figs. 2 and 2. Clearly, the
three-pulses sequence produces a wider top platform than the two-pulses
sequence because the fourth-order coefficient $\alpha_{34}$ vanishes. In this
situation, the transition probability is accurate to the sixth order in the
pulse area error $\epsilon$,

$P_{3}(\epsilon)=1+O(\epsilon^{6}).$ (41)

Table 1: Values of the coupling strength ratio to achieve population inversion. $N$ | 2 | 3 | 4 | 5 | 6 | 7
---|---|---|---|---|---|---
$\Theta_{1}$ | 0.3927 | 0.1944 | 0.0897 | 0.0477 | 0.0203 | 0.0137
$\Theta_{2}$ | 1.1781 | 0.7854 | 0.4688 | 0.2913 | 0.1508 | 0.1072
$\Theta_{3}$ | – | 1.3764 | 1.0613 | 0.7854 | 0.5004 | 0.3811
$\Theta_{4}$ | – | – | 1.4676 | 1.2795 | 0.9948 | 0.8233
$\Theta_{5}$ | – | – | – | 1.5231 | 1.3826 | 1.2459
$\Theta_{6}$ | – | – | – | – | 1.5431 | 1.4854
$\Theta_{7}$ | – | – | – | – | – | 1.5606

###  III.3 More than three pulses

Next, we briefly present the multiple-pulses design procedure. The derivations
of the sequence with larger pulse numbers are similar to the
two-(three)-pulses sequence. At first, the total propagator of the $N$-pulses
sequence is

$\displaystyle U^{(N)}=U_{N}(\Theta_{N})U_{N-1}(\Theta_{N-1})\dots
U_{1}(\Theta_{1}).$ (42)

In the $N$-pulses sequence, the transition probability could be expanded to
$(2N)\mathrm{th}$ order by the Taylor expansion

$P_{N}(\epsilon)=\alpha_{N,0}+\alpha_{N,2}\epsilon^{2}+\alpha_{N,4}\epsilon^{4}+\dots+O(\epsilon^{2N}).$
(43)

It is not hard to find that the zeroth coefficient $\alpha_{N,0}$ can be
expressed by (see Appendix A for details)

$\alpha_{N,0}=\sin^{2}{\bigg{[}2\sum^{N}_{k=1}(-1)^{k+1}\Theta_{k}\bigg{]}}.$
(44)

Therefore, the solution of the equation $\alpha_{N,0}=1$ can be written as

$\displaystyle\sum^{N}_{k=1}(-1)^{k+1}\Theta_{k}=\frac{\pi}{4}+\frac{m\pi}{2},$

where $m$ is an arbitrary integer.

Then, the remaining $(N-1)$ equations to be solved come from Eqs. (II). In
this situation, the maximum transition probability is accurate to the
$(2N)\mathrm{th}$ order in the pulse area error $\epsilon$,

$P_{N}(\epsilon)=1+O(\epsilon^{2N}).$ (45)

We present some numerical solutions of $\Theta_{n}$ in Table 1, and plot in
Figs. 2 and 2 the performance of the maximum transition probability produced
by the four pulses to the seven pulses. As expected, the more pulses we apply,
the broader region of the high fidelity against the pulse area error we will
get. For example, the inversion error $1-P_{N}(\epsilon)$ still remains below
$10^{-4}$ even though the pulse area error reaches 0.328 in the seven-pulses
sequence.

###  III.4 The influence of waveform deformation

In the ideal situation, the waveform of CPs is the perfect square wave.
Unfortunately, the limitations of experimental conditions would produce
imperfect pulse shapes in practice. The main problem is the waveform
deformation, which usually takes place at the moment of the square wave
switching, reflected in the strength jump delay, i.e., the rise and fall
edges. As a result, the quantum system is approximately driven by a time-
dependent pulse instead of constant pulses sequence. Here, taking the three-
pulses sequence for example, we study this issue by designing a function to
simulate the imperfect square pulse waveform, and the expression of the
coupling strength is

$\displaystyle\Omega(t)=\left\\{\begin{array}[]{ll}\displaystyle\Omega_{2}\\!-\\!\frac{\Omega_{2}\\!-\\!\Omega_{1}}{1\\!+\\!e^{\tau(t\\!-\\!t_{1})}},&\displaystyle~{}~{}~{}t\leq\frac{T\\!+\\!t_{1}\\!-\\!t_{3}}{2},\\\\[12.05553pt]
\displaystyle\Omega_{3}\\!-\\!\frac{\Omega_{3}\\!-\\!\Omega_{2}}{1\\!+\\!e^{\tau(t\\!-\\!t_{1}\\!-\\!t_{2})}},&\displaystyle~{}~{}~{}\frac{T\\!+\\!t_{1}\\!-\\!t_{3}}{2}<t\leq
T,{}{}\\\ \end{array}\right.$ (48)

where $\tau$ is called the waveform deformation parameter and
$t_{n}=2\pi/\sqrt{\Omega_{n}^{2}+\lambda_{n}^{2}}$, $n=1,2,3$.

Figure 3: (a) The waveform of the coupling strength $\Omega(t)$ with different
deformation degree in the three-pulses sequence. (b) The maximum transition
probability $P_{3}(\epsilon)$ vs the deformation parameter $\tau$ and the
pulse area error $\epsilon$. All $\Theta_{n}$ come from Table 1.

It is easily found from Eq. (48) that one can control the width of the rising
and falling edges by adjusting the dimensionless parameter $\tau$. When the
value of $\tau$ is larger, the sequence given by Eq. (48) is much closer to
the constant pulses sequence. For example, it is shown by the red solid curve
in Fig. 3 that the deformation degree of the square wave is already negligible
when $\tau=100$. When $\tau\rightarrow+\infty$, Eq. (48) describes the perfect
square wave. However, when the value of $\tau$ is sufficiently small, the
waveform described by Eq. (48) is quite different from the perfect square
wave, as shown by the blue dotted curve in Fig. 3. We plot in Fig. 3 the
transition probability as a function of the pulse area error and the waveform
deformation. The results demonstrate that a slight waveform deformation has
little influence on our scheme since the width of the high transition
probability region almost remains unchanged when $\tau>10$. A noticeable
change appears in the transition probability only when the waveform
deformation is very severe. However, the transition probability still
maintains at a high level. Hence, our scheme is also robust against the
waveform deformation.

###  III.5 Comparison with other composite pulses

In this section, we make a brief comparison of our CPs sequence with some
previous works. First, most previous works Husain  _et al._ (2013); Torosov
and Vitanov (2019); Jones (2013); Torosov and Vitanov (2011, 2020) are based
on phase modulation, in which the phase is a time-dependent variable
parameter. Different from the phase modulation, here we propose a distinct
modulation method: coupling strength modulation, where all phases of the
system remain unchanged during the evolution process. Second, the detuning
modulation is equivalent to the strength modulation, in principle, since two
categories of the Hamiltonian can be associated through a rotation
transformation Shevchenko  _et al._ (2010).

The common starting point of the current sequence and the sequence in Ref.
Torosov and Vitanov (2020) is to provide robustness against the pulse area
error in the three-level system. With the same total pulse area, the
robustness of our composite sequence agrees with that of the composite
sequence in Ref. Torosov and Vitanov (2020). However, there are some
differences between them. At first, the structure of the sequences is
distinct. We have only three pulses, while there are six pulses in Ref.
Torosov and Vitanov (2020) with the same pulse area. From the perspective of
the waveform, our sequence is a simpler one because it contains fewer pulse
jumps. Furthermore, the design procedures are quite different. The sequence in
Ref. Torosov and Vitanov (2020) is constructed by the combination of a
preexisting CPs sequence in two-level systems Torosov and Vitanov (2011),
ignoring the dynamics of the excited state. Here, we consider the dynamics of
all states in the three-level system, and directly design the composite pulses
sequence according to the propagator. Finally, the application systems are
different, since we use strength modulation here, while the sequence in Ref.
Torosov and Vitanov (2020) relies on phase modulation. To be specific, one
could adopt the phase modulation in some systems with easily adjustable phases
(e.g., the trapped ions system Ivanov and Vitanov (2011)) while the strength
modulation could be conveniently applied in systems with uncontrollable phases
(e.g., the integrated photonic circuits Politi  _et al._ (2008)).

For the strength modulation, the difference between the two-level system and
the three-level system is mainly embodied in the number of modulation
parameters. There are two adjustable coupling strengths in the three-level
system, while there is only single coupling strength in the two-level system.
When the pulse area is fixed, one cannot eliminate the error by only
modulating the coupling strength in a two-level system Torosov and Vitanov
(2019). In order to make this strength modulation valid, one must also
modulate the pulse area, as done in Ref. Wang  _et al._ (2014). However, this
situation is quite different in the three-level system. Without a change in
the pulse area, we can modulate the ratio of two coupling strengths to
effectively eliminate the systematic errors.

Moreover, the SMCPs in Ref. Wang  _et al._ (2014) focuses on altering the
exchange coupling to achieve robust quantum control in the singlet-triplet
qubit. This kind of pulse sequence needs to nest amounts of short pulses.
Hence, the pulse length would be long and may not maintain a high-fidelity
performance in the decoherence environment. In our scheme, the pulse length
can be appropriately chosen to equilibrate between pulse duration and
robustness. On the other hand, the SMCPs in Ref. Ghosh  _et al._ (2017) are
designed for suppressing leakage in the three-level system. This scheme Ghosh
_et al._ (2017) is required to simultaneously include two modulation
parameters, including the on-site potentials and the tunneling coupling. In
the current scheme, only one of the coupling strengths is sufficient to be
modulated to implement robust quantum control, while another coupling strength
remains unchanged during the evolution process.

##  IV Illustrative example

In this section, we mainly illustrate the applications of the current SMCPs
sequence for concrete physical models. This current sequence is suitable for
some systems where the phase is hardly to be modulated. For example, in
integrated photonic circuits Politi  _et al._ (2008), since the coupling
parameters between different waveguides are real valued, the traditional phase
modulation CPs is unapplicable. During the light transportation process from
one waveguide to another, one can alter the spatial distance between
waveguides, which is recognized as one kind of strength modulation between
different waveguides. And the amount of light coupling from one waveguide to
the other can be tuned by lithographically adjusting their width Politi  _et
al._ (2008). Thus, one can adopt strength modulation CPs to achieve robust
light transportation between different waveguides in integrated photonic
circuits Kyoseva  _et al._ (2019).

The current sequence is also feasible for the system where the strength
parameters are conveniently modulated. For instance, in the system consisting
of multiple quantum dots Reed  _et al._ (2016), there are two strength
parameters involved in the control of quantum dots at the same time. One is
the chemical potential of quantum dots, which is used for modulating the
detuning. The other is the tunneling barrier between adjacent quantum dots,
which alters the coupling of quantum dots. Note that both strength parameters
can be appropriately controlled by adjusting gate voltages, taking advantage
of short-range interaction and electrical controllability Reed  _et al._
(2016). Therefore, a proper value of the gate voltage would make two quantum
dots couple efficiently and produce a time-dependent exchange. So it is
convenient for employing the strength modulation in this system.

In the superconducting circuits system Xiang  _et al._ (2013), both the
strength modulation and the phase modulation can be used to achieve robust
quantum control by modulating the interaction of external fields. When
considering the weak dephasing in this system, the phase modulation method
Torosov  _et al._ (2021) is less robust against the phase noise than the
strength modulation method. Thus, systematic errors may not be effectively
eliminated due to the dephasing, and this problem can be elegantly
circumvented by the current SMCPs sequence. In the following, we choose this
system to illustrate how to achieve robust quantum control by our SMCPs
sequence.

The superconducting quantum interference devices (SQUIDs) are elementary
components in superconducting circuits Xia and Xie (2008); Yang  _et al._
(2003, 2004); Kang  _et al._ (2016); Guarcello  _et al._ (2020); Braumüller
_et al._ (2020), and the energy level of the SQUID qubits can be easily
modulated by changing the local bias fields Xiang  _et al._ (2013). As shown
in Fig. 4(a), the Hamiltonian of each SQUID qubit has the following form Yang
_et al._ (2003, 2004):

$H_{SQUID}=\frac{Q^{2}}{2C}+\frac{(\mathbf{\Phi}-\mathbf{\Phi}_{x})^{2}}{2L}-E_{J}\cos{(2\pi\frac{\mathbf{\Phi}}{\mathbf{\Phi}_{0}})},$
(49)

where $C$ is the junction capacitance, $L$ is the loop inductance, $Q$ is the
total charge on the capacitor, $\mathbf{\Phi}$ is the magnetic flux through
the loop, $\mathbf{\Phi}_{x}$ is the external flux applied to the loop,
$\mathbf{\Phi}_{0}=h/2e$ is the flux quantum, and
$E_{J}=I_{c}\mathbf{\Phi}_{0}/2\pi$ is the Josephson energy with $I_{c}$ being
the critical current of the junction.

Figure 4: (a) Schematic diagram of the SQUID qubit. (b) Four SQUID qubits
coupled to a single-mode cavity field. (c) The energy-level configuration of
every SQUID qubit.

Consider the system of four SQUID qubits coupled to a single-mode microwave
cavity field, as shown in Fig. 4(b). Every SQUID qubit has a three-level
structure with three states
$\\{|0\rangle_{s},|1\rangle_{s},|e\rangle_{s}\\}(s=1,2,3,4)$, as shown in Fig.
4(c), where $|0\rangle_{s}$ and $|1\rangle_{s}$ are two ground states and
$|e\rangle_{s}$ is the excited state. The ground state $|0\rangle_{s}$
resonantly couples to the excited state $|e\rangle_{s}$ by the cavity field
with the constant strength $g_{s}$. The classical field resonantly drives the
transition between the states $|1\rangle_{s}$ and $|e\rangle_{s}$ with the
coupling strength $\Omega_{s}$. The general forms of coupling strength
$\Omega_{s}$ and the coupling constant $g_{s}$ are given by Yang  _et al._
(2003, 2004)

$\displaystyle\Omega_{s}$ $\displaystyle=$
$\displaystyle\frac{1}{2L_{s}\hbar}\langle
1|\mathbf{\Phi}|e\rangle_{s}\int_{S_{s}}\mathbf{B}^{s}_{\mu
w}(\mathbf{r},t)\cdot d\mathbf{S},$ (50) $\displaystyle g_{s}$
$\displaystyle=$
$\displaystyle\frac{1}{L_{s}}\sqrt{\frac{\omega_{c}}{2\mu_{0}\hbar}}\langle
0|\mathbf{\Phi}|e\rangle_{s}\int_{S_{s}}\mathbf{B}^{s}_{c}(\mathbf{r})\cdot
d\mathbf{S},$ (51)

where $S_{s}$ is the surface surrounded by the ring of the $s\mathrm{th}$
SQUID, $L_{s}$ is the loop inductance of the $s\mathrm{th}$ SQUID,
$\omega_{c}$ is the angular frequency of the cavity, and $\mu_{0}$ is magnetic
permeability in a vacuum. $\mathbf{B}^{s}_{c}(\mathbf{r})$ is the magnetic
component of the normal cavity mode. $\mathbf{B}^{s}_{\mu w}(\mathbf{r},t)$ is
the magnetic component of the classical field.

From Eqs. (50) and (51), one finds that the coupling strength $\Omega_{s}$ and
the coupling constant $g_{s}$ can be controlled by the microwave magnetic
components $\mathbf{B}^{s}_{\mu w}(\mathbf{r},t)$ and the magnetic component
$\mathbf{B}^{s}_{c}(\mathbf{r})$ of the normal cavity mode, respectively. Due
to the spatial inhomogeneity of the magnetic fields, there may exist
uncertainty in $\Omega_{s}$ and $g_{s}$. On the other hand, the distorted
pulse shape or inexact evolution time could also produce uncertainty in the
pulse area. Consequently, the quantum operations would lose accuracy in this
system. To eliminate those uncertainties, we can employ the SMCPs scheme to
prepare the $W$ state with high fidelity, which is actually to control the
microwave magnetic components $\mathbf{B}^{s}_{\mu w}(\mathbf{r},t)$, in
practice. In the following, we adopt $\eta_{1}$ and $\eta_{2}$ to represent
the uncertainties from $\mathbf{B}^{s}_{\mu w}(\mathbf{r},t)$ and
$\mathbf{B}^{s}_{c}(\mathbf{r})$, respectively. The Hamiltonian of the whole
system in the interaction picture is $(\hbar=1)$ Kang  _et al._ (2016)

$\displaystyle H_{I}$ $\displaystyle=$ $\displaystyle H_{m}+H_{c},$ (52)
$\displaystyle H_{m}$ $\displaystyle=$
$\displaystyle\frac{1}{2}\sum^{4}_{s=1}\big{[}\Omega_{s}(1+\eta_{1})|e\rangle_{s}\langle
1|+\mathrm{H.c.}\big{]},$ (53) $\displaystyle H_{c}$ $\displaystyle=$
$\displaystyle\frac{1}{2}\sum^{4}_{s=1}\big{[}g_{s}(1+\eta_{2})|e\rangle_{s}\langle
0|\hat{a}+\mathrm{H.c.}\big{]},$ (54)

where $\hat{a}$ is the annihilation operator of the cavity field. $H_{c}$
($H_{m}$) describes the interaction between the cavity (classical) fields and
the SQUID qubits.

Note that the excited number is a conserved quantity in this system because
the excited number operator $\hat{N}_{e}=\sum_{s}(|e\rangle_{s}\langle
e|+|1\rangle_{s}\langle 1|)+\hat{a}^{{\dagger}}\hat{a}$ satisfies
$[\hat{N}_{e},H_{I}]=0$. Therefore, we can restrict the system dynamics in the
single excited subspace, namely,
$\langle\psi(t)|\hat{N}_{e}|\psi(t)\rangle=1$, where $|\psi(t)\rangle$ is the
evolution state of the system. The bases of the single excited subspace are

$\displaystyle|\psi_{1}\rangle$ $\displaystyle=$
$\displaystyle|0\rangle_{1}|0\rangle_{2}|0\rangle_{3}|1\rangle_{4}|0\rangle_{c},~{}~{}~{}|\psi_{2}\rangle=|0\rangle_{1}|0\rangle_{2}|0\rangle_{3}|e\rangle_{4}|0\rangle_{c},$
$\displaystyle|\psi_{3}\rangle$ $\displaystyle=$
$\displaystyle|0\rangle_{1}|0\rangle_{2}|0\rangle_{3}|0\rangle_{4}|1\rangle_{c},~{}~{}~{}|\psi_{4}\rangle=|e\rangle_{1}|0\rangle_{2}|0\rangle_{3}|0\rangle_{4}|0\rangle_{c},$
$\displaystyle|\psi_{5}\rangle$ $\displaystyle=$
$\displaystyle|0\rangle_{1}|e\rangle_{2}|0\rangle_{3}|0\rangle_{4}|0\rangle_{c},~{}~{}~{}|\psi_{6}\rangle=|0\rangle_{1}|0\rangle_{2}|e\rangle_{3}|0\rangle_{4}|0\rangle_{c},$
$\displaystyle|\psi_{7}\rangle$ $\displaystyle=$
$\displaystyle|1\rangle_{1}|0\rangle_{2}|0\rangle_{3}|0\rangle_{4}|0\rangle_{c},~{}~{}~{}|\psi_{8}\rangle=|0\rangle_{1}|1\rangle_{2}|0\rangle_{3}|0\rangle_{4}|0\rangle_{c},$
$\displaystyle|\psi_{9}\rangle$ $\displaystyle=$
$\displaystyle|0\rangle_{1}|0\rangle_{2}|1\rangle_{3}|0\rangle_{4}|0\rangle_{c}.$

For simplicity, we set $g_{1}=g_{2}=g_{3}=g$, $g_{4}=\sqrt{3}g$. This could be
realized by properly adjusting the parameters of the SQUIDs, such as $L_{s}$
or $S_{s}$. Then, we rewrite the Hamiltonian $H_{c}$ given by Eq. (54) in this
set of basis, i.e.,

$\begin{split}H_{c}=&g(1+\eta_{2})(|\psi_{4}\rangle+|\psi_{5}\rangle+|\psi_{6}\rangle)\langle\psi_{3}|\\\
&+\sqrt{3}g(1+\eta_{2})|\psi_{2}\rangle\langle\psi_{3}|+\mathrm{H.c.},\end{split}$
(55)

and the eigenstates of the Hamiltonian $H_{c}$ are

$\displaystyle|\phi_{1}\rangle$ $\displaystyle=$
$\displaystyle-\frac{1}{\sqrt{2}}\Big{[}|\psi_{2}\rangle\\!-\\!\frac{1}{\sqrt{3}}(|\psi_{4}\rangle\\!+\\!|\psi_{5}\rangle\\!+\\!|\psi_{6}\rangle)\Big{]},$
(56) $\displaystyle|\phi_{2}\rangle$ $\displaystyle=$
$\displaystyle\frac{1}{2}\Big{[}|\psi_{2}\rangle\\!+\\!\sqrt{2}|\psi_{3}\rangle\\!+\\!\frac{1}{\sqrt{3}}(|\psi_{4}\rangle\\!+\\!|\psi_{5}\rangle\\!+\\!|\psi_{6}\rangle)\Big{]},$
(57) $\displaystyle|\phi_{3}\rangle$ $\displaystyle=$
$\displaystyle\frac{1}{2}\Big{[}|\psi_{2}\rangle\\!-\\!\sqrt{2}|\psi_{3}\rangle\\!+\\!\frac{1}{\sqrt{3}}(|\psi_{4}\rangle\\!+\\!|\psi_{5}\rangle\\!+\\!|\psi_{6}\rangle)\Big{]},~{}~{}~{}~{}$
(58)

with the corresponding eigenvalues $E_{1}=0$, $E_{2}=\sqrt{6}g$, and
$E_{3}=-\sqrt{6}g$. Remarkably, $|\phi_{1}\rangle$ is the dark state of this
system. Next, we set $\Omega_{1}=\Omega_{2}=\Omega_{3}=\sqrt{2}\Omega_{a}$ and
$\Omega_{4}=\sqrt{2}\Omega_{b}$, and the Hamiltonian $H_{I}$ of the whole
system given by Eq. (52) can be rearranged into a new form,

$\displaystyle H^{\prime}_{I}$ $\displaystyle=$ $\displaystyle
H^{\prime}_{c}+H^{\prime}_{m},$ (59) $\displaystyle H^{\prime}_{c}$
$\displaystyle=$
$\displaystyle\sum^{3}_{i=1}E_{i}|\phi_{i}\rangle\langle\phi_{i}|,$ (60)
$\displaystyle H^{\prime}_{m}$ $\displaystyle=$
$\displaystyle\frac{\Omega_{a}}{\sqrt{2}}(1\\!+\\!\eta_{1})\big{(}|\phi_{2}\rangle\\!+\\!|\phi_{3}\rangle\big{)}\langle\psi_{1}|\\!+\\!\Omega_{a}(1\\!+\\!\eta_{1})|\phi_{1}\rangle\langle\psi_{1}|$
(61)
$\displaystyle\\!+\frac{\Omega_{b}}{\sqrt{2}}(1\\!+\\!\eta_{1})\big{(}|\phi_{2}\rangle\\!+\\!|\phi_{3}\rangle\big{)}\langle\psi_{f}|\\!+\\!\Omega_{b}(1\\!+\\!\eta_{1})|\phi_{1}\rangle\langle\psi_{f}|$
$\displaystyle\\!+\mathrm{H.c.},$

where
$|\psi_{f}\rangle=(|\psi_{7}\rangle+|\psi_{8}\rangle+|\psi_{9}\rangle)/\sqrt{3}$.
Finally, we perform the unitary transformation $U=e^{-iH^{\prime}_{c}t}$ on
this system, and the Hamiltonian $H^{\prime}_{I}$ given by Eq. (59) can be
rewritten as

$\displaystyle H^{\prime\prime}_{I}$ $\displaystyle=$
$\displaystyle\frac{\Omega_{a}}{\sqrt{2}}(1\\!+\\!\eta_{1})\Big{(}e^{i\sqrt{6}gt}|\phi_{2}\rangle+e^{-i\sqrt{6}gt}|\phi_{3}\rangle\Big{)}\langle\psi_{1}|$
(62)
$\displaystyle+\frac{\Omega_{b}}{\sqrt{2}}(1\\!+\\!\eta_{1})\Big{(}e^{i\sqrt{6}gt}|\phi_{2}\rangle+e^{-i\sqrt{6}gt}|\phi_{3}\rangle\Big{)}\langle\psi_{f}|$
$\displaystyle+\Omega_{a}(1\\!+\\!\eta_{1})|\phi_{1}\rangle\langle\psi_{1}|\\!+\\!\Omega_{b}(1\\!+\\!\eta_{1})|\phi_{1}\rangle\langle\psi_{f}|+\mathrm{H.c.}~{}~{}~{}~{}~{}~{}$

Figure 5: Schematic diagram of transition paths for the Hamiltonian
$H^{\prime\prime}_{I}$ given by Eq. (62).

Figure 5 shows the transition paths for the system governed by the Hamiltonian
$H^{\prime\prime}_{I}$. When the system parameters satisfy the condition
$\Omega_{a,b}\ll g$, the high-frequency oscillation terms $e^{\pm
i\sqrt{6}gt}$ can be safely ignored. Then, both the transition paths
$|\psi_{1}\rangle\\!\\!\leftrightarrow\\!\\!|\phi_{2}\rangle\\!\\!\leftrightarrow\\!\\!|\psi_{f}\rangle$
and
$|\psi_{1}\rangle\\!\\!\leftrightarrow\\!\\!|\phi_{3}\rangle\\!\\!\leftrightarrow\\!\\!|\psi_{f}\rangle$
are effectively suppressed, and only the transition path
$|\psi_{1}\rangle\\!\\!\leftrightarrow\\!\\!|\phi_{1}\rangle\\!\\!\leftrightarrow\\!\\!|\psi_{f}\rangle$
is feasible. As a result, the whole system is simplified into a $\Lambda$-type
physical model, and the effective Hamiltonian, in the basis
$\\{|\psi_{1}\rangle,|\psi_{f}\rangle,|\phi_{1}\rangle\\}$, reads Kang  _et
al._ (2016)

$\displaystyle
H_{e}=(1\\!+\\!\eta_{1})\left[\begin{array}[]{ccc}0&0&-\Omega_{a}\\\
0&0&\Omega_{b}\\\ -\Omega_{a}&\Omega_{b}&0\\\ \end{array}\right],$ (66)

where

$|\psi_{f}\rangle\\!=\\!\frac{1}{\sqrt{3}}(|1\rangle_{1}|0\rangle_{2}|0\rangle_{3}+|0\rangle_{1}|1\rangle_{2}|0\rangle_{3}+|0\rangle_{1}|0\rangle_{2}|1\rangle_{3})\otimes|0\rangle_{4}|0\rangle_{c}$

is the $W$ state. Therefore, the SMCPs scheme studied in Sec. III can be
directly employed to the effective Hamiltonian given by Eq. (66) to prepare
the $W$ state $|\psi_{f}\rangle$. Moreover, we can expect that the influence
caused by the uncertainty of the cavity field can be negligible during the
evolution process because $\eta_{2}$ does not appear in the effective
Hamiltonian given by Eq. (66).

Figure 6: The fidelity $F$ of the $W$ state vs two kinds of the uncertainties
$\eta_{1}$ and $\eta_{2}$ in the (a) two-pulses, (b) three-pulses, (c) four-
pulses, and (d) five-pulses sequence. $g/\Omega_{max}=7$. $\tau=10^{3}$. All
$\Theta_{n}$ come from Table 1.

In Figs. 6(a)–6(d), we plot the fidelity of the $W$ state as a function of the
uncertainties $\eta_{1}$ and $\eta_{2}$ with different pulses sequence, where
we set $\Omega_{max}=\max\\{\Omega_{na},\Omega_{nb}\\}$ and
$g/\Omega_{max}=7$. The numerical results verify that the uncertainty of the
coupling constant $g$ has a negligible effect on the final fidelity, since the
fidelity almost remains unchanged with the increase of $\eta_{2}$. Note that
the fidelity obviously changes when $\eta_{2}<-0.5$. This is because the
condition $g\gg\Omega_{a,b}$ is no longer satisfied very well. As a result,
the effective Hamiltonian given by Eq. (66) is invalid in this regime. On the
other hand, the fidelity is also robust against the uncertainty of the
coupling strength of the classical fields since the fidelity still keeps a
relatively high value even though $\eta_{1}=\pm 0.2$. Thus, we can achieve the
$W$ state in a robust way even when there are existing uncertainties in the
superconducting circuits.

Next, we turn to study the influence of waveform deformation on the fidelity
of the $W$ state. We take the three-pulses sequence as an example, and plot in
Fig. 7 the infidelity as a function of the pulse area error and the waveform
deformation. An inspection of Fig. 7 demonstrates that the scheme is
insensitive to waveform deformation, even in the case of severe deformation,
because the $W$ state still maintains a high fidelity in the region $\tau>3$.
Hence, our scheme can resist the imperfection of input waveform and thus give
more flexibility in terms of input waveform conditions.

Figure 7: The infidelity $1-F$ of the $W$ state vs the uncertainty $\eta_{1}$
and the deformation parameter $\tau$ in the three-pulses scheme, where
$g/\Omega_{max}=15$. The green dotted curve, the yellow dashed curve, the pink
dot-dashed curve, and the white solid curve correspond to the infidelity
$1-F=10^{-4}$, $10^{-3}$, $10^{-2}$, and $10^{-1}$. All $\Theta_{n}$ come from
Table 1.

Until now, we have not investigated the influence of decoherence on the
fidelity of the $W$ state. When considering the decoherence induced by the
cavity decay, the spontaneous emission, and the dephasing, the evolution of
this system can be governed by the Lindblad master equation,

$\dot{\rho}=i[\rho,H_{I}]+\mathcal{L}(\sqrt{\kappa}\hat{a})\rho+\sum_{l=0}^{1}\sum^{4}_{s=1}\Big{[}\mathcal{L}(\hat{\sigma}^{s}_{le})\rho+\mathcal{L}(\hat{\sigma}^{s}_{ee,ll})\rho\Big{]},\\\
$ (67)

where $\hat{\sigma}^{s}_{le}=\sqrt{\gamma_{ls}}|l\rangle_{s}\langle e|$,
$\hat{\sigma}^{s}_{ee,ll}=\sqrt{{\gamma^{\phi}_{ls}}/{2}}(|e\rangle_{s}\langle
e|-|l\rangle_{s}\langle l|)$, and the general form of the superoperator is

$\mathcal{L}(\hat{o})\rho=\hat{o}\rho\hat{o}^{{\dagger}}-\frac{1}{2}(\hat{o}^{{\dagger}}\hat{o}\rho+\rho\hat{o}^{{\dagger}}\hat{o}).$
(68)

Here, $\hat{o}$ denotes the standard Lindblad operator, $\kappa$ is the decay
rate of the cavity field, and $\gamma_{ls}$ ($\gamma^{\phi}_{ls}$) represent
the dissipation (dephasing) rate from the excited state $|e\rangle$ to the
ground state $|l\rangle$.

Figure 8: The fidelity of the $W$ state vs (a) $\gamma$/g and $\kappa$/g, (b)
$\gamma$/g and $\gamma^{\phi}$/g, and (c) $\gamma^{\phi}$/g and $\kappa$/g in
the three-pulses scheme. Here, the green dotted line, the yellow dashed line,
the pink dot-dashed line, and the white solid line correspond to the final
fidelity $F=0.99$, $0.97$, $0.95$, and $0.90$. (d) The population evolution of
the states $|\phi_{k}\rangle$, $|\psi_{1}\rangle$, and $|\psi_{f}\rangle$ in
the three pulses scheme, where the initial state is $|\psi_{1}\rangle$.
$g/\Omega_{max}=7$ and $\tau=10^{3}$. All $\Theta_{n}$ come from Table 1.

Figures 8(a)–8(c) demonstrate the relationship between the fidelity of the $W$
state and the decoherence parameters $\\{\gamma,\gamma^{\phi},\kappa\\}$,
where we set $\gamma_{ls}=\gamma$ and $\gamma^{\phi}_{ls}=\gamma^{\phi}$ for
simplicity. According to Fig. 8(a) and Fig. 8(c), the fidelity of the $W$
state is hardly affected by the cavity decay. Physically, it is due to the
fact that the single-photon state of the cavity field is almost decoupling
during the system evolution. We can observe from Fig. 8(a) and Fig. 8(b) that
the spontaneous emission of the SQUID qubits has a slight influence on the
fidelity of the $W$ state. The reason is as follows. Although
$|\phi_{k}\rangle$ ($k=1,2,3$) contain the excited states of the SQUID qubits,
as given by Eqs. (56)–(58), they are almost negligible in the evolution
process when satisfying the condition $\Omega_{a,b}\ll g$. This can be
verified by Fig. 8(d), which demonstrates that the population of the state
$|\phi_{k}\rangle$ ($k=1,2,3$) is suppressed within an extremely small range
during the system evolution. As a result, the current scheme can maintain a
high fidelity even when the dissipation rate of the excited state is large.

It is inevitable that the increase of the CPs sequence would prolong the total
interaction time and impact the final fidelity under the decoherence
environment. In the following, we study this issue. Figure 9(a) shows the
results obtained by various pulse lengths under different decoherence
environments. In the ideal situation, i.e., $\gamma/g=0$, the final fidelity
almost remains unchanged as the pulse number increases. When considering the
decoherence caused by the atomic spontaneous emission and dephasing, e.g.,
$\gamma/g=2\times 10^{-5}$, the final fidelity decreases slightly over an
increasing number of pulses. The reason is that a long pulse sequence takes a
long interaction time so that the decoherence becomes the main factor in
reducing the final fidelity. Hence, to obtain a relatively high fidelity in
the decoherence environment, we cannot choose an overlong pulse sequence.

Figure 9:  (a) The final fidelity of the $W$ state and (b) the benchmark width
$W_{b}$ vs different pulse number $N$. Here we set the benchmark value
$F_{b}=0.98$ and $\gamma_{p}=\gamma$. All $\Theta_{n}$ come from Table 1.
Other parameters are $\kappa/g=0.01$, $g/\Omega_{max}=20$, and $\tau=10^{3}$.

To quantify the robustness under the decoherence environment, we define a
benchmark width

$\displaystyle W_{b}$ $\displaystyle=$
$\displaystyle\eta_{1,max}-\eta_{1,min},$ (69)

where $\eta_{1,max(min)}$ is the positive (negative) maximum uncertainty in
the coupling strength $\Omega_{s}$ to make the final fidelity reach a
benchmark value $F_{b}$. Remarkably, this benchmark width can characterize the
robustness against uncertainty provided by different pulse sequences. Figure
9(b) reveals the benchmark widths achieved by different pulse sequences under
different decoherence environments. On the one hand, we can find that the too
short pulse sequence cannot obtain a remarkable benchmark width (e.g.,
$W_{b}=0.3828$ for $N=2$ and $\gamma/g=2\times 10^{-5}$). On the other hand,
the benchmark width becomes wider along with the pulse number increasing, and
the increment of $W_{b}$ becomes gradually limited when the pulse number is
large. In other words, when the system suffers from a serious decoherence, we
cannot obtain considerable robustness for an overlong sequence (e.g.,
$W_{b}=0.6899$ for $N=6$ and $\gamma/g=5\times 10^{-5}$). From this point of
view, to obtain excellent robustness, both a too short pulse sequence and an
overlong pulse sequence should not be considered preferentially.

Therefore, we need to choose a moderate pulse length to seek an optimal
tradeoff between the final fidelity and robustness (fault tolerance). In
recent experiments Xiang  _et al._ (2013); Stern  _et al._ (2014); Leek  _et
al._ (2010); Yan  _et al._ (2016); Ma  _et al._ (2020b); Xu  _et al._ (2020),
the following parameters are feasible: $\gamma=\gamma^{\phi}\sim
10\mathrm{kHz}$ and $\kappa\sim 625\mathrm{kHz}$. The coupling constant $g$
could reach $600\mathrm{MHz}$ Niemczyk  _et al._ (2010), and thus
$\gamma/g\sim 1.6667\times 10^{-5}$. According to these parameters, the final
fidelity of the $W$ state can reach 0.9955 by the four-pulses sequence, while
the benchmark width is $W_{b}=0.7307$ for $F_{b}=0.98$. As a result, our
scheme still retains good performance in a decoherence environment.

##  V Conclusion

In summary, we have put forward the SMCPs scheme to achieve robust population
inversion in a three-level system. The transition probability is directly
derived from a total propagator, and we mainly focus on reducing the influence
caused by the pulse area error on the transition probability. By means of the
Taylor expansion, we nullify the first few error terms in the transition
probability through designing the values of the coupling strengths. We have
given the analytical expressions of the coupling strengths in the two-pulses
sequence. For the three- and more-than-three-pulses (up to seven-pulses)
sequence, we present a set of optimal parameters by numerical methods. These
results indicate that the robust manner against the pulse area error can be
much better as the pulse number increases. Our scheme provides a feasibility
extension of a longer sequence. Meanwhile, this scheme provides a flexible
choice to design pulse length in a desirable manner for achieving high
robustness against errors and an appropriate total operating time. Moreover,
the pulse sequence interfered by waveform deformation is still valid to
achieve high transition probability. This property is very useful for the
resistance of pulse jump delay. As a result, our SMCPs scheme is robust
against the pulse area error and waveform deformations.

Then, we have extended this SMCPs scheme to prepare the $W$ state with high
fidelity in the superconducting circuits, where four SQUID qubits are coupled
with the cavity field and classical fields. The SMCPs that we design can
effectively compensate for the adverse impact caused by waveform deformations
and uncertainties in coupling strengths. Furthermore, we also discuss the
feasibility and the robustness of this scheme in a decoherence environment.
The simulation results show that our scheme can obtain a good performance
against the impact of the cavity decay and the spontaneous emission, and we
need to select a moderate pulse length to achieve robust and accurate quantum
control. We hope that this scheme will provide robust control for more
physical models in the future.

###### Acknowledgements.

This work is supported by the National Natural Science Foundation of China
under Grants No. 11874114 and No. 11805036, the Natural Science Foundation of
Fujian Province under Grant No. 2021J01575, the Natural Science Funds for
Distinguished Young Scholar of Fujian Province under Grant No. 2020J06011, and
the Project from Fuzhou University under Grant No. JG202001-2.

##  Appendix A Derivation of the propagator for the $N$ pulses sequence in the
two-photon resonance regime

In this appendix, we deduce the expression of the propagator without errors
when the three-level system works in the two-photon resonance regime, where
all parameters are time-independent. The general form of a three-level system
Hamiltonian for the $n\mathrm{th}$ pulse is given by

$H_{n}=\Delta_{n}|e\rangle\langle
e|+\frac{1}{2}\Big{(}\Omega_{n}e^{i\phi_{n}}|g\rangle\langle
e|+\lambda_{n}e^{i\varphi_{n}}|f\rangle\langle e|+\mathrm{H.c.}\Big{)},$ (70)

where $\Omega_{n}$ and $\phi_{n}$ ($\lambda_{n}$ and $\varphi_{n}$) are the
coupling strength and the phase of the transition
$|g\rangle\leftrightarrow|e\rangle$ ($|f\rangle\leftrightarrow|e\rangle$), and
$\Delta_{n}$ is the single-photon detuning. In this case, the propagator
becomes $U_{n}(\Theta_{n})=\exp\big{(}{-iH_{n}T_{n}}\big{)}$, where $T_{n}$ is
the pulse duration and $\Theta_{n}=\arctan(\Omega_{n}/\lambda_{n})$. By
selecting the pulse duration
$T_{n}=2\pi/\sqrt{\Omega_{n}^{2}+\lambda_{n}^{2}+\Delta_{n}^{2}}$, two ground
states are decoupled from the excited state $|e\rangle$ so that we only focus
on the subspace $\\{|g\rangle,|f\rangle\\}$. Thus, the matrix form of the
propagator $U_{n}(\Theta_{n})$ in the basis $\\{|g\rangle,|f\rangle\\}$ reads

$\displaystyle U_{n}(\Theta_{n})$ $\displaystyle=$
$\displaystyle\cos{x_{n}}\left[\begin{array}[]{cc}\cos{2\Theta_{n}}+i\tan{x_{n}}&-e^{i\Psi_{n}}\sin{2\Theta_{n}}\\\
-e^{-i\Psi_{n}}\sin{2\Theta_{n}}&-\cos{2\Theta_{n}}+i\tan{x_{n}}\\\
\end{array}\right],$ (73)

where the global phase is neglected, and

$\displaystyle
x_{n}=\frac{\pi\Delta_{n}}{2\sqrt{1+\Delta^{2}_{n}}},~{}~{}~{}~{}\Psi_{n}$
$\displaystyle=$ $\displaystyle\phi_{n}-\varphi_{n}.$

Note that the total propagator for the composite pulses sequence could be
obtained by nesting the propagator $U_{n}(\Theta_{n})$ in order. We first
start from the two-pulses sequence, and the propagator comes into

$\displaystyle U^{(2)}$ $\displaystyle=$ $\displaystyle
U_{2}(\Theta_{2})U_{1}(\Theta_{1})$ $\displaystyle=$
$\displaystyle\cos{x_{1}}\cos{x_{2}}\left[\begin{array}[]{cc}\cos{2\Theta_{2}}+i\tan{x_{2}}&-e^{i\Psi}\sin{2\Theta_{2}}\\\
-e^{-i\Psi}\sin{2\Theta_{2}}&-\cos{2\Theta_{2}}+i\tan{x_{2}}\\\
\end{array}\right]\left[\begin{array}[]{cc}\cos{2\Theta_{1}}+i\tan{x_{1}}&-e^{i\Psi}\sin{2\Theta_{1}}\\\
-e^{-i\Psi}\sin{2\Theta_{1}}&-\cos{2\Theta_{1}}+i\tan{x_{1}}\\\
\end{array}\right]$ $\displaystyle=$
$\displaystyle\cos{x_{1}}\cos{x_{2}}\left[\begin{array}[]{cc}\cos{2(\Theta_{1}-\Theta_{2})}+g_{2}&\big{[}-\sin{2(\Theta_{1}-\Theta_{2})}+h_{2}\big{]}e^{i\Psi}\\\\[1.29167pt]
\big{[}\sin{2(\Theta_{1}-\Theta_{2})}-h_{2}^{*}\big{]}e^{-i\Psi}&-\cos{2(\Theta_{1}-\Theta_{2})}+g_{2}^{*}\\\\[1.29167pt]
\end{array}\right],$

with

$\displaystyle g_{2}$ $\displaystyle=$
$\displaystyle-\tan{x_{1}}\tan{x_{2}}+i\left(\tan{x_{1}}\cos{2\Theta_{2}}+\tan{x_{2}}\cos{2\Theta_{1}}\right),$
$\displaystyle h_{2}$ $\displaystyle=$
$\displaystyle-i\left(\tan{x_{1}}\sin{2\Theta_{2}}+\tan{x_{2}}\sin{2\Theta_{1}}\right),$

where $g_{2}$ and $h_{2}$ are defined as functions associated with the
parameters $\\{x_{1},x_{2},\Theta_{1},\Theta_{2}\\}$. Similarly, the
propagator of the three-pulses sequence is

$\displaystyle U^{(3)}$ $\displaystyle=$ $\displaystyle
U_{3}(\Theta_{2})U_{2}(\Theta_{2})U_{1}(\Theta_{1})$ $\displaystyle=$
$\displaystyle\left(\prod_{n=1}^{3}\cos{x_{n}}\right)\left[\begin{array}[]{cc}\cos{2(\Theta_{1}-\Theta_{2}+\Theta_{3})+g_{3}}&\big{[}-\sin{2(\Theta_{1}-\Theta_{2}+\Theta_{3})}+h_{3}\big{]}e^{i\Psi}\\\\[1.29167pt]
\big{[}-\sin{2(\Theta_{1}-\Theta_{2}+\Theta_{3})}+h_{3}^{*}\big{]}e^{-i\Psi}&-\cos{2(\Theta_{1}-\Theta_{2}+\Theta_{3})}-g_{3}^{*}\end{array}\right],$

where

$\displaystyle g_{3}$ $\displaystyle=$
$\displaystyle-(\tan{x_{1}}\tan{x_{2}}\cos{2\Theta_{3}}+\tan{x_{1}}\tan{x_{3}}\cos{2\Theta_{2}}+\tan{x_{2}}\tan{x_{3}}\cos{2\Theta_{1}})$
$\displaystyle+i\big{[}\tan{x_{1}}\cos{2(\Theta_{2}-\Theta_{3})}+\tan{x_{2}}\cos{2(\Theta_{1}-\Theta_{3})}+\tan{x_{3}}\cos{2(\Theta_{1}-\Theta_{2})}-\tan{x_{1}}\tan{x_{2}}\tan{x_{3}}\big{]},$
$\displaystyle h_{3}$ $\displaystyle=$
$\displaystyle(\tan{x_{1}}\tan{x_{2}}\sin{2\Theta_{3}}+\tan{x_{1}}\tan{x_{3}}\sin{2\Theta_{2}}+\tan{x_{2}}\tan{x_{3}}\sin{2\Theta_{1}})$
$\displaystyle-i\big{[}\tan{x_{1}}\sin{2(\Theta_{2}-\Theta_{3})}+\tan{x_{2}}\sin{2(\Theta_{1}-\Theta_{3})}+\tan{x_{3}}\sin{2(\Theta_{1}-\Theta_{2})}\big{]}.$

Following the above rules, the propagator of the $N$-pulses sequence could be
obtained and reads

$\displaystyle U^{(N)}$ $\displaystyle=$ $\displaystyle
U_{N}(\Theta_{N})U_{N-1}(\Theta_{N-1})\cdots U_{1}(\Theta_{1})$ (85)
$\displaystyle=$
$\displaystyle\left(\prod_{n=1}^{N}\cos{x_{n}}\right)\left[\begin{array}[]{cc}\cos{2\Theta_{N}}+i\tan{x_{N}}&-e^{i\Psi}\sin{2\Theta_{N}}\\\
-e^{-i\Psi}\sin{2\Theta_{N}}&-\cos{2\Theta_{N}}+i\tan{x_{N}}\\\
\end{array}\right]\cdots\left[\begin{array}[]{cc}\cos{2\Theta_{1}}+i\tan{x_{1}}&-e^{i\Psi}\sin{2\Theta_{1}}\\\
-e^{-i\Psi}\sin{2\Theta_{1}}&-\cos{2\Theta_{1}}+i\tan{x_{1}}\\\
\end{array}\right]$ $\displaystyle=$
$\displaystyle\left(\prod_{n=1}^{N}\cos{x_{n}}\right)\left[\begin{array}[]{cc}\cos{(2\sum^{N}_{n=1}(-1)^{n+1}\Theta_{n})}+g_{N}&\left[-\sin{(2\sum^{N}_{n=1}(-1)^{n+1}\Theta_{n})}+h_{N}\right]e^{i\Psi}\\\\[1.29167pt]
(-1)^{N}\left[\sin{(2\sum^{N}_{n=1}(-1)^{n+1}\Theta_{n})}-h_{N}^{*}\right]e^{-i\Psi}&(-1)^{N}\left[\cos{(2\sum^{N}_{n=1}(-1)^{n+1}\Theta_{n})}+g_{N}^{*}\right]\end{array}\right]\\!,~{}~{}~{}$
(88)

where

$\displaystyle g_{N}$ $\displaystyle=$
$\displaystyle\sum^{N-1}_{a=1}\left[\prod^{N}_{1\leq j_{1}\leq\cdots\leq
j_{a}\leq
N}\left(i\tan{x_{j_{l}}}\right)\right]\cos{\left[\sum^{N}_{\begin{subarray}{c}1\leq
j^{\prime}_{1}\leq\cdots\leq j^{\prime}_{N-a}\leq N\\\ j_{l}\neq
j^{\prime}_{m}\end{subarray}}(-1)^{(m+1)}\left(2\Theta_{j^{\prime}_{m}}\right)\right]}+\prod^{N}_{n=1}\left(i\tan{x_{n}}\right),$
$\displaystyle h_{N}$ $\displaystyle=$
$\displaystyle-\sum^{N-1}_{a=1}\left[\prod^{N}_{1\leq j_{1}\leq\cdots\leq
j_{a}\leq
N}\left(i\tan{x_{j_{l}}}\right)\right]\sin{\left[\sum^{N}_{\begin{subarray}{c}1\leq
j^{\prime}_{1}\leq\cdots\leq j^{\prime}_{N-a}\leq N\\\ j_{l}\neq
j^{\prime}_{m}\end{subarray}}(-1)^{(m+1)}\left(2\Theta_{j^{\prime}_{m}}\right)\right]}.$

We observe from Eq. (85) that the transition probability of the target state
$|f\rangle$ is concurrently controlled by the coupling strengths and detuning.
When the system is always in the resonance regime ($\Delta_{n}=0$), it is
easily found that $\cos{x_{n}}=1$, $g_{N}=0$, and $h_{N}=0$. Then, the final
propagator in the absence of errors becomes very concise, which reads

$\displaystyle
U^{(N)}=\left[\begin{array}[]{cc}\cos{\left[2\sum^{N}_{n=1}(-1)^{n+1}\Theta_{n}\right]}&-\sin{\left[2\sum^{N}_{n=1}(-1)^{n+1}\Theta_{n}\right]}e^{i\Psi}\\\\[1.29167pt]
(-1)^{N}\sin{\left[2\sum^{N}_{n=1}(-1)^{n+1}\Theta_{n}\right]}e^{-i\Psi}&(-1)^{N}\cos{\left[2\sum^{N}_{n=1}(-1)^{n+1}\Theta_{n}\right]}\\\\[1.29167pt]
\end{array}\right].$

As a result, the expression of the transition probability of the target state
$|f\rangle$ is

$\displaystyle P_{N}=|U_{21}|^{2}$ $\displaystyle=$
$\displaystyle\sin^{2}{\left[2\sum^{N}_{n=1}(-1)^{n+1}\Theta_{n}\right]},$

which is Eq. (50) in the main text. In order to achieve complete population
inversion, the ratio of coupling strength for each pulse, obtained by solving
the equation $P_{N}=1$, must be satisfied, and the solution is

$\displaystyle\sum^{N}_{n=1}(-1)^{n+1}\Theta_{n}=\frac{\pi}{4}+\frac{m\pi}{2},$
(90)

where $m$ is an arbitrary integer. This is the first condition that
$\Theta_{n}$ should be satisfied.

##  Appendix B Detailed derivation process in the two pulses sequence

In this appendix, we present the detailed derivation process of the coupling
strengths in the two-pulses sequence. When $\Delta=0$, according to Eq. (6) in
the main text, the total propagator in the basis
$\\{|g\rangle,|f\rangle,|e\rangle\\}$ is given by

$\displaystyle U^{(2)}$ $\displaystyle=$ $\displaystyle
U_{2}(A_{2},\Theta_{2},\phi_{2},\varphi_{2})U_{1}(A_{1},\Theta_{1},\phi_{1},\varphi_{1})$
(98) $\displaystyle=$
$\displaystyle\left[\begin{array}[]{ccc}\displaystyle\cos^{2}{\Theta_{2}}+\sin^{2}{\Theta_{2}}\cos{\frac{A_{2}}{2}}&\displaystyle-e^{i\Psi_{2}}\sin^{2}{\frac{A_{2}}{4}}\sin{2\Theta_{2}}&\displaystyle-
ie^{i\phi_{2}}\sin{\Theta_{2}}\sin{\frac{A_{2}}{2}}\\\\[6.45831pt]
\displaystyle-e^{-i\Psi_{2}}\sin^{2}{\frac{A_{2}}{4}}\sin{2\Theta_{2}}&\displaystyle\sin^{2}{\Theta_{2}}+\cos^{2}{\Theta_{2}}\cos{\frac{A_{2}}{2}}&\displaystyle-
ie^{i\varphi_{2}}\cos{\Theta_{2}}\sin{\frac{A_{2}}{2}}\\\\[6.45831pt]
\displaystyle-
ie^{-i\phi_{2}}\sin{\Theta_{2}}\sin{\frac{A_{2}}{2}}&\displaystyle-
ie^{-i\varphi_{2}}\cos{\Theta_{2}}\sin{\frac{A_{2}}{2}}&\displaystyle\cos{\frac{A_{2}}{2}}\\\\[6.45831pt]
\end{array}\right]$
$\displaystyle\times\left[\begin{array}[]{ccc}\displaystyle\cos^{2}{\Theta_{1}}+\sin^{2}{\Theta_{1}}\cos{\frac{A_{1}}{2}}&\displaystyle-e^{i\Psi_{1}}\sin^{2}{\frac{A_{1}}{4}}\sin{2\Theta_{1}}&\displaystyle-
ie^{i\phi_{1}}\sin{\Theta_{1}}\sin{\frac{A_{1}}{2}}\\\\[6.45831pt]
\displaystyle-e^{-i\Psi_{1}}\sin^{2}{\frac{A_{1}}{4}}\sin{2\Theta_{1}}&\displaystyle\sin^{2}{\Theta_{1}}+\cos^{2}{\Theta_{1}}\cos{\frac{A_{1}}{2}}&\displaystyle-
ie^{i\varphi_{1}}\cos{\Theta_{1}}\sin{\frac{A_{1}}{2}}\\\\[6.45831pt]
\displaystyle-
ie^{-i\phi_{1}}\sin{\Theta_{1}}\sin{\frac{A_{1}}{2}}&\displaystyle-
ie^{-i\varphi_{1}}\cos{\Theta_{1}}\sin{\frac{A_{1}}{2}}&\displaystyle\cos{\frac{A_{1}}{2}}\\\\[6.45831pt]
\end{array}\right]$ $\displaystyle=$
$\displaystyle\left[\begin{array}[]{ccc}U^{(2)}_{11}&U^{(2)}_{12}&U^{(2)}_{13}\\\\[3.44444pt]
U^{(2)}_{21}&U^{(2)}_{22}&U^{(2)}_{23}\\\\[3.44444pt]
U^{(2)}_{31}&U^{(2)}_{32}&U^{(2)}_{33}\\\\[3.44444pt] \end{array}\right].$
(102)

When the system is in the ground state $|g\rangle$, the time evolution of the
system state reads
$|\psi(t)\rangle=U^{(2)}_{11}|g\rangle+U^{(2)}_{21}|f\rangle+U^{(2)}_{31}|e\rangle$.
Consider the error $\epsilon$ occurring in the pulse area:
$\tilde{A}_{n}=A_{n}(1+\epsilon)$, $n=1,2$. The transition amplitude of the
state $|f\rangle$ can be rewritten as

$\displaystyle U^{(2)}_{21}(\epsilon)$ $\displaystyle=$
$\displaystyle-e^{-i\Psi_{1}}\sin^{2}{\frac{\tilde{A}_{1}}{4}}\sin{2\Theta_{1}}(\sin^{2}{\Theta_{2}}+\cos{\frac{\tilde{A}_{2}}{2}}\cos^{2}{\Theta_{2}})-e^{-i\Psi_{2}}\sin^{2}{\frac{\tilde{A}_{2}}{4}}\sin{2\Theta_{2}}(\cos^{2}{\Theta_{1}}+\cos{\frac{\tilde{A}_{1}}{2}}\sin^{2}{\Theta_{1}})$
(104)
$\displaystyle-e^{-i(\phi_{1}-\varphi_{2})}\sin{\Theta_{1}}\cos{\Theta_{2}}\sin{\frac{\tilde{A}_{1}}{2}}\sin{\frac{\tilde{A}_{2}}{2}},$

where $\Psi_{n}=\phi_{n}-\varphi_{n}$.

The transition probability of the state $|f\rangle$ is equal to the modular
squaring of the probability amplitude, i.e., $|U^{(2)}_{21}|^{2}$. Hence, the
transition probability carries the erroneous pulse area $\tilde{A}_{n}$ as

$\displaystyle P_{2}(\epsilon)$ $\displaystyle=$
$\displaystyle\sin^{2}\\!\Theta_{1}\\!\cos^{2}\\!\Theta_{2}\\!\sin^{2}\\!\frac{\tilde{A}_{1}}{2}\\!\sin^{2}\\!\frac{\tilde{A}_{2}}{2}\\!+\\!8\big{[}\cos^{2}\\!\Theta_{2}\\!\sin\\!\Theta_{1}\\!\sin\\!\Theta_{2}\\!\sin^{3}\\!\frac{\tilde{A}_{2}}{4}\\!\sin\\!\frac{\tilde{A}_{1}}{2}\\!\cos\\!\frac{\tilde{A}_{2}}{4}\\!\cos(\phi_{1}\\!-\\!\phi_{2})\cos\\!2\Theta_{1}\sin^{2}\\!\frac{\tilde{A}_{1}}{4}\\!+\\!\cos^{2}\\!\frac{\tilde{A}_{1}}{4}$
(105)
$\displaystyle+\sin^{2}\\!\Theta_{1}\\!\cos\\!\Theta_{1}\\!\cos\\!\Theta_{2}\\!\sin^{3}\\!\frac{\tilde{A}_{1}}{4}\cos\\!\frac{\tilde{A}_{1}}{4}\sin\\!\frac{\tilde{A}_{2}}{2}\\!\cos(\varphi_{1}\\!-\\!\varphi_{2})\cos^{2}\Theta_{2}\\!\cos\\!\frac{\tilde{A}_{2}}{2}\\!+\\!\sin^{2}\\!\Theta_{2}\big{]}\\!$
$\displaystyle+\\!\sin^{2}\\!2\Theta_{1}\\!\sin^{4}\frac{\tilde{A}_{1}}{4}\Big{[}\cos^{2}\Theta_{2}\cos\frac{\tilde{A}_{2}}{2}\\!+\\!\sin^{2}\\!\Theta_{2}\Big{]}^{2}+\sin^{2}\\!2\Theta_{2}\sin^{4}\frac{\tilde{A}_{2}}{4}\Big{[}\sin^{2}\Theta_{1}\cos\frac{\tilde{A}_{1}}{2}\\!+\\!\cos^{2}\Theta_{1}\Big{]}^{2}\\!$
$\displaystyle+2\sin 2\Theta_{1}\sin
2\Theta_{2}\sin^{2}\\!\frac{\tilde{A}_{1}}{4}\\!\sin^{2}\frac{\tilde{A}_{2}}{4}\\!\Big{[}\cos
2\Theta_{1}\sin^{2}\frac{\tilde{A}_{1}}{4}\\!+\\!\cos^{2}\frac{\tilde{A}_{1}}{4}\Big{]}\Big{[}\cos^{2}\Theta_{2}\cos\frac{\tilde{A}_{2}}{2}+\sin^{2}\Theta_{2}\Big{]}\\!\cos\\!\Gamma_{1},$

where $\Gamma_{1}=\Psi_{1}-\Psi_{2}$. After substituting the pulse area
$A_{1}=A_{2}=2\pi$ into Eq. (105), we have

$\displaystyle P_{2}(\epsilon)$ $\displaystyle=$
$\displaystyle\sin^{2}\Theta_{1}\cos^{2}\Theta_{2}\sin^{4}\pi\epsilon\\!+\\!\cos^{4}\frac{\pi\epsilon}{2}\big{[}\sin^{2}2\Theta_{1}(\sin^{2}\Theta_{2}\\!-\\!\cos^{2}\Theta_{2}\cos\pi\epsilon)^{2}\\!+\\!\sin^{2}2\Theta_{2}(\cos^{2}\Theta_{1}\\!-\\!\sin^{2}\Theta_{1}\cos\pi\epsilon)^{2}\big{]}$
(106) $\displaystyle+2\sin 2\Theta_{1}\sin
2\Theta_{2}\cos^{4}\frac{\pi\epsilon}{2}(\cos
2\Theta_{1}\cos^{2}\frac{\pi\epsilon}{2}+\sin^{2}\frac{\pi\epsilon}{2})(\sin^{2}\Theta_{2}-\cos^{2}\Theta_{2}\cos\pi\epsilon)\cos\Gamma_{1}$
$\displaystyle+4\sin^{2}\pi\epsilon\cos^{2}\\!\frac{\pi\epsilon}{2}\big{[}\sin\Theta_{1}\sin\Theta_{2}\cos^{2}\Theta_{2}\cos(\phi_{1}-\phi_{2})(\cos
2\Theta_{1}\cos^{2}\frac{\pi\epsilon}{2}+\sin^{2}\frac{\pi\epsilon}{2})$
$\displaystyle+\cos\Theta_{1}\cos\Theta_{2}\sin^{2}\Theta_{1}\cos(\varphi_{1}-\varphi_{2})(\sin^{2}\Theta_{2}-\cos^{2}\Theta_{2}\cos\pi\epsilon)\big{]}.$

Then, by the Taylor expansion, $P_{2}(\epsilon)$ could be written as the
following form (up to the second-order term):

$P_{2}(\epsilon)=\alpha_{2,0}+\alpha_{2,2}\epsilon^{2}+O(\epsilon^{4}).$ (107)

Here, the zero-order is the precise transition probability, and the goal of
designing the composite pulses sequence is to eliminate as many of the high-
order coefficients as possible. Note that all odd-order coefficients are
eliminated due to the choice of $A_{1}=A_{2}=2\pi$. The expressions of the
zeroth-order and the second-order coefficients are

$\displaystyle\alpha_{2,0}$ $\displaystyle=$
$\displaystyle\frac{1}{2}\big{(}1-\cos{4\Theta_{1}}\cos{4\Theta_{2}}-\cos{\Gamma_{1}}\sin{4\Theta_{1}}\sin{4\Theta_{2}}\big{)},$
(108) $\displaystyle\alpha_{2,2}$ $\displaystyle=$
$\displaystyle\frac{\pi^{2}}{4}\big{[}\sin\\!4\Theta_{1}(\sin\\!2\Theta_{2}\\!+\\!2\\!\sin\\!4\Theta_{2})\\!-\\!\sin\\!2\Theta_{1}\\!\sin\\!4\Theta_{2}\big{]}\\!\cos\\!\Gamma_{1}\\!+\\!\frac{\pi^{2}}{2}(\cos\\!4\Theta_{1}\\!\cos\\!4\Theta_{2}\\!-\\!\sin^{2}\\!2\Theta_{1}\\!\cos\\!2\Theta_{2}\\!+\\!\cos\\!2\Theta_{1}\\!\sin^{2}\\!2\Theta_{2}\\!-\\!1)$
(110) $\displaystyle+2\pi^{2}\big{[}\sin\Theta_{2}\cos^{2}\Theta_{2}(\sin
3\Theta_{1}-\sin\Theta_{1})\cos(\phi_{1}-\phi_{2})-\sin^{2}\Theta_{1}\cos\Theta_{1}(\cos\Theta_{2}+\cos
3\Theta_{2})\cos(\varphi_{1}-\varphi_{2})\big{]}$

In the two-pulses sequence, there are only two parameters, $\Theta_{1}$ and
$\Theta_{2}$. So, only two equations can be satisfied, i.e.,

$\displaystyle\alpha_{2,0}=1,$ (111) $\displaystyle\alpha_{2,2}=0.$ (112)

Note that the expression given by Eq. (108) can be rearranged by
$\alpha_{2,0}=1/2\big{[}1-\mathcal{A}\cos{(4\Theta_{1}-\beta)}\big{]}$, where
$\mathcal{A}=\sqrt{\cos^{2}{4\Theta_{2}}+\cos^{2}{\Gamma_{1}}\sin^{2}{4\Theta_{2}}}$.
Remarkably, the following inequality is always satisfied:

$\displaystyle\mathcal{A}=\sqrt{\cos^{2}{4\Theta_{2}}+\cos^{2}{\Gamma_{1}}\sin^{2}{4\Theta_{2}}}\leq\sqrt{\cos^{2}{4\Theta_{2}}+\sin^{2}{4\Theta_{2}}}=1.$
(113)

Only when $\Gamma_{1}=m\pi$ could there be an equality, where $m$ is an
arbitrary integer. In other words, the equation $\alpha_{2,0}=1$ has real
solutions if and only if $\Gamma_{1}=m\pi$. For simplicity, we set
$\Gamma_{1}=0$ in the main text, and Eqs. (111) and (112) become

$\displaystyle\sin^{2}2(\Theta_{1}-\Theta_{2})=1,8ex.$ (114)
$\displaystyle\pi^{2}\sin 2(\Theta_{2}-\Theta_{1})\Big{[}\sin
2(\Theta_{1}-\Theta_{2})\\!+\\!\cos\Theta_{2}(2\sin\Theta_{1}\\!+\\!\sin\Theta_{2})\\!+\\!\sin\Theta_{1}\cos\Theta_{1}\Big{]}=0.$
(115)

One solution of Eqs. (114) and (115) can be written as

$\begin{split}\Theta_{1}=m\pi+\frac{\pi}{8},~{}\Theta_{2}=m\pi+\frac{3\pi}{8},\end{split}$
(116)

where $m$ is an arbitrary integer. With this group of solutions, the final
propagator without errors becomes (in the basis
$\\{|g\rangle,|f\rangle,|e\rangle\\}$)

$\displaystyle U^{(2)}=i\left[\begin{array}[]{ccc}0&1&0\\\ 1&0&0\\\ 0&0&-i\\\
\end{array}\right].$ (120)

This propagator is recognized as an $X$ gate with the global phase factor $i$.
For the initial state of the system $|\psi_{i}\rangle=|g\rangle$, the final
state would become $|\psi_{f}\rangle=i|f\rangle$. Similarly, when the phase
difference $\Psi=0$, by this group of solutions, we could also obtain a $Y$
gate,

$\displaystyle U^{(2)}=i\left[\begin{array}[]{ccc}0&-i&0\\\ i&0&0\\\ 0&0&-i\\\
\end{array}\right]$ (124)

with the same global phase factor. For the initial state of the system
$|\psi_{i}\rangle=|g\rangle$, the final state would become
$|\psi_{f}\rangle=-|f\rangle$. Note that both the $X$ gate and the $Y$ gate
could achieve the population inversion, and the phase $\pi/2$ ($\pi$) of the
final state for $\Psi=\pi/2$ ($\Psi=0$) is the global phase, which can be
ignored. As a result, the transition probability of the state $|f\rangle$ is
only associated with its amplitude, and the Taylor series is

$P_{2}(\epsilon)=1+O(\epsilon^{4}),\\\ $

which is accurate to the fourth order in the pulse area error $\epsilon$.

## References

  * Krause  _et al._ (1998) J. L. Krause, D. H. Reitze, G. D. Sanders, A. V. Kuznetsov, and C. J. Stanton, Quantum control in quantum wells, Phys. Rev. B 57, 9024 (1998).
  * Ramakrishna  _et al._ (2000) V. Ramakrishna, K. L. Flores, H. Rabitz, and R. J. Ober, Quantum control by decompositions of SU(2), Phys. Rev. A 62, 053409 (2000).
  * Vilela Mendes and Man’ko (2003) R. Vilela Mendes and V. I. Man’ko, Quantum control and the Strocchi map, Phys. Rev. A 67, 053404 (2003).
  * Brańczyk  _et al._ (2007) A. M. Brańczyk, P. E. M. F. Mendonça, A. Gilchrist, A. C. Doherty, and S. D. Bartlett, Quantum control of a single qubit, Phys. Rev. A 75, 012329 (2007).
  * Wang and Niu (2020) R. Wang and Y.-Y. Niu, Quantum control of the dissociation of LiH molecules with two intense laser pulses, Phys. Rev. A 102, 013113 (2020).
  * Zhang  _et al._ (2020a) R. Zhang, Y.-F. Liu, and T. Chen, Non-Hermiticity-induced quantum control of localization in quantum walks, Phys. Rev. A 102, 022218 (2020a).
  * Dalgaard  _et al._ (2020) M. Dalgaard, F. Motzoi, J. H. M. Jensen, and J. Sherson, Hessian-based optimization of constrained quantum control, Phys. Rev. A 102, 042612 (2020).
  * Chen and Xue (2020) T. Chen and Z.-Y. Xue, High-Fidelity and Robust Geometric Quantum Gates that Outperform Dynamical Ones, Phys. Rev. Appl. 14, 064009 (2020).
  * Nielsen and Chuang (2000) M. A. Nielsen and I. L. Chuang,  _Quantum Computation and Quantum Information_ (Cambridge University Press, Cambridge, 2000).
  * Avinadav  _et al._ (2014) C. Avinadav, R. Fischer, P. London, and D. Gershoni, Time-optimal universal control of two-level systems under strong driving, Phys. Rev. B 89, 245311 (2014).
  * Børkje  _et al._ (2013) K. Børkje, A. Nunnenkamp, J. D. Teufel, and S. M. Girvin, Signatures of Nonlinear Cavity Optomechanics in the Weak Coupling Regime, Phys. Rev. Lett. 111, 053603 (2013).
  * Doria  _et al._ (2011) P. Doria, T. Calarco, and S. Montangero, Optimal Control Technique for Many-Body Quantum Dynamics, Phys. Rev. Lett. 106, 190501 (2011).
  * Caneva  _et al._ (2011) T. Caneva, T. Calarco, and S. Montangero, Chopped random-basis quantum optimization, Phys. Rev. A 84, 022326 (2011).
  * Baksic  _et al._ (2016) A. Baksic, H. Ribeiro, and A. A. Clerk, Speeding up Adiabatic Quantum State Transfer by Using Dressed States, Phys. Rev. Lett. 116, 230503 (2016).
  * Wang  _et al._ (2015) X. Wang, M. Allegra, K. Jacobs, S. Lloyd, C. Lupo, and M. Mohseni, Quantum Brachistochrone Curves as Geodesics: Obtaining Accurate Minimum-Time Protocols for the Control of Quantum Systems, Phys. Rev. Lett. 114, 170501 (2015).
  * Bao  _et al._ (2018) S. Bao, S. Kleer, R. Wang, and A. Rahmani, Optimal control of superconducting gmon qubits using Pontryagin’s minimum principle: Preparing a maximally entangled state with singular bang-bang protocols, Phys. Rev. A 97, 062343 (2018).
  * Chen  _et al._ (2017) Y.-H. Chen, Z.-C. Shi, J. Song, Y. Xia, and S.-B. Zheng, Optimal shortcut approach based on an easily obtained intermediate Hamiltonian, Phys. Rev. A 95, 062319 (2017).
  * Leung  _et al._ (2017) N. Leung, M. Abdelhafez, J. Koch, and D. Schuster, Speedup for quantum optimal control from automatic differentiation based on graphics processing units, Phys. Rev. A 95, 042318 (2017).
  * Wu  _et al._ (2018) R.-B. Wu, B. Chu, D. H. Owens, and H. Rabitz, Data-driven gradient algorithm for high-precision quantum control, Phys. Rev. A 97, 042122 (2018).
  * Zhang and Yang (2021) F.-Y. Zhang and C.-P. Yang, Generation of generalized hybrid entanglement in cavity electro–optic systems, Quantum Sci. Technol. 6, 025003 (2021).
  * Lin and Li (2009) Q. Lin and J. Li, Quantum control gates with weak cross-Kerr nonlinearity, Phys. Rev. A 79, 022301 (2009).
  * Wu  _et al._ (2008) R. Wu, R. Chakrabarti, and H. Rabitz, Optimal control theory for continuous-variable quantum gates, Phys. Rev. A 77, 052303 (2008).
  * Kang  _et al._ (2018) Y.-H. Kang, Y.-H. Chen, Z.-C. Shi, B.-H. Huang, J. Song, and Y. Xia, Nonadiabatic holonomic quantum computation using Rydberg blockade, Phys. Rev. A 97, 042336 (2018).
  * Kosut  _et al._ (2013) R. L. Kosut, M. D. Grace, and C. Brif, Robust control of quantum gates via sequential convex programming, Phys. Rev. A 88, 052326 (2013).
  * Hegde  _et al._ (2020) S. S. Hegde, J. Zhang, and D. Suter, Efficient Quantum Gates for Individual Nuclear Spin Qubits by Indirect Control, Phys. Rev. Lett. 124, 220501 (2020).
  * Wu  _et al._ (2020) X. Wu, S. L. Tomarken, N. A. Petersson, L. A. Martinez, Y. J. Rosen, and J. L. DuBois, High-Fidelity Software-Defined Quantum Logic on a Superconducting Qudit, Phys. Rev. Lett. 125, 170502 (2020).
  * Guo  _et al._ (2020) C.-Y. Guo, L.-L. Yan, S. Zhang, S.-L. Su, and W. Li, Optimized geometric quantum computation with a mesoscopic ensemble of Rydberg atoms, Phys. Rev. A 102, 042607 (2020).
  * Dridi  _et al._ (2020a) G. Dridi, K. Liu, and S. Guérin, Optimal Robust Quantum Control by Inverse Geometric Optimization, Phys. Rev. Lett. 125, 250403 (2020a).
  * Kang  _et al._ (2020) Y.-H. Kang, Z.-C. Shi, B.-H. Huang, J. Song, and Y. Xia, Flexible scheme for the implementation of nonadiabatic geometric quantum computation, Phys. Rev. A 101, 032322 (2020).
  * Chen  _et al._ (2020) T. Chen, P. Shen, and Z.-Y. Xue, Robust and Fast Holonomic Quantum Gates with Encoding on Superconducting Circuits, Phys. Rev. Appl. 14, 034038 (2020).
  * Zhang  _et al._ (2020b) C. Zhang, T. Chen, S. Li, X. Wang, and Z.-Y. Xue, High-fidelity geometric gate for silicon-based spin qubits, Phys. Rev. A 101, 052302 (2020b).
  * Xue  _et al._ (2016) Z.-Y. Xue, J. Zhou, Y.-M. Chu, and Y. Hu, Nonadiabatic holonomic quantum computation with all-resonant control, Phys. Rev. A 94, 022331 (2016).
  * Ge  _et al._ (2020) X. Ge, H. Ding, H. Rabitz, and R.-B. Wu, Robust quantum control in games: An adversarial learning approach, Phys. Rev. A 101, 052317 (2020).
  * Zhong and Jing (2020) Y. Zhong and J. Jing, Enhancement of tripartite quantum correlation by coherent feedback control, Phys. Rev. A 101, 023813 (2020).
  * Ma  _et al._ (2020a) A. Ma, A. B. Magann, T.-S. Ho, and H. Rabitz, Optimal control of coupled quantum systems based on the first-order Magnus expansion: Application to multiple dipole-dipole-coupled molecular rotors, Phys. Rev. A 102, 013115 (2020a).
  * Takeda  _et al._ (2020) K. Takeda, A. Noiri, J. Yoneda, T. Nakajima, and S. Tarucha, Resonantly Driven Singlet-Triplet Spin Qubit in Silicon, Phys. Rev. Lett. 124, 117701 (2020).
  * Goldman  _et al._ (2020) M. L. Goldman, T. L. Patti, D. Levonian, S. F. Yelin, and M. D. Lukin, Optical Control of a Single Nuclear Spin in the Solid State, Phys. Rev. Lett. 124, 153203 (2020).
  * Hayes  _et al._ (2020) D. Hayes, D. Stack, B. Bjork, A. C. Potter, C. H. Baldwin, and R. P. Stutz, Eliminating Leakage Errors in Hyperfine Qubits, Phys. Rev. Lett. 124, 170501 (2020).
  * Allen and Eberly (1987) L. Allen and J. H. Eberly,  _Optical resonance and two-level atoms_ (Dover Publications, New York, 1987).
  * Melinger  _et al._ (1992) J. S. Melinger, S. R. Gandhi, A. Hariharan, J. X. Tull, and W. S. Warren, Generation of narrowband inversion with broadband laser pulses, Phys. Rev. Lett. 68, 2000 (1992).
  * Melinger  _et al._ (1994) J. S. Melinger, S. R. Gandhi, A. Hariharan, D. Goswami, and W. S. Warren, Adiabatic population transfer with frequency-swept laser pulses, J. Chem. Phys. 101, 6439 (1994).
  * Solá  _et al._ (1999) I. R. Solá, V. S. Malinovsky, B. Y. Chang, J. Santamaria, and K. Bergmann, Coherent population transfer in three-level $\Lambda$ systems by chirped laser pulses: Minimization of the intermediate-level population, Phys. Rev. A 59, 4494 (1999).
  * Band and Magnes (1994) Y. B. Band and O. Magnes, Chirped adiabatic passage with temporally delayed pulses, Phys. Rev. A 50, 584 (1994).
  * Vitanov  _et al._ (2001) N. Vitanov, T. Halfmann, B. Shore, and K. Bergmann, Laser-induced population transfer by adiabatic passage techniques, Annu. rev. Phys. Chem. 52, 763 (2001).
  * Goswami (2003) D. Goswami, Optical pulse shaping approaches to coherent control, Phys. Rep. 374, 385 (2003).
  * Gaubatz  _et al._ (1990) U. Gaubatz, P. Rudecki, S. Schiemann, and K. Bergmann, Population transfer between molecular vibrational levels by stimulated Raman scattering with partially overlapping laser fields. A new concept and experimental results, J. Chem. Phys. 92, 5363 (1990).
  * Marte  _et al._ (1991) P. Marte, P. Zoller, and J. L. Hall, Coherent atomic mirrors and beam splitters by adiabatic passage in multilevel systems, Phys. Rev. A 44, R4118 (1991).
  * Shore  _et al._ (1992) B. W. Shore, K. Bergmann, A. Kuhn, S. Schiemann, J. Oreg, and J. H. Eberly, Laser-induced population transfer in multistate systems: A comparative study, Phys. Rev. A 45, 5297 (1992).
  * Vitanov  _et al._ (2017) N. V. Vitanov, A. A. Rangelov, B. W. Shore, and K. Bergmann, Stimulated Raman adiabatic passage in physics, chemistry, and beyond, Rev. Mod. Phys. 89, 015006 (2017).
  * Guéry-Odelin  _et al._ (2019) D. Guéry-Odelin, A. Ruschhaupt, A. Kiely, E. Torrontegui, S. Martínez-Garaot, and J. G. Muga, Shortcuts to adiabaticity: Concepts, methods, and applications, Rev. Mod. Phys. 91, 045001 (2019).
  * Chen  _et al._ (2011) X. Chen, E. Torrontegui, and J. G. Muga, Lewis-Riesenfeld invariants and transitionless quantum driving, Phys. Rev. A 83, 062116 (2011).
  * Chen  _et al._ (2010) X. Chen, I. Lizuain, A. Ruschhaupt, D. Guéry-Odelin, and J. G. Muga, Shortcut to Adiabatic Passage in Two- and Three-Level Atoms, Phys. Rev. Lett. 105, 123003 (2010).
  * Lewis and Riesenfeld (1969) H. R. Lewis and W. B. Riesenfeld, An Exact Quantum Theory of the Time-Dependent Harmonic Oscillator and of a Charged Particle in a Time-Dependent Electromagnetic Field, J. Math. Phys. 10, 1458 (1969).
  * Demirplak and Rice (2003) M. Demirplak and S. A. Rice, Adiabatic Population Transfer with Control Fields, J. Phys. Chem. A 107, 9937 (2003).
  * Demirplak and Rice (2005) M. Demirplak and S. A. Rice, Assisted Adiabatic Passage Revisited, J. Phys. Chem. B 109, 6838 (2005).
  * Berry (2009) M. V. Berry, Transitionless quantum driving, J. Phys. A 42, 365303 (2009).
  * Lin  _et al._ (2020a) C. Lin, D. Sels, and Y. Wang, Time-optimal control of a dissipative qubit, Phys. Rev. A 101, 022320 (2020a).
  * Holland  _et al._ (2020) E. T. Holland, K. A. Wendt, K. Kravvaris, X. Wu, W. E. Ormand, J. L. DuBois, S. Quaglioni, and F. Pederiva, Optimal control for the quantum simulation of nuclear dynamics, Phys. Rev. A 101, 062307 (2020).
  * Lin  _et al._ (2020b) C. Lin, D. Sels, Y. Ma, and Y. Wang, Stochastic optimal control formalism for an open quantum system, Phys. Rev. A 102, 052605 (2020b).
  * Larocca  _et al._ (2020) M. Larocca, E. Calzetta, and D. A. Wisniacki, Exploiting landscape geometry to enhance quantum optimal control, Phys. Rev. A 101, 023410 (2020).
  * Tian  _et al._ (2020) J. Tian, H. Liu, Y. Liu, P. Yang, R. Betzholz, R. S. Said, F. Jelezko, and J. Cai, Quantum optimal control using phase-modulated driving fields, Phys. Rev. A 102, 043707 (2020).
  * An  _et al._ (2021) Z. An, H.-J. Song, Q.-K. He, and D. L. Zhou, Quantum optimal control of multilevel dissipative quantum systems with reinforcement learning, Phys. Rev. A 103, 012404 (2021).
  * Plantenberg  _et al._ (2007) J. H. Plantenberg, P. C. de Groot, C. J. P. M. Harmans, and J. E. Mooij, Demonstration of controlled-NOT quantum gates on a pair of superconducting quantum bits, Nature 447, 836 (2007).
  * Glaser  _et al._ (2015) S. J. Glaser, U. Boscain, T. Calarco, C. P. Koch, W. Köckenberger, R. Kosloff, I. Kuprov, B. Luy, S. Schirmer, T. Schulte-Herbrüggen, D. Sugny, and F. K. Wilhelm, Training Schrödinger’s cat: Quantum optimal control, Eur. Phys. J. D 69 (2015).
  * Wimperis (1994) S. Wimperis, Broadband, Narrowband, and Passband Composite Pulses for Use in Advanced NMR Experiments, J. Magn. Reson. 109, 221 (1994).
  * Levitt (1986) M. H. Levitt, Composite pulses, Prog. NMR Spectrosc. 18, 61 (1986).
  * Wimperis (1991) S. Wimperis, Iterative schemes for phase-distortionless composite 180° pulses, J. Magn. Reson. 93, 199 (1991).
  * Husain  _et al._ (2013) S. Husain, M. Kawamura, and J. A. Jones, Further analysis of some symmetric and antisymmetric composite pulses for tackling pulse strength errors, J. Magn. Reson. 230, 145 (2013).
  * Genov  _et al._ (2014) G. T. Genov, D. Schraft, T. Halfmann, and N. V. Vitanov, Correction of Arbitrary Field Errors in Population Inversion of Quantum Systems by Universal Composite Pulses, Phys. Rev. Lett. 113, 043001 (2014).
  * Ivanov and Vitanov (2011) S. S. Ivanov and N. V. Vitanov, High-fidelity local addressing of trapped ions and atoms by composite sequences of laser pulses, Opt. Lett. 36, 1275 (2011).
  * Torosov and Vitanov (2019) B. T. Torosov and N. V. Vitanov, Arbitrarily accurate variable rotations on the Bloch sphere by composite pulse sequences, Phys. Rev. A 99, 013402 (2019).
  * Torosov and Vitanov (2011) B. T. Torosov and N. V. Vitanov, Smooth composite pulses for high-fidelity quantum information processing, Phys. Rev. A 83, 053420 (2011).
  * Jones (2013) J. A. Jones, Designing short robust NOT gates for quantum computation, Phys. Rev. A 87, 052317 (2013).
  * Genov  _et al._ (2020) G. T. Genov, M. Hain, N. V. Vitanov, and T. Halfmann, Universal composite pulses for efficient population inversion with an arbitrary excitation profile, Phys. Rev. A 101, 013827 (2020).
  * Dridi  _et al._ (2020b) G. Dridi, M. Mejatty, S. J. Glaser, and D. Sugny, Robust control of a not gate by composite pulses, Phys. Rev. A 101, 012321 (2020b).
  * Torosov  _et al._ (2021) B. T. Torosov, B. W. Shore, and N. V. Vitanov, Coherent control techniques for two-state quantum systems: A comparative study, Phys. Rev. A 103, 033110 (2021).
  * Torosov  _et al._ (2020a) B. T. Torosov, S. S. Ivanov, and N. V. Vitanov, Narrowband and passband composite pulses for variable rotations, Phys. Rev. A 102, 013105 (2020a).
  * Kyoseva  _et al._ (2019) E. Kyoseva, H. Greener, and H. Suchowski, Detuning-modulated composite pulses for high-fidelity robust quantum control, Phys. Rev. A 100, 032333 (2019).
  * Wang  _et al._ (2012) X. Wang, L. S. Bishop, J. Kestner, E. Barnes, K. Sun, and S. D. Sarma, Composite pulses for robust universal control of singlet–triplet qubits, Nat. Commun. 3 (2012).
  * Wang  _et al._ (2014) X. Wang, L. S. Bishop, E. Barnes, J. P. Kestner, and S. D. Sarma, Robust quantum gates for singlet-triplet spin qubits using composite pulses, Phys. Rev. A 89, 022310 (2014).
  * Torosov and Vitanov (2020) B. T. Torosov and N. V. Vitanov, High-fidelity composite quantum gates for Raman qubits, Phys. Rev. Research 2, 043194 (2020).
  * Shi  _et al._ (2021) Z.-C. Shi, H.-N. Wu, L.-T. Shen, J. Song, Y. Xia, X. X. Yi, and S.-B. Zheng, Robust single-qubit gates by composite pulses in three-level systems, Phys. Rev. A 103, 052612 (2021).
  * Torosov  _et al._ (2020b) B. T. Torosov, M. Drewsen, and N. V. Vitanov, Chiral resolution by composite Raman pulses, Phys. Rev. Research 2, 043235 (2020b).
  * Mansourzadeh-Ashkani  _et al._ (2021) N. Mansourzadeh-Ashkani, M. Saadati-Niari, F. Zolfagharpour, and B. Nedaee-Shakarab, Nuclear-state population transfer using composite stimulated Raman adiabatic passage, Nucl. Phys. A 1007, 122119 (2021).
  * Ghosh  _et al._ (2017) J. Ghosh, S. N. Coppersmith, and M. Friesen, Pulse sequences for suppressing leakage in single-qubit gate operations, Phys. Rev. B 95, 241307 (2017).
  * Friesen  _et al._ (2017) M. Friesen, J. Ghosh, M. A. Eriksson, and S. N. Coppersmith, A decoherence-free subspace in a charge quadrupole qubit, Nat. Commun. 8, 15923 (2017).
  * Shi  _et al._ (2022) Z.-C. Shi, C. Zhang, D. Ran, Y. Xia, R. Ianconescu, A. Friedman, X. X. Yi, and S.-B. Zheng, Composite pulses for high fidelity population transfer in three-level systems, New J. Phys. 24, 023014 (2022).
  * Pillet  _et al._ (1993) P. Pillet, C. Valentin, R.-L. Yuan, and J. Yu, Adiabatic population transfer in a multilevel system, Phys. Rev. A 48, 845 (1993).
  * Featonby  _et al._ (1996) P. D. Featonby, G. S. Summy, J. L. Martin, H. Wu, K. P. Zetie, C. J. Foot, and K. Burnett, Adiabatic transfer for atomic interferometry, Phys. Rev. A 53, 373 (1996).
  * Král and Shapiro (2001) P. Král and M. Shapiro, Cyclic Population Transfer in Quantum Systems with Broken Symmetry, Phys. Rev. Lett. 87, 183002 (2001).
  * Lacour  _et al._ (2006) X. Lacour, N. Sangouard, S. Guérin, and H. R. Jauslin, Arbitrary state controlled-unitary gate by adiabatic passage, Phys. Rev. A 73, 042321 (2006).
  * Klein  _et al._ (2007) J. Klein, F. Beil, and T. Halfmann, Robust Population Transfer by Stimulated Raman Adiabatic Passage in a ${\mathrm{Pr}}^{3+}\mathrm{\text{:}}{\mathrm{Y}}_{2}{\mathrm{SiO}}_{5}$ Crystal, Phys. Rev. Lett. 99, 113003 (2007).
  * Klein  _et al._ (2008) J. Klein, F. Beil, and T. Halfmann, Experimental investigations of stimulated Raman adiabatic passage in a doped solid, Phys. Rev. A 78, 033416 (2008).
  * Longhi (2009) S. Longhi, Quantum-optical analogies using photonic structures, Laser Photonics Rev. 3, 243 (2009).
  * Wang and Clerk (2012) Y.-D. Wang and A. A. Clerk, Using dark modes for high-fidelity optomechanical quantum state transfer, New J. Phys. 14, 105010 (2012).
  * Scully and Zubairy (1997) M. O. Scully and M. S. Zubairy,  _Quantum Optics_ (Cambridge University Press, Cambridge, 1997).
  * Brion  _et al._ (2007) E. Brion, L. H. Pedersen, and K. Mølmer, Adiabatic elimination in a lambda system, J. Phys. A: Math. Theor. 40, 1033 (2007).
  * Shevchenko  _et al._ (2010) S. Shevchenko, S. Ashhab, and F. Nori, Landau-Zener-Stückelberg interferometry, Phys. Rep. 492, 1 (2010).
  * Politi  _et al._ (2008) A. Politi, M. J. Cryan, J. G. Rarity, S. Yu, and J. L. O’Brien, Silica-on-Silicon Waveguide Quantum Circuits, Science 320, 646 (2008).
  * Reed  _et al._ (2016) M. D. Reed, B. M. Maune, R. W. Andrews, M. G. Borselli, K. Eng, M. P. Jura, A. A. Kiselev, T. D. Ladd, S. T. Merkel, I. Milosavljevic, E. J. Pritchett, M. T. Rakher, R. S. Ross, A. E. Schmitz, A. Smith, J. A. Wright, M. F. Gyure, and A. T. Hunter, Reduced Sensitivity to Charge Noise in Semiconductor Spin Qubits via Symmetric Operation, Phys. Rev. Lett. 116, 110402 (2016).
  * Xiang  _et al._ (2013) Z.-L. Xiang, S. Ashhab, J. Q. You, and F. Nori, Hybrid quantum circuits: Superconducting circuits interacting with other quantum systems, Rev. Mod. Phys. 85, 623 (2013).
  * Xia and Xie (2008) L.-X. Xia and Q.-T. Xie, Unconventional geometric quantum phase gates with two SQUIDs in a cavity, Opt. Commun. 281, 2700 (2008).
  * Yang  _et al._ (2003) C.-P. Yang, S.-I. Chu, and S. Han, Possible realization of entanglement, logical gates, and quantum-information transfer with superconducting-quantum-interference-device qubits in cavity QED, Phys. Rev. A 67, 042311 (2003).
  * Yang  _et al._ (2004) C.-P. Yang, S.-I. Chu, and S. Han, Quantum Information Transfer and Entanglement with SQUID Qubits in Cavity QED: A Dark-State Scheme with Tolerance for Nonuniform Device Parameter, Phys. Rev. Lett. 92, 117902 (2004).
  * Kang  _et al._ (2016) Y.-H. Kang, Y.-H. Chen, Z.-C. Shi, J. Song, and Y. Xia, Fast preparation of $W$ states with superconducting quantum interference devices by using dressed states, Phys. Rev. A 94, 052311 (2016).
  * Guarcello  _et al._ (2020) C. Guarcello, R. Citro, O. Durante, F. S. Bergeret, A. Iorio, C. Sanz-Fernández, E. Strambini, F. Giazotto, and A. Braggio, rf-SQUID measurements of anomalous Josephson effect, Phys. Rev. Research 2, 023165 (2020).
  * Braumüller  _et al._ (2020) J. Braumüller, L. Ding, A. P. Vepsäläinen, Y. Sung, M. Kjaergaard, T. Menke, R. Winik, D. Kim, B. M. Niedzielski, A. Melville, J. L. Yoder, C. F. Hirjibehedin, T. P. Orlando, S. Gustavsson, and W. D. Oliver, Characterizing and Optimizing Qubit Coherence Based on SQUID Geometry, Phys. Rev. Appl. 13, 054079 (2020).
  * Stern  _et al._ (2014) M. Stern, G. Catelani, Y. Kubo, C. Grezes, A. Bienfait, D. Vion, D. Esteve, and P. Bertet, Flux Qubits with Long Coherence Times for Hybrid Quantum Circuits, Phys. Rev. Lett. 113, 123601 (2014).
  * Leek  _et al._ (2010) P. J. Leek, M. Baur, J. M. Fink, R. Bianchetti, L. Steffen, S. Filipp, and A. Wallraff, Cavity Quantum Electrodynamics with Separate Photon Storage and Qubit Readout Modes, Phys. Rev. Lett. 104, 100504 (2010).
  * Yan  _et al._ (2016) F. Yan, S. Gustavsson, A. Kamal, J. Birenbaum, A. P. Sears, D. Hover, T. J. Gudmundsen, D. Rosenberg, G. Samach, S. Weber, J. L. Yoder, T. P. Orlando, J. Clarke, A. J. Kerman, and W. D. Oliver, The flux qubit revisited to enhance coherence and reproducibility, Nat. Commun. 7, 12964 (2016).
  * Ma  _et al._ (2020b) Y. Ma, X. Pan, W. Cai, X. Mu, Y. Xu, L. Hu, W. Wang, H. Wang, Y. P. Song, Z.-B. Yang, S.-B. Zheng, and L. Sun, Manipulating Complex Hybrid Entanglement and Testing Multipartite Bell Inequalities in a Superconducting Circuit, Phys. Rev. Lett. 125, 180503 (2020b).
  * Xu  _et al._ (2020) Y. Xu, Y. Ma, W. Cai, X. Mu, W. Dai, W. Wang, L. Hu, X. Li, J. Han, H. Wang, Y. P. Song, Z.-B. Yang, S.-B. Zheng, and L. Sun, Demonstration of Controlled-Phase Gates between Two Error-Correctable Photonic Qubits, Phys. Rev. Lett. 124, 120501 (2020).
  * Niemczyk  _et al._ (2010) T. Niemczyk, F. Deppe, H. Huebl, E. P. Menzel, F. Hocke, M. J. Schwarz, J. J. Garcia-Ripoll, D. Zueco, T. Hümmer, E. Solano, A. Marx, and R. Gross, Circuit quantum electrodynamics in the ultrastrong-coupling regime, Nat. Phys. 6, 772 (2010).