\section{Introduction}

Hybrid quantum systems have attracted great attentions due to their diversified applications in quantum computing~\cite{quantumcomputing}, quantum communications~\cite{quantumcommunication}, and quantum sensing~\cite{quantumsense}. During the past decade, the hybrid quantum systems consisting of the collective spin excitations in the ferromagnetic crystals have been applied to novel quantum technologies~\cite{magnon,magnon2,magnon4,magnon5,magnon6,magnoncavity} by virtue of the long coherent-time of the magnet-spin ensemble. The strong dipole-transition for efficiently coupling to the microwave photons and phonons allows the construction of the hybrid magnon-photon, magnon-phonon, and magnon-photon-phonon systems in both theoretical proposals~\cite{yigcavity1,yigcavity2,mppentang} and experimental demonstrations~\cite{yigcavity3,yigcavity4,yigcavity5,yigcavity6}. In more recent experiments~\cite{magnonqubit,magnonqubit2}, a hybrid qubit-photon-magnon system has been realized, opening a new avenue to studying the intermediate transitions by the interaction Hamiltonian of different constituents. It is expected to have many interesting applications, including providing a platform to generate hybrid entangled states.

The entangled states~\cite{quantumentanglement} are essential ingredients for extracting quantum advantage in the quantum communication protocols~\cite{quantumcommun}, such as quantum key distribution~\cite{keydistribute}, quantum secret sharing~\cite{quantumsecretsharing}, and quantum secure direct communication~\cite{quantumsecurecommun}. Many quantum platforms~\cite{trapion,ghzstate,superconduct2,superconduct3,Schrodingerstate} and protocols for preparing and measuring the entangled states have therefore been intensively pursued for a long time~\cite{GHZ,ghzstate,GHZstate2,bellstate3,multiphoton,photonic,noon1,noon2} and are still under active investigation. The simplest and the most popular maximally-entangled-states are Bell states involving only two discrete systems~\cite{quantumentanglement,bellstate}, classified to two groups, i.e., $|\psi_\pm\rangle=(|01\rangle\pm|10\rangle)/\sqrt{2}$ and $|\varphi_\pm\rangle=(|00\rangle\pm|11\rangle)/\sqrt{2}$. Here $|0\rangle$ and $|1\rangle$ are respectively the ground and the excited states of a two-level system. The Bell states are fundamental and important in both quantum cryptography~\cite{cryptography,cryptography2} and quantum teleportation~\cite{teleportation}. In a recent experiment~\cite{yigcavity3}, the single-excitation Bell states $|\psi_\pm\rangle$ of the photon-magnon system have been steadily generated under the parity-time-symmetry broken condition. However, how to generate the double-excitation Bell states $|\varphi_\pm\rangle$ in a hybrid photon-magnon system is still an open question. More generally, the maximally entangled states for the $N$-partite discrete quantum systems, such as the Greenberger-Horne-Zeilinger (GHZ) states~\cite{GHZ,ghzstate}, i.e., $(|0\rangle^{\otimes N}+|1\rangle^{\otimes N})/\sqrt{2}$ and the Werner (W) states~\cite{wstate,wernerstate}, i.e., $(|100\cdots0\rangle+|010\cdots0\rangle+\cdots+|000\cdots1\rangle)/\sqrt{N}$, are also of great practical importance. Note that the multiple GHZ states are natural extensions of the double-excitation Bell states $|\varphi_\pm\rangle$ and cannot be converted to the W-states with nonzero probability~\cite{quantumentanglement}.

This work focuses on generating the local Bell states $|\varphi_\pm\rangle$ and the global GHZ states in a tripartite hybrid system consisting of a superconducting qubit and a magnon mode embedded in a microwave cavity (see the diagram in Fig.~\ref{diagram}), taking advantage of the tunability of the transition frequencies of qubit~\cite{cqed} and magnon~\cite{magnon} and their respective strong or ultrastrong couplings to the cavity mode (resonator). Due to the counter-rotating terms with no conservation of the excitation-number in the Rabi interactions~\cite{cq,ultrastrong,ultrastrong2}, it is possible to design an effective Hamiltonian as well as the generation protocol for connecting states in the whole Hilbert space nearby the desired multi-excitation resonant points. The effective Hamiltonian for the interested indirect transitions can be extracted using the high-order Fermi's Golden rule~\cite{fermigolden,secondorder,noon3,ae,ae1}. It states that if the shortest path connecting $|i\rangle$ and $|f\rangle$ (two eigenstates of the free Hamiltonian) is an $n$th-order process, then the effective coupling strength (in the leading order) between them reads,
\begin{equation}\label{Fermi}
g_{\rm eff}=\sum_{m_1 m_2...m_{n-1}}\frac{V_{fm_{n-1}}...V_{m_2m_1}V_{m_1i}}{(E_i-E_{m_{n-1}})...(E_i-E_{m_1})},
\end{equation}
where $V_{m_2m_1}\equiv\langle m_2|V|m_1\rangle$ is the transition matrix element of the interaction Hamiltonian and $E_i$ is the eigenvalue of the $i$th eigenstate of the free Hamiltonian (with no degeneracy).

The rest part of the work is organized as follows. In Sec.~\ref{model}, we introduce the total Hamiltonian of the hybrid qubit-photon-magnon system. In Sec.~\ref{bell}, we show that at a desired point, where the transition frequency of the qubit is nearly resonant with the frequency sum of the cavity and magnon modes, an effective Hamiltonian can be constructed to describe a photon and a magnon simultaneously excited by annihilating an atomic excitation. This ``three-wave-mixing'' Hamiltonian is applied to generate the double-excitation Bell states of the photon-magnon subsystem in Sec.~\ref{Bellstep}. And then in Sec.~\ref{Bellnum}, the fidelity of generation is numerically estimated with experimental parameters. In Sec.~\ref{GHZ}, the working point of our protocol moves to the near-resonant point for the transition frequency of the qubit and the detuning of the photon and magnon modes, where one can construct an effective Hamiltonian to simultaneously excite a qubit and a magnon by one photon. Then we show that the GHZ state of the whole hybrid system can be concisely generated by a three-step scheme in Sec.~\ref{GHZstep}, whose fidelity under dissipation is evaluated in Sec.~\ref{GHZnum}. The protocol in Sec.~\ref{GHZstep} is slightly modified to generate the Bell state of the qubit-magnon subsystem in Sec.~\ref{discussion}. We finally conclude this work in Sec.~\ref{conclusion}. The derivation details for the effective Hamiltonians in Secs.~\ref{bell} and \ref{GHZ} can be found respectively in Appendices~\ref{appa} and \ref{appb}.

\section{Model Hamiltonian of hybrid system}\label{model}

The hybrid quantum model we considered in this work consists of a superconducting qubit, a microwave cavity/resonator and a ferromagnetic crystal in the Kittle mode, as shown in Fig.~\ref{diagram}. The interaction between the qubit and the resonator is described by a general Rabi model, and simultaneously the magnetostatic mode (magnon) is coupled to the resonator via the magnetic dipole interaction. Thus the Hamiltonian for the whole system~\cite{magnonqubit,magnonqubit2,ultrastrong2,secondorder} $(\hbar\equiv1)$ can be written as
\begin{equation}\label{Hamiltonian}
\begin{aligned}
&H=H_0+V,\\
&H_0=\omega_aa^\dag a+\omega_mm^\dag m+\omega_q\sigma_+\sigma_-,\\
&V=g(a+a^\dag)(m+m^\dag)+G(a+a^\dag)(\sigma_x \cos\theta+\sigma_z\sin\theta).
\end{aligned}
\end{equation}
Here $a (a^\dag)$ and $m (m^\dag)$ are the annihilation (creation) operators of the photon and magnon modes, respectively. $\omega_a$ and $\omega_m$ are their respective eigen-frequencies. $g$ is the coupling strength between the magnon and the photon mode. The two levels of the qubit are labeled by $|g\rangle$ and $|e\rangle$, indicating the ground and the excited states, respectively. The Pauli operators for the qubit are then written as $\sigma_+\equiv|e\rangle\langle g|$, $\sigma_-\equiv|g\rangle\langle e|$, $\sigma_x\equiv |e\rangle\langle g|+|g\rangle\langle e|$, and $\sigma_z\equiv|e\rangle\langle e|-|g\rangle\langle g|$. We do not apply the rotating-wave approximation to the interaction Hamiltonian $V$. The angle $\theta$ parameterizes the amount of the longitudinal and the transversal couplings between the qubit and the resonator, which is adjustable and independent of the coupling strength $G$. Arbitrary mixture of the longitudinal and the transversal couplings has been realized in the circuit-QED experiments~\cite{cq,ultrastrong}, where the coupling strengths can reach the ultrastrong coupling regime.

\begin{figure}[htbp]
\centering
\includegraphics[width=0.45\textwidth]{diagram.eps}
\caption{A tripartite hybrid system: a single-mode cavity (resonator) is coupled to a superconducting qubit and a YIG sphere (in the Kittel mode, a uniform magnetostatic mode) with the coupling strengths $G$ and $g$, respectively. }\label{diagram}
\end{figure}

Our approach works in the dispersive regime of the hybrid model in Eq.~(\ref{Hamiltonian}), where the coupling strengths are much smaller than the transition frequencies of the three subsystems and their detunings, i.e., $g,G\ll\omega_a$, $\omega_m$, $\omega_q$, $|\omega_a-\omega_m|$, $|\omega_a-\omega_q|$, $|\omega_m-\omega_q|$. However, it is crucial to know precisely their respective ranges of validity. We would apply the high-order Fermi's Golden rule in Eq.~(\ref{Fermi}) based on the standard perturbation theory to obtain the effective Hamiltonians at the near-resonant points in charge of the desired Rabi oscillations, through which we can generate the Bell states of the photon-magnon subsystem and the GHZ states of the whole qubit-photon-magnon system. For each effective Hamiltonian, a pair of the effective coupling strength $g_{\rm eff}$ and the energy shift $\delta$ in the leading order of $g$ and $G$ can be analytically derived. In Appendices~\ref{appa} and \ref{appb}, we will benchmark the ranges of validity of the coupling strengths $g,G$ by comparing the analytical results of $g_{\rm eff}$ and $\delta$ by the effective Hamiltonian and their numerical simulation over the whole Hilbert space.

\section{Generating Bell state of photon-magnon subsystem}

\subsection{Effective Hamiltonian}\label{bell}

In this section, we propose a protocol to generate the double-excitation Bell state of photon-magnon subsystem. The basic mechanism in this protocol is similar to the three-wave mixing schemes~\cite{threewave1,threewave2} in the nonlinear quantum optics. At a special point in the parametric space, the de-excitation of the superconducting qubit gives rise to a magnon-photon pair. State transfer occurs by the Rabi oscillation  $|e00\rangle\equiv|e\rangle_q|0\rangle_a|0\rangle_m\leftrightarrow|g11\rangle$, where $|e00\rangle$ and $|g11\rangle$ are two eigenstates of the free Hamiltonian $H_0$ in Eq.~(\ref{Hamiltonian}).

\begin{figure}[htbp]
\centering
\includegraphics[width=0.45\textwidth]{e1p1m1eigen.eps}
\caption{Energy levels and the avoided level crossing for the states $|e00\rangle$ and $|g11\rangle$, normalized by $\omega_a$ and plotted as functions of the qubit transition frequency $\omega_q$. The avoided level crossings of these two eigenstates are distinguished by the dark circle. Here the parameters are fixed as $\omega_m=1.7\omega_a$, $\theta=\pi/4$, and $g=G=0.1\omega_a$. }\label{e1p1m1eigen}
\end{figure}

As shown in Fig.~\ref{e1p1m1eigen}, the interaction Hamiltonian $V$ in Eq.~(\ref{Hamiltonian}) shifts the point of the avoided level crossing for $|e00\rangle$ and $|g11\rangle$ from the exact double-resonant point by $\delta$, which satisfies $\omega_q=\omega_a+\omega_m+\delta$. In Fig.~\ref{e1p1m1eigen}, the eigenvalues $E_n$ of both states are obtained by a standard numerical diagonalization over the full Hamiltonian in a truncated Hilbert space. A sufficiently large number of energy eigenstates have been used to ensure that the result is not significantly affected by the truncation. The avoided level crossing is distinguished with a black circle, demonstrating the nonlinear resonance between the states $|e00\rangle$ and $|g11\rangle$ with an effective transition rate $g_{\rm eff}$. This phenomenon can be well described by the effective Hamiltonian up to the leading-order process involving all the coupling strengths, provided that $g,G\ll\omega_a$, $\omega_m$, $|\omega_a-\omega_m|$. The derivation details to obtain the effective Hamiltonian are provided in Appendix~\ref{appa}. In the subspace spanned by $\{|e00\rangle, |g11\rangle\}$, we have
\begin{equation}\label{Heff1}
H_{\rm eff}=g_{\rm eff}\left(|e00\rangle\langle g11|+|g11\rangle\langle e00|\right),
\end{equation}
at the resonant point $\omega_q=\omega_a+\omega_m+\delta$. We find that the magnitudes of the energy shift
\begin{equation}\label{delta1}
\delta=-2G^2\cos^2\theta\left(\frac{1}{\omega_m}+\frac{1}{2\omega_a+\omega_m}\right)-\frac{2g^2}{\omega_a+\omega_m},
\end{equation}
and the coupling strength $g_{\rm eff}$
\begin{equation}\label{geff1}
g_{\rm eff}=\frac{2G^2g\sin(2\theta)}{\omega_m(\omega_a-\omega_m)},
\end{equation}
are in the second and third orders of the coupling strengths $g$ and $G$, respectively. $g_{\rm eff}$ achieves the maximum value when $\theta=\pi/4$. Then we stick to this sweet spot in the following.

\begin{figure}[http]
\centering
\includegraphics[width=0.45\linewidth]{e1p1m1Ut.eps}
\includegraphics[width=0.45\linewidth]{e1p1m1ut2.eps}
\includegraphics[width=0.45\linewidth]{e1p1m1ut3.eps}
\includegraphics[width=0.45\linewidth]{e1p1m1Ut4.eps}
\caption{Time evolution of the state population of the initial state $|e00\rangle$ and the target state $|g11\rangle$ under various coupling strengths to show the difference between the effective Hamiltonian and the full Hamiltonian at the avoided-level-crossing point shown in Fig.~\ref{e1p1m1eigen}. (a) $g=G=0.1\omega_a$, (b) $g=0.1\omega_a, G=0.05\omega_a$, (c) $g=0.05\omega_a, G=0.1\omega_a$, and (d) $g=G=0.05\omega_a$. }\label{e1p1m1ut}
\end{figure}

Under the effective Hamiltonian in Eq.~(\ref{Heff1}), a completed Rabi oscillation between the states $|e00\rangle$ and $|g11\rangle$ could be accurately established with a period of $\tau=\pi/|g_{\rm eff}|$ as demonstrated by the yellow dotted lines and the purple dash-dotted lines, respectively, in Fig.~\ref{e1p1m1ut}. We also plot the time evolutions under the full Hamiltonian $H$ in Eq.~(\ref{Hamiltonian}), where the blue solid lines and the red dashed lines represent the state-populations of $|e00\rangle$ and $|g11\rangle$, respectively. Note all the Rabi oscillations are ideal by the effective Hamiltonian, while the high-order effect emerges in the numerical evaluations by the full Hamiltonian.

Figures~\ref{e1p1m1ut}(a), (b), (c), and (d) vary with decreasing coupling strengths of $g$ and $G$ while fixing the other parameters. In the first period of the Rabi oscillation, one can observe that the maximum population of the state $|g11\rangle$ is about $P_{\rm max}=0.88$, $0.91$, $0.95$, and $0.98$, respectively under these four situations (see the red-dashed lines). The practical period of the Rabi oscillation between $|e00\rangle$ and $|g11\rangle$ becomes longer with smaller coupling strengths of $g$ and $G$ and it is slightly greater than that determined by the effective Hamiltonian. Similar to the maximum population, it is found that the generation time of the target state $\tau/2$ approaches also the ideal result by the effective Hamiltonian with decreasing coupling strengths. The relative errors of the period are about $2.7\%$, $4.3\%$, $0.4\%$, and $0.3\%$ for Figs.~\ref{e1p1m1ut}(a), (b), (c), and (d), respectively. Roughly, the reduction of $g$ has more impacts than that of $G$. In general, we can estimate from the four sub-figures that the analytical results from $H_{\rm eff}$ become gradually close to the numerical results from $H$ by decreasing the original coupling strengths $g$ and $G$. It is interesting to find that the range of validity of the effective Hamiltonian has approached the ultrastrong regime with $g=G=0.05\omega_a$.

\subsection{Generation protocol for Bell state}\label{Bellstep}

With the effective Hamiltonian in Eq.~(\ref{Heff1}), one can generate the Bell state of the photon-magnon subsystem by the following two-step protocol.

Step-$1$: The qubit frequency is adjusted to be far-off-resonant with both the resonator and the magnon modes (Note the latter two modes has already been set to be off-resonant). The whole hybrid system is prepared in the ground state of the free Hamiltonian, i.e., $|g\rangle_q|0\rangle_a|0\rangle_m=|g\rangle|00\rangle$. Then we perform a single-qubit gate operation $U$ on the qubit~\cite{onegate} and the operator reads,
\begin{equation}\label{U}
U=e^{i\frac{\pi}{4}\vec{\sigma}\cdot\vec{n}}=\begin{bmatrix}
1 & ie^{-i\phi}\\
ie^{i\phi} & 1
\end{bmatrix},
\end{equation}
where $\vec{\sigma}\equiv(\sigma_x, \sigma_y, \sigma_z)$ and $\vec{n}=(\cos\phi, \sin\phi, 0)$. It rotates the qubit into a superposed state
\begin{equation}\label{psi0s}
\frac{1}{\sqrt{2}}\left(|g\rangle+ie^{i\phi}|e\rangle\right),
\end{equation}
where the phase $\phi$ is tunable as desired and determines the final local phase of the double-excitation Bell states. For example, when $\phi=0,\pi$, the generated Bell state is $(|00\rangle\pm|11\rangle)/\sqrt{2}$.

Step-$2$: Thus the state of the whole system is now
\begin{equation}\label{psit0}
|\varphi(0)\rangle=\frac{1}{\sqrt{2}}\left(|g00\rangle+ie^{i\phi}|e00\rangle\right).
\end{equation}
Then we tune the qubit frequency to be nearly-resonant with the sum of the frequencies of the photon and magnon modes in an adiabatic way. As we shown in Eq.~(\ref{Heff1}), it will create an effective transition rate $g_{\rm eff}$ between the states $|e00\rangle$ and $|g11\rangle$, and the state $|g00\rangle$ is unaffected. The system state then evolves with time as
\begin{equation}\label{psit}
\begin{aligned}
|\varphi(t)\rangle&=\frac{1}{\sqrt{2}}\big[|g00\rangle+ie^{i\phi}\cos(g_{\rm eff}t)|e00\rangle\\
&+e^{i\phi}\sin(g_{\rm eff}t)|g11\rangle\big].
\end{aligned}
\end{equation}
After a time $T=\pi/(2|g_{\rm eff}|)$ (one half of the Rabi oscillation), the state evolves to
\begin{equation}\label{psiTT}
\begin{aligned}
|\varphi(T)\rangle&=\frac{1}{\sqrt{2}}\left(|g00\rangle+e^{i\phi}|g11\rangle\right)\\
&=|g\rangle\otimes\frac{1}{\sqrt{2}}\left(|00\rangle+e^{i\phi}|11\rangle\right).
\end{aligned}
\end{equation}
Then we tune the qubit-frequency faraway from the preceding point of the avoided level crossing and the Bell state $(|00\rangle+e^{i\phi}|11\rangle)/\sqrt{2}$ of the photon-magnon subsystem can be survival due to the large detuning between the resonator and the magnon mode.

\subsection{The fidelity of Bell state under dissipation}\label{Bellnum}

The fidelity of the generated state can be studied using the master equation approach by taking into account the dissipations from all parts of the hybrid system. By applying the standard Markovian approximation to the individual external environments (assumed to be at the vacuum states), we arrive at the Lindblad master equation for the density operator $\rho(t)$ of the hybrid system consisting of the qubit, the resonator and the magnon mode~\cite{bellstate},
\begin{equation}\label{lindblad}
\begin{aligned}
\dot{\rho}(t)&=-i[H_{\rm diag}, \rho]\\
&+\kappa_a\mathcal{L}[X_a]\rho(t)+\kappa_m\mathcal{L}[X_m]\rho(t)+\gamma\mathcal{L}[S_-]\rho(t).
\end{aligned}
\end{equation}
Here $H_{\rm diag}$ indicates that the full Hamiltonian $H$ is now expressed by its eigenvectors $|E_n\rangle$'s. $\kappa_a$, $\kappa_m$, and $\gamma$ are the dissipation rates for the resonator, the magnon, and the qubit, respectively. And the superoperator $\mathcal{L}[O]$ is defined as
\begin{equation}\label{L}
\mathcal{L}[O]\rho\equiv\frac{1}{2}\left(2O\rho O^\dag-O^\dag O\rho-\rho O^\dag O\right).
\end{equation}
Here $O=X_a, X_m, S_-$ are the dressed lowering operators, defined respectively in terms of their bare counterparts $o=a, m, \sigma_-$ as
\begin{equation}\label{O}
O\equiv\sum_{E_n>E_m}\langle E_m|(o+o^{\dag})|E_n\rangle|E_m\rangle\langle E_n|.
\end{equation}
To simplify the robustness estimation of our protocol but with no loss of generality, we assume all of the decoherence rates to be in the same order of magnitude $\kappa_a=\kappa_m=\gamma=\kappa$. It is consistent with the relative decoherence rates obtained in recent experiments~\cite{magnonqubit,magnonqubit2,ultrastrong}, i.e., $\kappa_a/\omega_a\sim 10^{-6}-10^{-5}$, $\kappa_m/\omega_m\sim 10^{-6}-10^{-5}$, and $\gamma/\omega_q\sim 10^{-6}-10^{-5}$.

\begin{figure}[htbp]
\centering
\includegraphics[width=0.4\textwidth]{Bellfidelity1.eps}
\includegraphics[width=0.4\textwidth]{Bellfidelity2.eps}
\caption{The final fidelity $F(T)$ in the parametric space of the coupling strengths $g$ and $G$ under (a) $\kappa_a=\kappa_m=\gamma=0$ and (b) $\kappa_a=\kappa_m=\gamma=10^{-5}\omega_a$. The other parameters are fixed as $\omega_m=1.7\omega_a$ and $\theta=\pi/4$.}\label{e1p1m1master}
\end{figure}

The state-generation fidelity $F$ is defined as $F(t)=\langle\varphi(t)|\rho|\varphi(t)\rangle$, where $|\varphi(T)\rangle$ is the target state. Here the phase $\phi$ is set as zero, then $|\varphi(T)\rangle=(|00\rangle+|11\rangle)/\sqrt{2}$.

In Fig.~\ref{e1p1m1master}(a) and (b), we demonstrate and compare the generation fidelities $F(T)$ of the Bell state with no decoherence (under the full Hamiltonian rather than the effective Hamiltonian) and that with dissipation (under the Lindblad master equation). To be consistent with the dynamics in Fig.~\ref{e1p1m1ut}, it is shown in Fig.~\ref{e1p1m1master}(a) that the generation fidelity is enhanced when the coupling strengths $g$ and $G$ are reduced. One can see that the fidelity approaches $0.95$ when both the coupling strengths $g$ and $G$ are reduced to about $0.05\omega_a$. This observation supports again the validity of our effective Hamiltonian. In contrast, a smaller coupling strength does not always give rise to a higher fidelity. As shown in Fig.~\ref{e1p1m1master}(b), it is found that the generation fidelity is about $0.77$ under the coupling strengths $g=G=0.05\omega_a$, and about $0.90$ when $g=G=0.1\omega_a$. A compromise in terms of the coupling strength in fidelity is expected due to the fact that the period of the desired Rabi oscillation is inversely proportional to the coupling strengths and the dissipation becomes more destructive under a longer evolution time. A complete scanning over the parametric space demonstrates a remarkable working regime for the generation of the Bell state: $0.07<g/\omega_a<0.13$ and $0.05<G/\omega_a<0.1$, where we have $F>0.9$.

\begin{figure}[htbp]
\centering
\includegraphics[width=0.45\textwidth]{bellfidelity.eps}
\caption{Time evolution of the fidelity by the master equation~(\ref{lindblad}) under different dissipation rates. The other parameters are fixed as $\omega_m=1.7\omega_a$, $G=g=0.1\omega_a$, and $\theta=\pi/4$.}\label{e1p1m1lindblad}
\end{figure}

We then pick up a pair of $g$ and $G$ to plot the fidelity dynamics under different dissipation rates in Fig.~\ref{e1p1m1lindblad}. One can observe that the protocol works well for $\kappa\le10^{-5}\omega_a$, producing the desired Bell state with a fidelity over $0.90$, close to $0.93$ in the situation with no decoherence. The fidelity $F$ can be maintained above $0.78$ even when $\kappa$ is enhanced to $10^{-4}\omega_a$. It means that our protocol is robust even when all the decoherence channels of the hybrid system are simultaneously switched on.

\section{Generating GHZ state of qubit-photon-magnon system}

\subsection{Effective Hamiltonian}\label{GHZ}

This section is devoted to generating the GHZ state of the whole hybrid system, which shares the same key step or basic mechanism with the protocol for generating the Bell state in Sec.~\ref{bell}. At the desired point of the avoided level crossing, one can prepare the ``excited'' state $|e11\rangle\equiv|e\rangle_q|1\rangle_a|1\rangle_m$ from $|g20\rangle$, where the annihilation of one photon excites the qubit and the magnon mode simultaneously.

\begin{figure}[htbp]
\centering
\includegraphics[width=0.45\textwidth]{GHZeigen.eps}
\caption{Energy levels and the avoided level crossing for $|g20\rangle\leftrightarrow|e11\rangle$, normalized by $\omega_m$ and plotted as functions of the qubit transition frequency $\omega_q$. Here the parameters are fixed as $\omega_a=2.4\omega_m$, $g=G=0.1\omega_m$, and $\theta=\pi/4$.}\label{GHZeigen}
\end{figure}

In Fig.~\ref{GHZeigen}, the avoided level crossing for $|g20\rangle$ and $|e11\rangle$ is distinguished in the dark circle, demonstrating an effective transition rate $|g'_{\rm eff}|$. The energy shift $\Delta$ between the avoided-level-crossing point and the exact double-resonant point arises from the interaction Hamiltonian $V$ and then defined by $\omega_q=\omega_a-\omega_m+\Delta$. According to Appendix~\ref{appb}, the effective Hamiltonian reads,
\begin{equation}\label{Heff2}
H'_{\rm eff}=g'_{\rm eff}\left(|e11\rangle\langle g20|+|g20\rangle\langle e11|\right),
\end{equation}
with the effective coupling strength
\begin{equation}\label{geffghz}
g'_{\rm eff}=-\frac{\sqrt{2}G^2g\sin(2\theta)(\omega_a+3\omega_m)}{\omega_a\omega_m(\omega_a+\omega_m)}.
\end{equation}
And the energy shift reads
\begin{equation}\label{deltaghz}
\Delta=4G^2\cos^2\theta\left(\frac{1}{\omega_m}-\frac{1}{2\omega_a-\omega_m}\right)+\frac{2g^2}{\omega_a-\omega_m}.
\end{equation}

\begin{figure}[htbp]
\centering
\includegraphics[width=0.45\linewidth]{GHZut.eps}
\includegraphics[width=0.45\linewidth]{GHZut2.eps}
\includegraphics[width=0.45\linewidth]{GHZut3.eps}
\includegraphics[width=0.45\linewidth]{GHZut4.eps}
\caption{Time evolution of the state population of the initial state $|g20\rangle$ and the target state $|e11\rangle$ under various coupling strengths to show the difference between the effective Hamiltonian and the full Hamiltonian at the avoided-level-crossing point shown in Fig.~\ref{GHZeigen}. (a) $g=G=0.1\omega_m$, (b) $g=0.1\omega_m, G=0.05\omega_m$, (c) $g=0.05\omega_m, G=0.1\omega_m$, and (d) $g=G=0.05\omega_m$. }\label{p1e1m1ut}
\end{figure}

Under the effective Hamiltonian in Eq.~(\ref{Heff2}), a completed Rabi oscillation between states $|g20\rangle$ and $|e11\rangle$ could be accurately established with a period of $\tau=\pi/|g'_{\rm eff}|$ as demonstrated by the yellow dotted lines and the purple dash-dotted lines, respectively in Fig.~\ref{p1e1m1ut}. To estimate the range of validity of  the effective Hamiltonian, we also present the time evolutions under the full Hamiltonian $H$ in Eq.~(\ref{Hamiltonian}), where the blue solid lines and the red dashed lines indicate the states populations for $|g20\rangle$ and $|e11\rangle$, respectively. In Figs.~\ref{p1e1m1ut}(a), (b), (c), and (d), one can observe that the maximum population of the state $|e11\rangle$ in the first period of Rabi oscillation is gradually enhanced with the decreasing coupling strengths of $g$ and $G$ while fixing the other parameters, i.e., $P_{\rm max}=0.94$, $0.95$, $0.96$, and $0.98$, respectively. While the relative errors of the period between the effective Hamiltonian and the total Hamiltonian have no clear dependence on the coupling strengths. It indicates that our estimation over the effective coupling rate $|g'_{\rm eff}|$ for the transition $|e11\rangle\leftrightarrow|g20\rangle$ in Fig.~\ref{GHZeigen} is not as accurate as that for the transition $|e00\rangle\leftrightarrow|g11\rangle$ in Fig.~\ref{e1p1m1eigen}.

\subsection{Generation protocol for GHZ state}\label{GHZstep}

With the effective Hamiltonian in Eq.~(\ref{Heff2}), one can generate the GHZ state for the whole qubit-photon-magnon system by the following three-step protocol.

Step-$1$: The transition frequency of the superconducting qubit is tuned to be far-off-resonant with those for both photon and magnon modes. And the whole system is initially at the ground state of the free Hamiltonian, i.e., in the state $|g\rangle_q|0\rangle_a|0\rangle_m=|g\rangle|00\rangle$. Then we rotate the qubit state into a superposed state
\begin{equation}\label{psi0sGHZ}
\frac{1}{\sqrt{2}}\left(|g\rangle-e^{i\phi}|e\rangle\right),
\end{equation}
with a single-qubit gate
\begin{equation}\label{UGHZ}
U=e^{i\frac{\pi}{4}\vec{\sigma}\cdot\vec{n}}=\begin{bmatrix}
1 & -e^{-i\phi}\\
-e^{i\phi} & 1
\end{bmatrix},
\end{equation}
where $\vec{n}=(\sin\phi, \cos\phi, 0)$. The phase $\phi$ is tunable as desired and determines the final local phase in the GHZ states $(|g00\rangle+e^{i\phi}|e11\rangle)/\sqrt{2}$.

Step-$2$: The state of the whole system is thus written as
\begin{equation}\label{psi0s}
|\varphi'(0)\rangle=\frac{1}{\sqrt{2}}\left(|g00\rangle-e^{i\phi}|e00\rangle\right).
\end{equation}
Then we tune the qubit frequency adiabatically into the nearly-two-photon resonance with the resonator mode, namely, $\omega_q=2\omega_a+\tilde{\delta}$. Note the magnon mode is far off-resonant from them. In this case, the full Hamiltonian of the system in Eq.~(\ref{Hamiltonian}) turns out to describe a two-photon Jaynes-Cummings model~\cite{bellstate},
\begin{equation}\label{Hefftwophoton}
\tilde{H}_{\rm eff}=\tilde{g}_{\rm eff}\left(|e00\rangle\langle g20|+|g20\rangle\langle e00|\right),
\end{equation}
where $\tilde{g}_{\rm eff}=-\sqrt{2}G^2\sin(2\theta)/\omega_a$ and $\tilde{\delta}=-2G^2/\omega_a-2g^2/(\omega_a+\omega_m)+2g^2/(\omega_a-\omega_m)$ could also be directly obtained by the high-order Fermi's Golden rule given in Eq.~(\ref{Fermi}) or the standard perturbation method in Eq.~(\ref{secondp}). The state $|g00\rangle$ is not influenced by $\tilde{H}_{\rm eff}$, and the state $|e00\rangle$ evolves with time as
\begin{equation}\label{psitau}
\begin{aligned}
|\varphi'(t)\rangle&=\frac{1}{\sqrt{2}}\big[|g00\rangle-e^{i\phi}\cos(\tilde{g}_{\rm eff}t)|e00\rangle\\
&+ie^{i\phi}\sin(\tilde{g}_{\rm eff}t)|g20\rangle\big].
\end{aligned}
\end{equation}
After a time $T'=\pi/|2\tilde{g}_{\rm eff}|$, the state becomes
\begin{equation}\label{psitaus}
|\varphi'(T')\rangle=\frac{1}{\sqrt{2}}\left(|g00\rangle+ie^{i\phi}|g20\rangle\right).
\end{equation}

Step-$3$: We then tune the qubit frequency into the near-resonance with the detuning between the photon and magnon modes, i.e., the avoided-level-crossing point shown in Fig.~\ref{GHZeigen}: $\omega_q=\omega_a-\omega_m+\Delta$. Then driven by Eq.~(\ref{Heff2}), the state of the hybrid system evolves as
\begin{equation}\label{psit}
\begin{aligned}
|\varphi'(T'+t)\rangle&=\frac{1}{\sqrt{2}}\big[|g00\rangle+ie^{i\phi}\cos(g'_{\rm eff}t)|g20\rangle\\
&+e^{i\phi}\sin(g'_{\rm eff}t)|e11\rangle\big].
\end{aligned}
\end{equation}
After a time $T=\pi/(2|g'_{\rm eff}|)$, it turns out to be
\begin{equation}\label{psiTT}
|\varphi'(T'+T)\rangle=\frac{1}{\sqrt{2}}\left(|g00\rangle+e^{i\phi}|e11\rangle\right).
\end{equation}
Then we tune the qubit faraway from the resonant point and the GHZ state $(|g00\rangle+e^{i\phi}|e11\rangle)/\sqrt{2}$ can be maintained for the transition frequencies of the three subsystems are now far-off-detuning from each other.

\subsection{The fidelity of GHZ state under dissipation}\label{GHZnum}

The generation fidelity of the GHZ state can be also studied using the standard Lindblad master equation~(\ref{lindblad}). The fidelity is defined as $F(t)=\langle\varphi(T'+t)|\rho|\varphi(T'+t)\rangle$ with $|\varphi(T'+T)\rangle=(|g00\rangle+|e11\rangle)/\sqrt{2}$ in Eq.~(\ref{psiTT}) and $\phi=0$.

\begin{figure}[htbp]
\centering
\includegraphics[width=0.4\textwidth]{GHZfidelity.eps}
\includegraphics[width=0.4\textwidth]{GHZfidelity2.eps}
\caption{The final fidelity $F(T)$ in the parametric space of the coupling strengths $g$ and $G$ under (a) $\kappa_a=\kappa_m=\gamma=0$ and (b) $\kappa_a=\kappa_m=\gamma=10^{-5}\omega_m$. The other parameters are fixed as $\omega_a=2.4\omega_m$ and $\theta=\pi/4$.}\label{fidelity}
\end{figure}

In Fig.~\ref{fidelity}(a) and (b), we demonstrate and compare the final fidelity $F(T)$ of the GHZ state with no decoherence and that under the simultaneous dissipations from qubit, resonator, and magnon. In Fig.~\ref{fidelity}(a), one can observe that a high-fidelity of state-generation can be maintained when the coupling strengths $g$ and $G$ are reduced, being consistent with the Rabi oscillations in Fig.~\ref{p1e1m1ut}. The fidelity is greater than $0.95$ for the coupling strengths $g, G<0.2\omega_m$ and is even close to $0.98$ when $g=0.05\omega_m$ and $G=0.05\omega_m$. This justifies our effective Hamiltonian in Eq.~(\ref{Heff2}). In the presence of the dissipation, however, the dependence of the fidelity on the coupling strengths in Fig.~\ref{fidelity}(b) is not monotonic as shown in Fig.~\ref{fidelity}(a). For example, when $g=G=0.05\omega_m$, the generation fidelity is about $0.70$; in contrast, when $g=G=0.1\omega_m$, it is over $0.95$. During the long-time evolution induced by small coupling strengths, the effect from decoherence on the fidelity will destroy the fidelity of the final GHZ state. When $0.1\le g/\omega_m\le0.17$ and $0.1\le G/\omega_m\le0.17$, we have an optimized regime with $F>0.9$.

\begin{figure}[htbp]
\centering
\includegraphics[width=0.45\textwidth]{GHZlindblad.eps}
\caption{Time evolution of the fidelity of the hybrid GHZ state under the master equation~(\ref{lindblad}) under different decoherence rates. Here the parameters are fixed as $\omega_a=2.4\omega_m$, $G=g=0.1\omega_m$ and $\theta=\pi/4$.}\label{GHZlindblad}
\end{figure}

In Fig.~\ref{GHZlindblad}, we present the dynamics of the generation fidelity with a fixed pair of $g$ and $G$. One can observe that the protocol for the hybrid system works well for $\kappa\le10^{-5}\omega_m$, producing a global GHZ state with a fidelity over $0.95$, close to $0.98$ in the ideal situation with no decoherence. It decreases to $0.62$ when $\kappa$ is enhanced to $10^{-4}\omega_m$. The global GHZ state is therefore less robust than the local Bell state.

\section{Discussion}\label{discussion}

After a minor modification, the protocol for GHZ state in Sec.~\ref{GHZ} can be immediately applied to generate the double-excitation Bell state of the qubit-magnon subsystem. It is still of a three-step scheme. Step-$1$ remains invariant, i.e., the whole system again starts from Eq.~(\ref{psi0sGHZ}):
\begin{equation}
|\varphi'(0)\rangle=\frac{1}{\sqrt{2}}\left(|g\rangle_q-e^{i\phi}|e\rangle_q\right)|0\rangle_a|0\rangle_m.
\end{equation}

Step-$2$: We now tune the qubit frequency to be single-photon resonant with the resonator mode instead of the double-photon resonance as in Sec.~\ref{GHZstep}. Around this point, the ground state $|g00\rangle$ holds and the state $|e00\rangle$ undergoes a Rabi oscillation with an effective transition rate $G$. Then after an evolution time $T'=\pi/(2G)$, the system state evolves into
\begin{equation}\label{psitaus}
|\varphi'(T')\rangle=\frac{1}{\sqrt{2}}\left(|g00\rangle+ie^{i\phi}|g10\rangle\right).
\end{equation}

Step-$3$: We then again tune the qubit frequency to satisfy $\omega_q=\omega_a-\omega_m+\Delta$, which is the avoided-level-crossing point for the states $|g10\rangle$ and $|e01\rangle$. The effective Hamiltonian becomes $g'_{\rm eff}(|g10\rangle\rangle e01|+|e01\rangle\rangle g10|)$. Although formally it seems almost the same as Eq.~(\ref{Heff2}), here $g'_{\rm eff}$ is yet not the same as Eq.~(\ref{geffghz}) since the transition paths as well as the transition rates connecting $|g10\rangle$ and $|e01\rangle$ are not the same as those connecting $|g20\rangle$ and $|e11\rangle$. In the current case, it is found that
\begin{equation}\label{deltabell2}
g'_{\rm eff}=-\frac{2G^2g\sin(2\theta)}{\omega_m(\omega_a+\omega_m)}.
\end{equation}
And after a time $T=\pi/(2|g'_{\rm eff}|)$, the state in Eq.~(\ref{psitaus}) becomes
\begin{equation}\label{psibell2}
|\varphi'(T'+T)\rangle=\frac{1}{\sqrt{2}}\left(|g00\rangle+e^{i\phi}|e01\rangle\right),
\end{equation}
which is the desired Bell state for the qubit-magnon subsystem, since the middle state for the resonator is now separable. At this moment, one can tune the qubit frequency faraway from the resonant point and then the Bell state $(|g0\rangle+e^{i\phi}|e1\rangle)/\sqrt{2}$ can be maintained.

\section{Conclusion}\label{conclusion}

The protocols for generating local and global entangled states we proposed can be performed in a hybrid setup consisting of a single superconducting qubit, a microwave resonator, and a YIG sphere (magnon)~\cite{magnonqubit,magnonqubit2}. The resonator is simultaneously strongly coupled with the magnon via the magnetic dipole interaction, and with the qubit via a general Rabi interaction. In recent experiments~\cite{magnonqubit,magnonqubit2,ultrastrong}, the coupling strength between photon and qubit $G/2\pi\approx120$ MHz, the coupling strength between photon and magnon $g/2\pi\approx20$ MHz, and the transition frequencies of photon mode, magnon mode and qubit are almost in the same order of GHz. Thus the generation time $T$ of our protocols is nearly about $0.1\sim10\mu s$. Note our target entangled state is of a ``discrete-variable'' type rather than a ``continuous-variable'' one in Ref.~\cite{mppentang}. Our study is of interested in pursuit of the entangled states with the counter-rotating interaction and of importance to control the quantum state in a level-resolved hybrid system.

In conclusion, we have presented a concise protocol for the deterministic generation of local Bell state of the photon-magnon or the qubit-magnon subsystems, and global GHZ state of the whole qubit-photon-magnon system. Our protocol relies on the effective Hamiltonian at the avoided-level-crossing points, which reserves the effects of the counter-rotating interactions and the leading-order contributions of the state transitions. By properly tuning the transition frequency of the superconducting qubit, various scenarios of ``three-wave-mixing'' relevant to either double-excitation Bell state or GHZ state are constructed. Moreover, the generation fidelities of these entangled states are numerically estimated with the standard Lindblad master equation and our protocol is found to be robust against the external dissipation noises.

\section*{Acknowledgments}

We acknowledge grant support from the National Science Foundation of China (Grants No. 11974311 and No. U1801661), and the Zhejiang Provincial Natural Science Foundation of China under Grant No. LD18A040001.

