Kerr-magnon-assisted asymptotic stationary photon-phonon squeezing

Hybrid quantum systems consisting of collective spin excitations in ferromagnetic crystals have recently attracted intensive attention [1, 2, 3, 4, 5, 6, 7, 8]. They offer promising avenues for advancements in quantum computing [9], quantum communication [10], and quantum sensing [11]. Similar to cavity-QED [12] and cavity optomechanics [13], cavity magnomechanics [14, 15, 16, 17] developed rapidly as an alternative candidate for quantum information processing in both theoretical and experimental aspects. In particular, a cavity magnomechanical system comprises a single-crystal yttrium iron garnet (YIG) sphere inside a microwave cavity. The magnon mode arises from excitations of the collective angular momentum within the magnetic-material sphere. It couples with the cavity photon through magnetic dipole interaction and the sphere deformation phonon mode via magnetostrictive force. Typical applications based on such hybrid system include quantum entanglement [18, 19] and steering [20, 21], quantum squeezed states [22, 23, 24], and quantum memory [25, 26, 27].

Two-mode squeezed states (TMSS), also named Einstein-Podolsky-Rosen (EPR) states, are crucial in quantum computation [28], information [29], teleportation [30], and metrology [31]. Bosonic TMSS can be generated by mixing two single-mode squeezed states on a beam splitter [32] or via a nonlinear interaction [33, 34] such as spontaneous parametric down conversion [35]. A nondegenerate optical parametric amplifier is often used to generate optical TMSS [36, 37, 38]. Cavity optomechanics [13] provides an alternative model for creating optical [39, 40, 41] or mechanical [42, 43] TMSS. In addition, TMSS are well established in thermal gases [44], Bose-Einstein condensates of ultracold atoms [45, 46, 47, 48], spin ensembles in cavities [49, 50, 51, 52], antiferromagnet magnons [53], and superconducting circuits based on Josephson junctions [54, 55].

In this work, we propose an approach to generate photon-phonon TMSS in the cavity magnomechanics [14, 15, 16, 17] on account of the fundamental interest in a level-resolved process. The squeezing generation is governed by an effective Hamiltonian that describes photon-phonon squeezing interaction, assisted by exploiting strong or even ultrastrong magnon-photon and magnon-phonon interactions. To confirm the validity of the effective Hamiltonian, we employ the diagonalization of the Liouvillian superoperator of the whole system. This approach can effectively address the squeezing Hamiltonian that does not conserve the whole excitations. It is thus distinct from the previous method [26, 56] involving a standard numerical diagonalization of the system Hamiltonian. Additionally, our approach can be extended to arbitrary bosonic systems, such as cavity optomechanics [13] and cavity optomagnomechanics [57], to evaluate the two-mode squeezing induced by virtual processes.

Such a two-mode squeezing naturally leads to entanglement without reservoir engineering. Under the constraint of the stable condition, the squeezing level cannot go beyond $3$ ${\rm dB}$ below the vacuum limit [39]. However, by examining the system’s dynamics within the open-quantum-system framework, we find that the stability of the Gaussian system, i.e., the covariance matrix (CM) becomes invariant when $t\to\infty$, is a sufficient but not necessary condition for the stationary generation of TMSS. We find that an asymptotic stationary TMSS can be obtained in unstable evolutions, displaying an enhanced squeezing level exceeding the steady limit. Environmental noises alter the optimized quadrature operator of two-mode squeezing while simultaneously stabilizing TMSS in an asymptotic way. The coupling between two modes can change the squeezing level, but it does not influence the squeezing stationarity, even if it is beyond the stability threshold [35] of the CM. Through analysis of the system’s dynamics under noise, we found that strong coupling essentially implies strong and stationary quantum entanglement.

The rest of this work is organized as follows. In Sec. II, we introduce a hybrid cavity magnomechanical system and provide an effective Hamiltonian for photon-phonon squeezing mediated by the magnon. In Sec. III, we confirm the effective Hamiltonian by comparing the effective coupling strength and energy shift with numerical results obtained by diagonalization of the system’s Liouvillian superoperator. Section IV phenomenologically analyzes the generation process of the photon-phonon TMSS by the quantum Langevin equation. We find that the asymptotic stationary two-mode squeezing can be obtained even in an unstable dynamic regime. Finally, we discuss the experimental feasibility and summarize the work in Sec. V.

Consider a hybrid cavity magnomechanical system as shown in Fig. 1, where a YIG sphere is inserted into a three-dimensional microwave cavity. The system is constituted by the microwave-mode photons, the magnons provided by the YIG sphere, and the vibrational modes (phonons) of the sphere. The magnons are coupled to photons via the Zeeman interaction and to phonons by the magnetostrictive interaction. The hybrid system has been experimentally realized in recent works [14, 15, 16]. The full system Hamiltonian thus reads ($\hbar\equiv 1$)

	$\displaystyle H_{s}$	$\displaystyle=\omega_{a}a^{\dagger}a+\omega_{m}m^{\dagger}m-K_{m}m^{\dagger}mm^{\dagger}m+\omega_{b}b^{\dagger}b$		(1)
		$\displaystyle+g_{ma}(a^{\dagger}m+am^{\dagger})+g_{mb}m^{\dagger}m(b+b^{\dagger})$
		$\displaystyle+\epsilon_{d}(a^{\dagger}e^{-i\omega_{d}t}+ae^{i\omega_{d}t}),$

where $a(a^{\dagger})$, $m(m^{\dagger})$, and $b(b^{\dagger})$ are the annihilation (creation) operators of the photon mode, the magnon of the ground Kittel mode [1], and the phonon mode with transition frequencies $\omega_{a}$, $\omega_{m}$, and $\omega_{b}$, respectively. The magnon-mode frequency is determined by $\omega_{m}=\gamma h$, where $\gamma$ is the gyromagnetic ratio and $h$ is the external bias magnetic field. Thus, it can be appropriately tuned by the external magnetic field. $K_{m}$ is the nonlinear coefficient for the Kerr effect due to the ensued magnetocrystalline anisotropy, which can be either positive or negative by adjusting the crystallographic axis of the YIG sphere along the bias magnetic field and is inversely proportional to the volume of the YIG sphere. The Kerr effect cannot be neglected when the magnon excitation number is sufficiently large [15, 58, 8, 3]. $g_{ma}$ is the photon-magnon coupling strength, entering into the strong-coupling regime. The single-magnon magnomechanical coupling strength $g_{mb}$ is typically small, considering the large frequency mismatch between the magnon and the phonon modes, yet it can be compensated by a strong drive. The last term in Eq. (1) describes the external driving of the photon mode, where $\epsilon_{d}$ is the Rabi frequency and $\omega_{d}$ is the driving frequency.

The magnon mode under strong driving is assumed to have a large expectation value $|\langle m\rangle|\gg 1$, which allows us to linearize the system dynamics. Following the standard linearization approach [13, 18], the full system Hamiltonian turns out to be

	$\displaystyle H$	$\displaystyle=\Delta_{a}a^{\dagger}a+\Delta^{\prime}_{m}m^{\dagger}m+\omega_{b}b^{\dagger}b+g\cosh r(am^{\dagger}+a^{\dagger}m)$		(2)
		$\displaystyle+g\sinh r(ame^{-i\theta}+a^{\dagger}m^{\dagger}e^{i\theta})$
		$\displaystyle+Ge^{r}(me^{-i\frac{\theta}{2}}+m^{\dagger}e^{i\frac{\theta}{2}})(b+b^{\dagger}),$

where $\Delta^{\prime}_{m}=\Delta_{m}/\cosh(2r)$ is the modified magnon detuning, $\Delta_{m}=\omega_{m}-\omega_{d}-2|K|$ and $K=K_{m}\langle m\rangle^{2}$ is the driving-enhanced Kerr parameter. $r$ is the squeezing parameter induced by the linearization of the Kerr effects and $\tanh(2r)=2|K|/\Delta_{m}$. $g=g_{ma}$ for simplicity, $G=g_{mb}|\langle m\rangle|$ is the driving-enhanced magnomechanical coupling strength. $\theta$ is a phase associated with $\langle m\rangle=|\langle m\rangle|e^{i\theta/2}$. The details can be found in Appendix A.

At the large detuning regime, i.e., $g\cosh r,g\sinh r,Ge^{r}\ll|\Delta^{\prime}_{m}-\omega_{b}|,|\Delta^{\prime}_{m}-\Delta_{a}|$, and under the near-resonant condition $\Delta_{a}=-\omega_{b}+\delta$, we can extract an effective Hamiltonian describing the photon-phonon squeezing by perturbative theory [56, 59]. The effective Hamiltonian is found to be

$$H_{\rm eff}=g_{\rm eff}(e^{i\frac{\theta}{2}}a^{\dagger}b^{\dagger}+e^{-i\frac{\theta}{2}}ab),$$

(3)

where the effective coupling strength is

$$g_{\rm eff}=gG\cosh(2r)\frac{\omega_{b}\cosh(2r)-\Delta_{m}e^{2r}}{\Delta^{2}_{m}-\omega^{2}_{b}\cosh^{2}(2r)},$$

(4)

and the energy shift is

$$\delta=\frac{2G^{2}\Delta_{m}e^{2r}\cosh(2r)+g^{2}(\Delta_{m}-\omega_{b})\cosh^{2}(2r)}{\Delta^{2}_{m}-\omega^{2}_{b}\cosh^{2}(2r)}.$$

(5)

The details can be found in Appendix B. The phase $\theta$ specifies the squeezing quadrature operator but does not affect the squeezing level. Therefore, we set $\theta=\pi$ in the following for simplicity and with no loss of generality. The corresponding squeezing operators can be written as

$$X(t)=\frac{1}{\sqrt{2}}[X_{a}(t)+X_{b}(t)],$$

(6)

where $X_{o}=(o+o^{\dagger})/\sqrt{2},o=a,b$. When the initial state is a vacuum state, its variance turns into $\Delta X(t)=\langle X^{2}(t)\rangle-\langle X(t)\rangle^{2}=e^{2g_{\rm eff}t}/2$ under the time-evolution of Hamiltonian (3) [35, 43]. Obviously, the quadrature operator $X$ is squeezed, i.e., $\Delta X(t)<1/2$, when $g_{\rm eff}<0$.

In this section, we check the applicability range of the effective Hamiltonian in Eq. (3) in terms of the coupling strengths and squeezing parameters. In previous works [56, 26], diagonalizing the full system Hamiltonian in a truncated finite-dimensional Hilbert space is employed to confirm the effective Hamiltonian constructed by virtual transitions. It is essential to observe a desired avoided level crossing between two eigenstates in the eigenenergies diagram, and the energy splitting at the avoided level crossing point precisely equals twice the effective coupling strength. However, the two-mode squeezing effective Hamiltonian (3) is not conserved in the excitation number, as shown by the non-commutativity $[H_{\rm eff},\hat{N}]\neq 0$ with $\hat{N}=a^{\dagger}a+b^{\dagger}b$ the excitation-number operator. Thus, the effective Hamiltonian (3) cannot be rigorously diagonalized within an appropriate truncated Hilbert space. Additionally, an explicit avoided-level crossing between two eigenstates does not appear for the two-mode squeezing Hamiltonian. Consequently, we propose a distinct method to validate the two-mode squeezing Hamiltonian in Eq. (3).

Under the evolution of the full system Hamiltonian in Eq. (2), the time-evolved quadrature operators in the Heisenberg picture can be written as

$$\dot{u}(t)=i[H,u(t)]=i\mathcal{L}u(t),$$

(7)

where $u(t)=[X_{a}(t),Y_{a}(t),X_{b}(t),Y_{b}(t),X_{m}(t),Y_{m}(t)]^{T}$ is the vector of quadrature operators, and $X_{o}=(o+o^{\dagger})/\sqrt{2},Y_{o}=(o-o^{\dagger})/i\sqrt{2},o=a,b,m$. $\mathcal{L}$ represents the Liouvillian superoperator,

$\displaystyle\mathcal{L}=-i\begin{bmatrix}0&\Delta_{a}&0&0&0&g_{+}\\ -\Delta_{a}&0&0&0&g_{-}&0\\ 0&0&0&\omega_{b}&0&0\\ 0&0&-\omega_{b}&0&0&-G_{r}\\ 0&g_{+}&G_{r}&0&0&\Delta^{\prime}_{m}\\ g_{-}&0&0&0&-\Delta^{\prime}_{m}&0\end{bmatrix},$

(8)

where $g_{\pm}=g\sinh r\pm g\cosh r$ and $G_{r}=2Ge^{r}$. The Heisenberg equation in Eq. (7) can be regarded as a discrete Schrödinger equation, where $u(t)$ is conceptualized as an effective operator wave function [60]. The superoperator $\mathcal{L}$ then can be analogously regarded as the full system Hamiltonian, and its diagonalization values are the system’s eigenvalues.

We now analyze the distinct phenomenon observed in the energy diagram of the Liouvillian superoperator at two-mode squeezing. Rotating the effective Hamiltonian (3) with $\theta=\pi$ into the lab frame, it becomes

$$H_{\rm ab}=\Delta_{a}a^{\dagger}a+\omega_{b}b^{\dagger}b+g_{\rm eff}(ia^{\dagger}b^{\dagger}-iab).$$

(9)

The corresponding Heisenberg equation is

$$\dot{u}^{\rm eff}(t)=i[H_{\rm ab},u^{\rm eff}(t)]=i\mathcal{L}_{\rm ab}u^{\rm eff}(t),$$

(10)

where $u^{\rm eff}(t)=[X_{a}(t),Y_{a}(t),X_{b}(t),Y_{b}(t)]^{T}$. Four eigenvalues of the Liouvbilian superoperator $\mathcal{L}_{\rm ab}$ can be derived as

	$\displaystyle E_{\pm}$	$\displaystyle=\frac{\omega_{b}-\Delta_{a}\pm\sqrt{(\omega_{b}+\Delta_{a})^{2}-4g_{\rm eff}^{2}}}{2},$		(11)
	$\displaystyle E^{\prime}_{\pm}$	$\displaystyle=-E_{\mp}.$

The real parts of the eigenvalues $E_{\pm}$ ($E^{\prime}_{\pm}$) converge, while the imaginary parts split as the detuning $\Delta_{a}$ gradually approaches $-\omega_{b}$. Until $\Delta_{a}=-\omega_{b}$, the real parts of $E_{\pm}$ ($E^{\prime}_{\pm}$) become identical, while the imaginary parts reach their extreme values of $\pm g_{\rm eff}$. Then, in the energy-level diagram of the whole superoperator $\mathcal{L}$ as a function of $\Delta_{a}$, one can demonstrate the two-mode squeezing interaction through the level attractions of the real parts and the maximal splittings of the imaginary parts.

We plot the energy levels (all six real and two relevant imaginary parts) in Figs. 2(a) and 2(b), where the eigenvalues $\{E_{n}\}$ are obtained by the standard numerical diagonalization on the whole superoperator $\mathcal{L}$ in Eq. (8). Figure 2(a) shows the real parts of all six eigenvalues. The light-blue and dark-blue lines describe the energies of the magnon and are not relevant to photon-phonon squeezing. For the other four eigenvalues, two level attractions appear simultaneously as the detuning frequency of photon $\Delta_{a}$ approaches (but does not exactly equal) the opposite frequency of phonon $-\omega_{b}$. This level attraction is highlighted by a dark circle, and the inset further emphasizes it. The imaginary parts of the two relevant eigenvalues (red and orange lines) are illustrated in Fig. 2(b). As the real parts of the two eigenvalues gradually converge, their imaginary parts progressively increase, reaching a maximum absolute value $|g_{\rm eff}|$ at $\Delta_{a}=-\omega_{b}+\delta$. The shift $\delta$ is induced by the mutual interaction between the photon (phonon) and the magnon.

The maximal splitting $|g_{\rm eff}|$ of the imaginary parts of the two eigenvalues [see Fig. 2(b)] is presented in Figs. 3(a) and 3(c) as a function of the original coupling strengths $g$ and $G$ in Eq. (2), respectively. The analytical result in Eq. (4) is compared to the numerical simulation over the superoperator $\mathcal{L}$ in Eq. (8). Blue dots and purple squares represent the numerical results at squeezing parameters $r=0$ and $r=0.25$, respectively. The red solid and black dashed lines are the analytical results at $r=0$ and $r=0.25$, respectively. In Fig. 3(a), the analytical $g_{\rm eff}$ do match well with their numerical results for the coupling strength $g\leq 0.3\omega_{b}$, regardless of whether $r$ is $0$ or $0.25$. The valid range has entered the ultrastrong coupling regime, $g/\omega_{b}\geq 0.1$ [61]. The value distinguished by the black box corresponds to Fig. 2(b). In Fig. 3(c), $g_{\rm eff}$ is valid at the range of $G/\omega_{b}\leq 0.3$ when $r=0$. As $r$ increases to $0.25$, the valid range reduces to $G/\omega_{b}\leq 0.24$. Similarly, the energy shift $\delta$ in Eq. (5) can also be justified by Figs. 3(b) and 3(d). In Fig. 3(b), the analytical $\delta$ shows a slight deviation from the numerical results, but this deviation gradually decreases as $g$ and $r$ increase. In Fig. 3(d), it is found that the energy shift $\delta$ is valid at $G/\omega_{b}\leq 0.3$, whatever the parameter $r$ is $0$ or $0.25$.

Similar results about $g_{\rm eff}$ and $\delta$ are plotted as a function of the squeezing parameter $r$ in Figs. 3(e) and (f), respectively. It is evident that both the analytical results for $g_{\rm eff}$ and $\delta$ match well with the numerical results at $\Delta_{m}=3\omega_{b}$. The effective coupling $g_{\rm eff}$ is significantly amplified as $r$ increases, which can improve the photon-phonon squeezing level discussed later. However, at $\Delta_{m}=0.5\omega_{b}$, there exists a slight distinction between the numerical and analytical results for both $g_{\rm eff}$ or $\delta$. This distinction is mainly attributed to the constraint of the parameter setting, which yields that the detuning between magnon and phonon, $|\Delta_{m}/\cosh(2r)-\omega_{b}|$, is not sufficiently larger than the coupling strengths. Furthermore, $g_{\rm eff}$ decreases as $r$ increases, leading to a low two-mode squeezing level.

Using the effective Hamiltonian in Eq. (3), one can generate naturally and directly the photon-phonon TMSS. In this section, we take the open-quantum-system framework to discuss the dynamics of TMSS generation and explain why the asymptotic two-mode squeezing can be obtained even in an unstable dynamical evolution. Under the standard assumptions, i.e., Markovian approximation and structure-free environment at zero temperature, the dynamics of the quantum system are governed by the quantum Langevin equation (QLE), written in a matrix form

$$\dot{u}^{\rm eff}(t)=A_{\rm eff}u^{\rm eff}(t)+n^{\rm eff}(t),$$

(12)

where $u^{\rm eff}(t)$ is the same as Eq. (10). The transition matrix $A_{\rm eff}(t)$ is

$\displaystyle A_{\rm eff}=\begin{bmatrix}-\kappa_{a}&0&g_{\rm eff}&0\\ 0&-\kappa_{a}&0&-g_{\rm eff}\\ g_{\rm eff}&0&-\kappa_{b}&0\\ 0&-g_{\rm eff}&0&-\kappa_{b}\end{bmatrix}$

(13)

where $\kappa_{a}$ and $\kappa_{b}$ are the decay rates of the modes $a$ and $b$, respectively. $n^{\rm eff}(t)=[\sqrt{2\kappa_{a}}X^{in}_{a}(t),\sqrt{2\kappa_{a}}Y^{in}_{a}(t),\sqrt{2\kappa_{b}}X^{in}_{b}(t),\sqrt{2\kappa_{b}}Y^{in}_{b}(t)]^{T}$ is the vector of Gaussian noise operators, and $X^{in}_{o}=(o_{in}+o^{\dagger}_{in})/\sqrt{2},Y^{in}_{o}=(o_{in}-o^{\dagger}_{in})/i\sqrt{2},o=a,b$. $a_{in}$ and $b_{in}$ are characterized by their covariance functions, $\langle o_{in}(t)o^{\dagger}_{in}(t^{\prime})\rangle=[N_{o}+1]\delta(t-t^{\prime})$ and $\langle o^{\dagger}_{in}(t)o_{in}(t^{\prime})\rangle=N_{o}\delta(t-t^{\prime})$, where $N_{o}$ is the mean population of mode $o$ at the thermal equilibrium state.

The input zero-mean quantum Gaussian noises yield the quantum state as a zero-mean Gaussian state, which can be completely characterized by a $4\times 4$ covariance matrix (CM) $V^{\rm eff}(t)$. By virtue of the QLE in Eq (12), the dynamics of the CM $V^{\rm eff}(t)$ satisfies

$$\dot{V}^{\rm eff}(t)=A_{\rm eff}V^{\rm eff}(t)+V^{\rm eff}(t)A_{\rm eff}^{T}+D^{\rm eff}.$$

(14)

The elements of $V^{\rm eff}(t)$ are defined as

$$V^{\rm eff}_{ij}(t)=\frac{\langle u^{\rm eff}_{i}(t)u^{\rm eff}_{j}(t)+u^{\rm eff}_{j}(t)u^{\rm eff}_{i}(t)\rangle}{2},\\$$

(15)

where $u^{\rm eff}_{i}(t)$ is the $i$-term of $u^{\rm eff}(t)$ and $i=1,2,3,4$. $D^{\rm eff}=Diag[\kappa_{a}(2N_{a}+1),\kappa_{a}(2N_{a}+1),\kappa_{b}(2N_{b}+1),\kappa_{b}(2N_{b}+1)]$ is the diffusion matrix, which is defined through $D^{\rm eff}_{ij}(t)=\langle n^{\rm eff}_{i}(t)n^{\rm eff}_{j}(t)+n^{\rm eff}_{j}(t)n^{\rm eff}_{i}(t)\rangle/2$. In the stable condition, the CM is invariant under time evolution, i.e, $\dot{V}^{\rm eff}=0$ in Eq. (14), which requires $g^{2}_{\rm eff}<\kappa_{a}\kappa_{b}$.

Assume the photon and phonon are both in vacuum states initially, which can be realized by precooling them to their respective ground states [62]. Then the initial CM can be written as $V^{\rm eff}(0)=I_{4}/2$, $I_{4}$ is an identity matrix with four dimensions. Under this initial condition, the non-zero matrix elements in $V^{\rm eff}(t)$ can be solved as

	$\displaystyle V^{\rm eff}_{11}(t)$	$\displaystyle=C_{+}(1-\sin\varphi)e^{(\Omega-\kappa_{a}-\kappa_{b})t}-C_{0}\cos\varphi e^{-(\kappa_{a}+\kappa_{b})t}$		(16)
		$\displaystyle+C_{-}(1+\sin\varphi)e^{-(\Omega+\kappa_{a}+\kappa_{b})t}+c_{a};$
	$\displaystyle V^{\rm eff}_{33}(t)$	$\displaystyle=C_{+}(1+\sin\varphi)e^{(\Omega-\kappa_{a}-\kappa_{b})t}+C_{0}\cos\varphi e^{-(\kappa_{a}+\kappa_{b})t}$
		$\displaystyle+C_{-}(1-\sin\varphi)e^{-(\Omega+\kappa_{a}+\kappa_{b})t}+c_{b};$
	$\displaystyle V^{\rm eff}_{13}(t)$	$\displaystyle=C_{+}\cos\varphi e^{(\Omega-\kappa_{a}-\kappa_{b})t}-C_{0}\sin\varphi e^{-(\kappa_{a}+\kappa_{b})t}$
		$\displaystyle-C_{-}\cos\varphi e^{-(\Omega+\kappa_{a}+\kappa_{b})t}+c,$

and $V^{\rm eff}_{22}(t)=V^{\rm eff}_{11}(t),V^{\rm eff}_{44}(t)=V^{\rm eff}_{33}(t),V^{\rm eff}_{24}(t)=-V^{\rm eff}_{13}(t)$. The parameters are defined as

	$\displaystyle\Omega$	$\displaystyle=\sqrt{4g^{2}_{\rm eff}+(\kappa_{a}-\kappa_{b})^{2}},\quad\tan\varphi=\frac{\kappa_{a}-\kappa_{b}}{2g_{\rm eff}},$		(17)
	$\displaystyle C_{\pm}$	$\displaystyle=\pm\frac{\kappa_{+}\mp\sin\varphi\kappa_{-}}{4[\Omega\mp(\kappa_{a}+\kappa_{b})]}+\frac{1}{4},\quad C_{0}=\frac{\cos\varphi\kappa_{-}}{2(\kappa_{a}+\kappa_{b})},$
	$\displaystyle\kappa_{\pm}$	$\displaystyle=\kappa_{a}(2N_{a}+1)\pm\kappa_{b}(2N_{b}+1).$

And

	$\displaystyle c_{o}$	$\displaystyle=N_{o}+\frac{1}{2}+\frac{g_{\rm eff}}{\kappa_{o}}c,\quad o=a,b$		(18)
	$\displaystyle c$	$\displaystyle=\frac{g_{\rm eff}\kappa_{a}\kappa_{b}(N_{a}+N_{b}+1)}{(\kappa_{a}\kappa_{b}-g^{2}_{\rm eff})(\kappa_{a}+\kappa_{b})},$

which are also the solutions of the matrix elements $V^{\rm eff}_{11}$, $V^{\rm eff}_{33}$, and $V^{\rm eff}_{13}$ in the stable regime, respectively, obtained by setting $\dot{V}^{\rm eff}=0$ in Eq. (14). These stable CM elements are the asymptotic values as $t\to\infty$.

With the CM definition in Eq. (15) and its solution in Eq. (16), the variance of the quadrature operator $X$ in Eq. (6) can be derived as

	$\displaystyle\Delta X(t)$	$\displaystyle=\frac{1}{2}[V^{\rm eff}_{11}(t)+V^{\rm eff}_{33}(t)+2V^{\rm eff}_{13}(t)]$		(19)
		$\displaystyle=(1+\cos\varphi)C_{+}e^{(\Omega-\kappa_{a}-\kappa_{b})t}-\sin\varphi C_{0}e^{-(\kappa_{a}+\kappa_{b})t}$
		$\displaystyle+(1-\cos\varphi)C_{-}e^{-(\Omega+\kappa_{a}+\kappa_{b})t}+C,$

where

$\displaystyle C$

$\displaystyle=\frac{1}{2}(N_{a}+N_{b}+1)\frac{\kappa_{a}\kappa_{b}(2g_{\rm eff}+\kappa_{a}+\kappa_{b})}{(\kappa_{a}\kappa_{b}-g^{2}_{\rm eff})(\kappa_{a}+\kappa_{b})}.$

(20)

In the stable regime, $g^{2}_{\rm eff}<\kappa_{a}\kappa_{b}$, one can demonstrate that the exponent factor $\Omega-\kappa_{a}-\kappa_{b}$ in Eq. (19) is negative through the definition of $\Omega$ in Eq. (17). That leads to $\Delta X(\infty)=C$, consistent with the result obtained by $\dot{V}^{\rm eff}(t)=0$. When the photon decay rate is larger than the phonon decay rate, i.e., $\kappa_{a}>\kappa_{b}$, the minimal value of $C$ can be obtained as

$$C_{\rm min}=\frac{1}{2}(N_{a}+N_{b}+1)\frac{\kappa_{a}}{\kappa_{a}+\kappa_{b}}$$

(21)

at $g_{\rm eff}=-\kappa_{b}$. Even at the zero temperature, $N_{a}=N_{b}=0$, the minimum $C_{\rm min}>0.25$. The value $0.25$ corresponds to the upper bound of the squeezing level, $S=3$dB, under the stable condition. The squeezing level $S$ in the decibel unit is defined by $S=-10\log_{10}(\Delta X/\Delta X_{zp})$ [63], where $\Delta X_{zp}=0.5$ is the standard fluctuation in the zero-point level. The decay rates satisfy $\kappa_{a}\gg\kappa_{b}$ in the recent cavity magnomechanical system [14], resulting in $C_{\rm min}\approx 0.5$ even when $N_{a}=N_{b}=0$. It implies that $X$ cannot be squeezed under stable conditions in this specific experimental platform.

In the unstable regime, $g^{2}_{\rm eff}>\kappa_{a}\kappa_{b}$, all the CM elements $V^{\rm eff}_{11}$, $V^{\rm eff}_{33}$, and $V^{\rm eff}_{13}$ in Eq. (16) exhibit exponential divergence due to the exponential factor $\Omega-\kappa_{a}-\kappa_{b}>0$. These are clearly illustrated by their respective numerical results, shown by a blue-solid line, a red-dashed line, and a black dash-dotted line in Fig. 4(a). The corresponding variance $\Delta X(t)$ in Eq. (19) is depicted by a blue solid line in Fig. 4(b). One can observe that it initially decreases until it reaches its minimum value $\Delta X(\tau)$, where the moment $\tau$ can be analytically determined by setting the derivation $\dot{\Delta X}(\tau)=0$. After reaching its minimum, $\Delta X(t)$ increases exponentially, and $\Delta X(\infty)\to+\infty$.

However, both the CM and the variance $\Delta X(t)$ instabilities do not imply nonstationary TMSS. To find a stationary TMSS with a higher squeezing level, we define a general two-mode squeezing operator $X_{\phi}=\cos\phi X_{a}+\sin\phi X_{b}$, where $\phi$ is an angle to optimize. With the CM elements (16), its variance $\Delta X_{\phi}=\langle X_{\phi}^{2}\rangle-\langle X_{\phi}\rangle^{2}$ can be described as

	$\displaystyle\Delta X_{\phi}(t)$	$\displaystyle=\cos^{2}\phi V^{\rm eff}_{11}(t)+\sin^{2}\phi V^{\rm eff}_{33}(t)+\sin(2\phi)V^{\rm eff}_{13}(t)$		(22)
		$\displaystyle=C_{+}(1-\sin\tilde{\varphi})e^{(\Omega-\kappa_{a}-\kappa_{b})t}-C_{0}\cos\tilde{\varphi}e^{-(\kappa_{a}+\kappa_{b})t}$
		$\displaystyle+C_{-}(1+\sin\tilde{\varphi})e^{-(\Omega+\kappa_{a}+\kappa_{b})t}+C_{\phi},$

where $\tilde{\varphi}=\varphi-2\phi$ and $C_{\phi}=\cos^{2}\phi c_{a}+\sin^{2}\phi c_{b}+\sin(2\phi)c$, $c_{a},c_{b},c$ are constants in Eq. (18).

From Eq. (22), one can find that the exponential divergence term of $\Delta X_{\phi}(t)$ can be canceled at an optimized angle $\tilde{\varphi}=\pi/2$, i.e., the angle $\phi$ satisfies

$$\tan(2\phi)=-\cot(\varphi)=\frac{2g_{\rm eff}}{\kappa_{b}-\kappa_{a}}.$$

(23)

Specifically, $\Delta X_{\phi}(t)=\Delta X(t)$ at $\kappa_{a}=\kappa_{b}$. Under this optimized angle, $\Delta X_{\phi}(t)$ turns into

$$\Delta X_{\phi}(t)=\frac{1}{2}+2C_{-}e^{-(\Omega+\kappa_{a}+\kappa_{b})t}-2C_{-}.$$

(24)

It is evident that $\Delta X_{\phi}(0)=0.5$, corresponding to the standard fluctuation in the zero-point level [63]. The condition $\Delta X_{\phi}(t)<\Delta X_{\phi}(0)$ signifies the occurrence of two-mode squeezing during the evolution [63]. The value of $\Delta X_{\phi}(t)$ quantifies the level of the squeezing, a smaller $\Delta X_{\phi}(t)$ yielding a stronger squeezing. Equation (24) also shows an asymptotic stationary squeezing over a long evolution, i.e.,

$$\Delta X_{\phi}(\infty)=\frac{\Omega\kappa_{+}+(\kappa_{a}-\kappa_{b})\kappa_{-}}{2\Omega(\Omega+\kappa_{a}+\kappa_{b})}.$$

(25)

Given the definition of $\Omega$ in Eq. (17), it follows that $\Delta X_{\phi}(\infty)$ decreases as well as the two-mode squeezing enhances as $g_{\rm eff}$ increases.

It is a natural result of the system’s evolution under the two-mode squeezing Hamiltonian (3) in the open quantum system framework. The two-mode squeezing interaction generates the photon-phonon squeezing, gradually increasing the squeezing level over time. In contrast, Markovian noises reduce the squeezing as time progresses. Both factors constitute a competitive mechanism in the open-quantum-system framework, leading to an asymptotic stationary photon-phonon squeezing with an optimized squeezing operator.

In Fig. 4 (b), we plot $\Delta X_{\phi}(t)$ using the effective Hamiltonian by blue dash-dotted line. After a long time evolution $\omega_{b}t\geq 450\approx 1.5\tau$, it tends to stabilize a certain value of $0.05$, and the corresponding squeezing level $S$ is about $10$dB below vacuum fluctuation, which is larger than the upper bound $3$dB in the stable dynamic condition.

We also simulate the logarithmic negativity (LN) to quantify the photon-phonon squeezing (EPR entanglement) [64, 65, 66], which is defined as

$$E_{N}=\rm{Max}\left\{0,-\frac{1}{2}\rm{ln}\left\{2[\mathcal{P}-(\mathcal{P}^{2}-4\det V^{\rm eff})^{1/2}]\right\}\right\}.$$

(26)

$V^{\rm eff}=[V_{a},V_{ab};V^{T}_{ab},V_{b}]$, with $V_{a}$, $V_{b}$, and $V_{ab}$ being the $2\times 2$ blocks of $V^{\rm eff}$, and $\mathcal{P}\equiv\det V_{a}+\det V_{b}-2\det V_{ab}$. Then, the time-dependent LN can be derived as

$$E_{N}(t)=\rm{Max}\left\{0,-\rm{ln}\left[\eta\left(1-\sqrt{1+\delta^{\prime}}\right)\right]\right\}.$$

(27)

where $\eta=V^{\rm eff}_{11}(t)+V^{\rm eff}_{33}(t)$ and $\delta^{\prime}=\{4[V^{\rm eff}_{13}(t)]^{2}-4V^{\rm eff}_{11}(t)V^{\rm eff}_{33}(t)\}/\eta^{2}$. In the unstable dynamical regime $g^{2}_{eff}>\kappa_{a}\kappa_{b}$, using the solutions in Eq. (16), one can demonstrate that $\delta^{\prime}\to 0$ when $t\to\infty$. With the Taylor expansion $\sqrt{1+\delta^{\prime}}\approx 1+\delta^{\prime}/2$ up to the first order of $\delta^{\prime}$, one can finally obtain the LN at infinity

$$E_{N}(\infty)=\rm{Max}\left\{0,-\ln[2\Delta X_{\phi}(\infty)]\right\}.$$

(28)

When asymptotic stationary two-mode squeezing emerges during unstable dynamics, as indicated by $\Delta X_{\phi}(\infty)<0.5$, asymptotic stationary bipartite entanglement between the photon and phonon also arises. The LN increases with the value of $\Delta X_{\phi}(\infty)$ decreases, indicating that stronger two-mode squeezing corresponds to greater entanglement.

The above results obtained by the effective Hamiltonian in Eq. (3) can be further confirmed by the whole system’s dynamics. Similar as Eq. (14), using the full system Hamiltonian $H$ in Eq. (2), the dynamics of the whole system CM $V(t)$ is determined by

$$\dot{V}(t)=AV(t)+V(t)A^{T}+D.$$

(29)

The elements of $V(t)$ are defined as

$$V_{ij}(t)=\frac{\langle u_{i}(t)u_{j}(t)+u_{j}(t)u_{i}(t)\rangle}{2},i,j=1,2\cdot\cdot\cdot 6,$$

(30)

where $u(t)$ is shown in Eq. (7). The transition matrix $A=i\mathcal{L}+\tilde{A}$, where $\mathcal{L}$ is the superoperator in Eq. (8) and $\tilde{A}=Diag[-\kappa_{a},-\kappa_{a},-\kappa_{b},-\kappa_{b},-e^{2r}\kappa_{m},-e^{2r}\kappa_{m}]$. The magnon decay rate is exponentially enlarged due to the Kerr effect, i.e., $\kappa_{m}\rightarrow e^{2r}\kappa_{m}$ [67]. $D=Diag[\kappa_{a}(2N_{a}+1),\kappa_{a}(2N_{a}+1),\kappa_{b}(2N_{b}+1),\kappa_{b}(2N_{b}+1),e^{2r}\kappa_{m}(2N_{m}+1),e^{2r}\kappa_{m}(2N_{m}+1)]$ is the matrix of noises covariance. Then, the dynamics of $\Delta X(t)$ and $\Delta X_{\phi}(t)$ can be obtained by numerically calculating the CM $V(t)$,

	$\displaystyle\Delta X(t)$	$\displaystyle=\frac{1}{2}[V_{11}(t)+V_{33}(t)+2V_{13}(t)],$		(31)
	$\displaystyle\Delta X_{\phi}(t)$	$\displaystyle=\cos^{2}\phi V_{11}(t)+\sin^{2}\phi V_{33}(t)+\sin(2\phi)V_{13}(t).$

The initial condition is $V(0)=I_{6}/2$, and $I_{6}$ is a six-dimensional identity matrix.

Numerical results are shown in Figs. 4(a) and (b). All of the matrix elements in Fig. 4(a), $V_{11}(t)$ (blue solid line with circles), $V_{33}(t)$ (red dashed line with squares), and $V_{13}(t)$ (black dash-dotted line with diamonds), along with the variances in Fig. 4(b), $\Delta X(t)$ (red dashed line with circles) and $\Delta X_{\phi}(t)$ (purple dotted line with squares) obtained using Eq. (29) do match well with the corresponding results via the effective Hamiltonian (3).

We also numerically analyze the asymptotic stationary LN, $E_{N}(\infty)$, using the whole system’s dynamics to evaluate our protocol. For the whole system, the $E_{N}$ is defined by replacing the $V^{\rm eff}$ in Eq. (26) with $V(1:4,1:4)$ in Eq. (29). Besides, we approximate $E_{N}(\infty)$ by evaluating the value after a sufficient long time, specifically $E_{N}(\infty)=E_{N}(2\tau)$, where $\tau$ is the moment when $\Delta X(t)$ in Eq. (19) reaches its minimum, as shown in Fig. 4 (b).

In Fig. 5(a), we illustrate the $E_{N}(\infty)$ in the parameter space of coupling strengths $g$ and $G$. One can observe that a high entanglement LN can be obtained when the coupling strengths $g$ and $G$ are increased. The LN $E_{N}(\infty)$ is greater than $2.5$ for $g,G\geq 0.2\omega_{b}$. In Fig. 5(b), we present $E_{N}(\infty)$ in the parameter space defined by magnon detuning $\Delta_{m}$ and the squeezing parameter $r$. A low entanglement regime $E_{N}(\infty)\leq 0.5$ around $\Delta_{m}\approx\omega_{b}$ is observed, and the low-region enhances with increasing $r$, which is attributed to the invalidity of the effective Hamiltonian (3) constructed in perturbation condition. In the region where $\Delta_{m}<\omega_{b}$, the effective coupling strength $g_{\rm eff}$ decreases as $r$ increases [see Fig. 3(e)], which leads to a reduction in LN with increasing $r$. In contrast, when $\Delta_{m}>\omega_{b}$, the effective coupling $g_{\rm eff}$ increases as $r$ increases, leading consequently to a enhancement of the LN with increasing $r$. However, under the constraint of the large detuning condition, i.e., $|\Delta_{m}/\cosh(2r)-\omega_{b}|\gg g,G$, an even larger $r$ does not always yield better results. A larger $\Delta_{m}$ allows a broader range of $r$ to achieve a high LN. An optimal parameter space $E_{N}(\infty)\geq 2$ exists at $\Delta_{m}\leq 3.5\omega_{b},r\geq 0.15$. From these findings, we conclude that our protocol significantly enhances the asymptotic stationary LN, $E_{N}$, from the previously reported values of approximately $0.1-0.3$ [18]. It can also exceed the maximum theoretical limit at the stable dynamic condition, i.e., $0.69$ [39].

Our protocol is primarily centered on the cavity magnomechanical system, focusing on constructing the magnon-assisted photon-phonon squeezing. In recent experiments [14, 15, 16], the coupling strength between photon and magnon, $g\sim 1-10$MHz, is about $0.1\omega_{b}$, where the phonon frequency $\omega_{b}\sim 10-100$MHz. The single-excitation magnon-phonon coupling $g_{mb}$ is related to the volume of the YIG sphere and the biased magnetic field directions. For a YIG sphere with a diameter $100\mu$m, the coupling strength is approximately $g_{mb}\sim 0.1$Hz [14]. The Rabi frequency of driving is defined as $\Omega\equiv\sqrt{\kappa P_{d}/(\hbar\omega_{d})}$ [15, 16], where $\kappa$ is the cavity decay rate associated with the driving point, $P_{d}$ and $\omega_{d}$ represent the power and frequency of the microwave drive, respectively. In recent experiments [15], $\kappa/2\pi\sim 0.1$MHz and $P_{d}\sim 20-30$dBm ($100-1000$mW). They give rise a Rabi frequency $\Omega\sim 10^{14}-10^{15}$Hz, corresponding to $|\langle m\rangle|\approx 10^{6}-10^{7}$. The enhanced coupling between magnon and phonon is given by $G\equiv g_{mb}|\langle m\rangle|\sim 0.1\omega_{b}$. Consequently, the effective coupling strength satisfies $|g_{\rm eff}|\sim 0.01\omega_{b}$. These values are consistent with the parameters used in Figs. 2 and 3. The magnon Kerr effect can be enhanced by reducing the volume of the YIG sphere [8, 68, 58]. For specifical parameters $K_{m}\sim 10$nHz, $\Delta_{m}=3\omega_{b}$, and $|\langle m\rangle|=10^{7}$, the squeezing parameter $r$ roughly equals $0.05$. The decay rates of phonon and photon modes are $\kappa_{b}\sim 10^{-5}-10^{-4}\omega_{b}$ and $\kappa_{a}\sim 100\kappa_{b}$, respectively. Thus, the cooperativity $C\equiv g^{2}_{\rm eff}/\kappa_{a}\kappa_{b}\gtrsim 100$. The unstable dynamical regime in our protocol, characterized by $g^{2}_{\rm eff}>\kappa_{a}\kappa_{b}$, can be successfully realized. Besides, the magnon decay rate, $\kappa_{m}\sim 1\rm{MHz}$, is challenging to reduce further due to intrinsic damping [2]. However, its impact is insignificant as the magnon primarily serves as an interface. At a low temperature of $T\sim 10\rm{mK}$, the thermal occupations of photon, magnon, and phonon are respective $N_{a}\approx 0$, $N_{m}\approx 0$, and $N_{b}\approx 10$, consistent with the parameters used in Figs. 4 and 5.

In summary, we have presented a protocol for generating photon-phonon squeezing in the cavity magnomechanics, where the magnon in the YIG sphere is coupled to both microwave photons and the mechanical vibration modes in the same sphere. Our protocol offers significant advantages in terms of controllability within the system, and the Kerr effect of magnon can further enhance the squeezing level in a proper regime. This magnon-assisted protocol relies on the effective two-mode squeezing Hamiltonian for coupling photons and phonons. We apply an interesting method by diagonalizing the whole system’s Liouvillian superoperator to numerically confirm the validity of the effective Hamiltonian, which is beneficial to more physics with nonconservative excitations. In the open-quantum-system framework, we derive the process for generating TMSS with the effective Hamiltonian. Our analysis demonstrates that the asymptotic stationary TMSS with high squeezing levels can be achieved even in unstable dynamics. Our work provides an important implementation of TMSS generation in a solid system under realistic noises. It extends the application of cavity magnomechanics as a promising hybrid platform for quantum information processing.

In addition to the cavity magnomechanical system, our protocol can be extended to other quantum systems. For instance, we can utilize a mechanical interface to realize the microwave-optical photon squeezed state [69, 70] or create the photon-magnon squeezed state [57]. Our scheme presents an extendable framework to create TMSS, which will be widely applied in quantum information processing and quantum metrology using bosonic systems.

We acknowledge financial support from the National Science Foundation of China (Grant No. 12404405) and the Science Foundation of Hebei Normal University of China (Grant No. L2024B10).