# Efficient quantum algorithms for solving quantum linear system problem

Hefeng Wang1<EMAIL_ADDRESS>Hua Xiang2<EMAIL_ADDRESS>1Department
of Applied Physics, School of Physics, Xi’an Jiaotong University and Shaanxi
Province Key Laboratory of Quantum Information and Quantum Optoelectronic
Devices, Xi’an, 710049, China
MOE Key Laboratory for Nonequilibrium Synthesis and Modulation of Condensed
Matter, Xi’an Jiaotong University, Xi’an, 710049, China  2School of
Mathematics and Statistics, Wuhan University, Wuhan, 430072, China

###### Abstract

We transform the problem of solving linear system of equations
$A\mathbf{x}=\mathbf{b}$ to a problem of finding the right singular vector
with singular value zero of an augmented matrix $C$, and present two quantum
algorithms for solving this problem. The first algorithm solves the problem
directly by applying the quantum eigenstate filtering algorithm with query
complexity of $O\left(s\kappa\log\left(1/\epsilon\right)\right)$ for a
$s$-sparse matrix $C$, where $\kappa$ is the condition number of the matrix
$A$, and $\epsilon$ is the desired precision. The second algorithm uses the
quantum resonant transition approach, the query complexity scales as
$O\left[s\kappa+\log\left(1/\epsilon\right)/\log\log\left(1/\epsilon\right)\right]$.
Both algorithms meet the optimal query complexity in $\kappa$, and are simpler
than previous algorithms.

##  I Introduction

Solving linear system of equations (LSE) is one of the fundamental problems in
scientific computation. Given an $N\times N$ matrix $A$ and a vector
$\mathbf{b}$, the task is to find a vector $\mathbf{x}$ such that
$A\mathbf{x}=\mathbf{b}$. Solving high-dimensional LSE is expensive on a
classical computer. Classical linear solvers can be categorized into the
direct methods and the iterative methods Golub0 ; demmel . The direct methods
such as Gaussian elimination solve LSE with runtime scales as $O(N^{3})$.
There exists more efficient classical linear system solver that scales as
$O(N^{\nu})$ where $\nu\leqslant 2.373$ CW ; Gall , but it is difficult to
utilize in practice due to numerical instability. The iterative methods show
great advantages when they converge quickly, the iteration number is an
indicator for the efficiency of these methods. E.g., for a symmetric positive-
definite problem, the steepest descent method needs
$O(\kappa\log(1/\epsilon))$ iterations and the conjugate gradient method needs
$O(\sqrt{\kappa}\log(1/\epsilon))$, where $\kappa$ is the condition number of
the matrix $A$ defined as the ratio between the largest and the smallest
singular values of $A$, or $\left\|A\right\|\left\|A^{-1}\right\|$, where
$\left\|\cdot\right\|$ denotes vector or matrix $2$-norm, and $\epsilon$ is
the desired precision of the solution.

Quantum computation provides an efficient way of solving quantum linear system
problem (QLSP), which aims to prepare a quantum state that is proportional to
the solution vector of a given LSE. Quantum algorithms for solving the QLSP
either apply the matrix inversion operator $A^{-1}$ directly on the state
$|b\rangle$ to obtain $|x\rangle$, e.g. the HHL algorithm HHL ; childs , or
transform the LSE to an eigenvalue problem where the ground state of the
problem Hamiltonian is the state $|x\rangle$ AQCLSE ; Lin . Recently, a
quantum eigenstate filtering (QEF) algorithm Lin  was proposed for solving the
QLSP by combining with adiabatic quantum computing (AQC) or quantum Zeno
effect (QZE), and achieves near optimal query complexity of
$O\left(s\kappa\text{poly}\log\left(\kappa\right)\log\left(1/\epsilon\right)\right)$,
where $s$ is the sparsity of the matrix $A$. The QEF algorithm approximates a
spectral projection operator by using the quantum signal processing (QSP) QSP
method and project out the quantum state $|x\rangle$ from an initial state
prepared through an AQC procedure. Such a procedure can also be realized
through QZE by applying a sequence of QEF procedures following the same path
as that of the AQC method.

In this work, we solve the QLSP in a simpler way. We first prove that the
solution vector to the LSE is contained in a vector proportional to the right
singular vector of the augmented matrix
$C=\left(\begin{array}[]{cc}A&\beta^{-1}\mathbf{b}\end{array}\right)$ with
corresponding singular value $0$, where $\beta$ is a parameter. Then we
construct a Hermitian matrix

$B=\left(\begin{array}[]{cc}0&C\\\ C^{{\dagger}}&0\end{array}\right),$ (1)

and transform the QLSP to a problem of finding the eigenstate of the matrix
$B$ with eigenvalue $0$. In our previous work, we proposed a quantum algorithm
for finding an eigenstate with known eigenvalue of a Hermitian matrix through
quantum resonant transition (QRT) whf0 ; whf1 . Here, we set a simple initial
state whose overlap with the desired eigenstate of the matrix $B$ is
$O\left(1\right)$, then both the QEF and the QRT algorithms can be applied for
solving the QLSP. The QEF algorithm has query complexity of
$O\left(s\kappa\log\left(1/\epsilon\right)\right)$, and the query complexity
of the QRT algorithm scales as
$O\left[s\kappa+\log\left(1/\epsilon\right)/\log\log\left(1/\epsilon\right)\right]$
for solving the problem. Both algorithms achieve the optimal query complexity
in $\kappa$. In previous quantum algorithms, the complexity for preparing the
state $|b\rangle$ and querying the matrix $A$ are considered separately, while
they are combined together in our work. The algorithms do not need to use
complex procedures such as phase estimation, amplitude amplification, AQC or
QZE, thus they are easier to implement experimentally than previous
algorithms.

##  II LSE and singular value decomposition

The LSE $A\mathbf{x}=\mathbf{b}$ can be written in form of

$C\left(\begin{array}[]{c}\mathbf{x}\\\ -\beta\end{array}\right)=0.$ (2)

By performing singular value decomposition (SVD), we have $C=SDV^{{\dagger}}$,
where $S$ and $V$ are unitary matrices of dimension $N$ and $(N+1)$,
respectively, and
$D=\left(\begin{array}[]{cc}D^{\prime}&\mathbf{0}\end{array}\right)$,
$D^{\prime}=$diag$(\sigma_{1},\ldots,\sigma_{N})$, and $\mathbf{0}$ is the
zero column vector of dimension $N$,
$\sigma_{1}\geqslant\cdots\geqslant\sigma_{N}$ are the singular values of the
matrix $C$. Correspondingly, $\mathbf{s}_{1},\ldots,\mathbf{s}_{N}$ are column
vectors of $S$, and $\mathbf{v}_{1},\ldots,\mathbf{v}_{N+1}$ are column
vectors of $V$.

Theorem 1.–Suppose $\mathbf{x}$ is the solution to the equation
$A\mathbf{x}=\mathbf{b}$, where $A$ is an $N\times N$ nonsingular matrix,
$\mathbf{v}_{N+1}$ is the right singular vector corresponding to
$C\mathbf{v}_{N+1}=\mathbf{0}$, where
$C=\left(\begin{array}[]{cc}A&\beta^{-1}\mathbf{b}\end{array}\right)$, then
the vector
$\left(\begin{array}[]{cc}\mathbf{x^{\text{{T}}},}&-\beta\end{array}\right)^{\text{{T}}}$
that satisfies Eq. ($1$) is proportional to the singular vector
$\mathbf{v}_{N+1}$ of the matrix $C$.

Proof.–Since the rank of the matrix $A$ is $N$, let $\sigma_{N+1}=0$, the
singular values of the matrix $C$ satisfy
$\sigma_{1}\geqslant\cdots\geqslant\sigma_{N}>\sigma_{N+1}=0$. The matrix $C$
can be written as
$C\mathbf{=}\sum_{i=1}^{N}\sigma_{i}\mathbf{s}_{i}\mathbf{v}_{i}^{{\dagger}}$.
The column vectors $\mathbf{v}_{1},\ldots,\mathbf{v}_{N+1}$ of the unitary
matrix $V$ are orthogonal to each other, therefore to satisfy Eq. ($1$), the
vector
$\left(\begin{array}[]{cc}\mathbf{x^{\text{{T}}},}&-\beta\end{array}\right)^{\text{{T}}}$
must be orthogonal to the vectors $\mathbf{v}_{1},\ldots,\mathbf{v}_{N}$, thus
it is proportional to the vector $\mathbf{v}_{N+1}$. $\square$

Let $\left\|A\right\|\leq 1$, we perform SVD for the matrix $A$,
$\bar{S}^{{\dagger}}A\bar{V}=$diag$\left(\bar{\sigma}_{1},\bar{\sigma}_{2},\ldots,\bar{\sigma}_{N}\right)$.
The singular values of the matrices $A$ and $C$ satisfy the following relation
(Golub0, , Corollary 8.6.3, p487):

$\sigma_{1}\geqslant\bar{\sigma}_{1}\geqslant\sigma_{2}\geqslant\bar{\sigma}_{2}\ldots\geqslant\bar{\sigma}_{N}>\sigma_{N+1}=0.$
(3)

Then the energy gap between the singular states $\mathbf{v}_{N+1}$ and
$\mathbf{v}_{N}$ is

$\Delta=\sigma_{N}-\sigma_{N+1}\geqslant\Delta^{\ast}=\bar{\sigma}_{N}=\frac{1}{\kappa}.$
(4)

Suppose $\left\|\mathbf{b}\right\|=1$, then
$\left\|\mathbf{x}\right\|=\left\|A^{-1}\right\|\left\|\mathbf{b}\right\|$ is
in the range of $\left[1,\kappa\right]$. The vector
$\left(\begin{array}[]{cc}\mathbf{x^{\text{{T}}},}&-\beta\end{array}\right)^{\text{{T}}}$
can be normalized as
$\mathbf{v}_{N+1}=d_{0}\left(\begin{array}[]{cc}\mathbf{x^{\text{{T}}}/\left\|\mathbf{x}\right\|,}&0\end{array}\right)^{\text{{T}}}+d_{1}\left(\begin{array}[]{cc}\mathbf{0^{\text{{T}}},}&1\end{array}\right)^{\text{{T}}}$,
where $d_{0}/d_{1}=\mathbf{\left\|\mathbf{x}\right\|/}\beta$.

In the algorithm, the solution state $\mathbf{x/\left\|\mathbf{x}\right\|}$ is
obtained in two steps, we first obtain the state $\mathbf{v}_{N+1}$, then
extract the solution state $\mathbf{x/\left\|\mathbf{x}\right\|}$ from
$\mathbf{v}_{N+1}$. The choice of the parameter $\beta$ has dual influence on
both the overlap of $\mathbf{v}_{N+1}$ with the initial state and the post-
selection of solution state. We first set the parameter $\beta$ in the order
of $O\left(\kappa\right)$, thus the component $d_{1}$ is guaranteed to be in
the order of $O\left(1\right)$, and run the algorithm to obtain the state
$\mathbf{v}_{N+1}$. We run the algorithm and perform measurement on the state
$\mathbf{v}_{N+1}$ for a number of times to estimate the value of $d_{1}$.
This can be achieved by measuring the $(N+1)$-th component of the vector
$\mathbf{v}_{N+1}$ to obtain the probability of its outcome being $1$, and
obtain $d_{0}=\sqrt{1-d_{1}^{2}}$. After that we reset the parameter $\beta$
such that $\beta$ is in the order of $\mathbf{\left\|\mathbf{x}\right\|}$ and
run the algorithm again. Then both the components $d_{0}$ and $d_{1}$ are in
the order of $O\left(1\right)$. The state
$\mathbf{x^{\text{{T}}}/\left\|\mathbf{x}\right\|}$ can be extracted
efficiently from the state $\mathbf{v}_{N+1}$.

In the following we describe the algorithms for obtaining the state
$\mathbf{v}_{N+1}$. The matrix $B$ has eigenvalues and eigenstates as follows:
$B|\varphi_{j}\rangle=E_{j}|\varphi_{j}\rangle$, $j=1,\ldots,2N+1$, and
$E_{1}=-\sigma_{1}$, $\ldots$, $E_{N}=-\sigma_{N}$, $E_{N+1}=0$,
$E_{N+2}=\sigma_{N}$, $\ldots$, $E_{2N+1}=\sigma_{1}$, and
$|\varphi_{1}\rangle=(\mathbf{s}_{1},-\mathbf{v}_{1})$, $\ldots$,
$|\varphi_{N}\rangle=(\mathbf{s}_{N},-\mathbf{v}_{N})$,
$|\varphi_{N+1}\rangle=|\mathbf{v}_{N+1}\rangle$,
$|\varphi_{N+2}\rangle=(\mathbf{s}_{N},\mathbf{v}_{N})$, $\ldots$,
$|\varphi_{2N+1}\rangle=(\mathbf{s}_{1},\mathbf{v}_{1})$. The problem of
solving the LSE is transformed to finding the eigenstate
$|\mathbf{v}_{N+1}\rangle=(\mathbf{0},\mathbf{v}_{N+1})$ of the matrix $B$
with eigenvalue $0$. The gap between the eigenstate $|\mathbf{v}_{N+1}\rangle$
and its nearest neighboring eigenstate $|\varphi_{N+2}\rangle$ is in the order
of $O(\kappa^{-1})$ in the worst case. We set the initial state of the quantum
circuit as $|\mathbf{1}\rangle=(\mathbf{0},\mathbf{0},1)$, whose overlap with
the desired eigenstate $|\mathbf{v}_{N+1}\rangle$ of the matrix $B$ is
$d_{1}=\langle\mathbf{1}|\mathbf{v}_{N+1}\rangle$ which scales as $O(1)$ in
the algorithm. We then apply the QEF and the QRT algorithms for obtaining the
state $|\mathbf{v}_{N+1}\rangle$ from the initial state $|\mathbf{1}\rangle$,
and extract the solution state $\mathbf{x/\left\|\mathbf{x}\right\|}$ from the
state $\mathbf{v}_{N+1}$.

##  III Solving QLSP through QEF

The QEF algorithm projects out the desired eigenstate of a Hermitian matrix
from an initial state by implementing an eigenstate-filtering function using
the QSP method, which is a powerful algorithm for implementing a polynomial
function of matrices on a quantum computer with minimal number of ancilla
qubits. The matrix $B$ is encoded in a unitary matrix by using the block-
encoding technique, then the QSP method is applied to implement the QEF
function to project out the eigenstate $|\mathbf{v}_{N+1}\rangle$ from the
initial state.

The matrix $B$ is represented on an $n$-qubit quantum register, an
($m+n$)-qubit unitary operator $U_{B}$ is called a
$(\alpha,m,\epsilon)$-block-encoding of the matrix $B$ Qubitize , if

$\left\|B-\alpha\left(\langle 0^{m}|\otimes
I_{n}\right)U_{B}\left(|0^{m}\rangle\otimes I_{n}\right)\right\|\leq\epsilon,$
(5)

where $I_{n}$ is an $n$-qubit identity matrix. Block-encoding of $B$ can also
be written in form of

$U_{B}=\left(\begin{array}[]{cc}B/\alpha&\cdot\\\
\cdot&\cdot\end{array}\right),$ (6)

where $\alpha$ is a parameter such that $\left\|B/\alpha\right\|\leq 1$. It is
determined by the largest eigenvalue of the matrix $B$, or the largest
singular value $\sigma_{1}$ of the matrix $C$. We estimate the upper bound of
$\sigma_{1}$ as follows:

$\displaystyle\sigma_{\max}^{2}(C)$ $\displaystyle=$
$\displaystyle\max_{\mathbf{y}\neq
0}\frac{\left\|C\mathbf{y}\right\|^{2}}{\left\|\mathbf{y}\right\|^{2}}$ (14)
$\displaystyle=$
$\displaystyle\max_{\mathbf{y}=\left(\begin{array}[]{c}y_{1}\\\
y_{2}\end{array}\right)\neq
0}\frac{\left\|\left(\begin{array}[]{cc}A&\beta^{-1}\mathbf{b}\end{array}\right)\left(\begin{array}[]{c}y_{1}\\\
y_{2}\end{array}\right)\right\|^{2}}{\left\|\left(\begin{array}[]{c}y_{1}\\\
y_{2}\end{array}\right)\right\|^{2}}$ $\displaystyle=$
$\displaystyle\max_{\left(\begin{array}[]{c}y_{1}\\\
y_{2}\end{array}\right)\neq
0}\frac{\left\|Ay_{1}+\beta^{-1}\mathbf{b}y_{2}\right\|^{2}}{\left\|y_{1}\right\|^{2}+\left\|y_{2}\right\|^{2}}$
(17) $\displaystyle\leq$ $\displaystyle\max_{\left(\begin{array}[]{c}y_{1}\\\
y_{2}\end{array}\right)\neq
0}\frac{\left(\left\|Ay_{1}\right\|+\beta^{-1}\left\|\mathbf{b}\right\|\left\|y_{2}\right\|\right)^{2}}{\left\|y_{1}\right\|^{2}+\left\|y_{2}\right\|^{2}}.$
(20)

Then we have:

$\displaystyle\sigma_{\max}\left(C\right)$ $\displaystyle\leq$
$\displaystyle\max_{\left(\begin{array}[]{c}y_{1}\\\
y_{2}\end{array}\right)\neq
0}\frac{\left\|Ay_{1}\right\|+\beta^{-1}\left\|\mathbf{b}\right\|\left\|y_{2}\right\|}{\sqrt{\left\|y_{1}\right\|^{2}+\left\|y_{2}\right\|^{2}}}$
(23) $\displaystyle\leq$ $\displaystyle\max_{\left(\begin{array}[]{c}y_{1}\\\
y_{2}\end{array}\right)\neq
0}\frac{\left\|Ay_{1}\right\|}{\sqrt{\left\|y_{1}\right\|^{2}+\left\|y_{2}\right\|^{2}}}+\max_{\left(\begin{array}[]{c}y_{1}\\\
y_{2}\end{array}\right)\neq
0}\frac{\beta^{-1}\left\|\mathbf{b}\right\|\left\|y_{2}\right\|}{\sqrt{\left\|y_{1}\right\|^{2}+\left\|y_{2}\right\|^{2}}}$
$\displaystyle\leq$ $\displaystyle\max_{y_{1}\neq
0}\frac{\left\|Ay_{1}\right\|}{\left\|y_{1}\right\|}+\max_{y_{2}\neq
0}\frac{\beta^{-1}\left\|\mathbf{b}\right\|\left\|y_{2}\right\|}{\left\|y_{2}\right\|}$
$\displaystyle=$ $\displaystyle\sigma_{\max}\left(A\right)+\beta^{-1}$
$\displaystyle=$ $\displaystyle\bar{\sigma}_{1}+\beta^{-1}.$ (29)

The parameter $\beta$ is in the range of $\left[1,\kappa\right]$, thus the
maximum singular value of the matrix $C$ is no more than $2$.

The polynomial eigenvalue transformation of $B/\alpha$ can be implemented
based on the following theorem Lin ; QSVT :

Theorem 2 (Polynomial eigenvalue transformation with definite parity via
quantum signal processing): Let $U_{B}$ be an ($\alpha$, $m$, $0$)-block-
encoding of the Hermitian matrix $B/\alpha$ and $P\in\mathbb{R}[w]$ be a
degree-$l$ even or odd real polynomial and $\left|P\left(w\right)\right|\leq
1$ for any $w\in[-1,1]$. Then there exists a ($1,m+1,0$)-block-encoding
$U_{\vec{\varphi}}$ of $P(B/\alpha)$ using $l$ queries of $U_{B}$,
$U{}_{B}^{{\dagger}}$, and $O((m+1)l)$ other primitive quantum gates.

Let
$B/\alpha=\sum_{\lambda}\lambda\left|\lambda\right\rangle\left\langle\lambda\right|$,
and define $\Pi:=\left|0\right\rangle\left\langle 0\right|$ acting on an
auxiliary qubit, and $\Pi_{\phi}:=e^{i\phi\left(2\Pi-I\right)}$. For $l$ is
even, $U_{\vec{\varphi}}$ is in form of grandunif

$U_{\vec{\varphi}}=\prod_{k=1}^{l/2}\left(\Pi_{\varphi_{2k-1}}U_{B}^{\dagger}\Pi_{\varphi_{2k}}U_{B}\right)=\left(\begin{array}[]{cc}P\left(B/\alpha\right)&\cdot\\\
\cdot&\cdot\end{array}\right),$ (30)

where
$P\left(B/\alpha\right)=\sum_{\lambda}P\left(\lambda\right)\left|\lambda\right\rangle\left\langle\lambda\right|$
is a polynomial transform of the eigenvalues of $B/\alpha$. The phase factors
$(\varphi_{1},\ldots,\varphi_{l})$ can be calculated efficiently phase1 ;
phase2 ; phase3 ; phase4 ; phase5 . The eigenvalue transform can be used to
project out the desired eigenstate with known eigenvalue and filter out other
unrelated states. For a Hermitian matrix with eigenvalue $\lambda$ that is
known to be separated from other eigenvalues by a gap $\Delta>0$, it has been
shown that the following degree-($l=2k$) polynomial

$R_{k}\left(w;\Delta\right)=\frac{T_{k}\left(-1+2\frac{w^{2}-\Delta^{2}}{1-\Delta^{2}}\right)}{T_{k}\left(-1+2\frac{-\Delta^{2}}{1-\Delta^{2}}\right)}$
(31)

is an optimal polynomial for filtering out the unwanted eigenstates Lin ,
where $T_{k}\left(w\right)$ is the $k$th Chebyshev polynomial of the first
kind. By using this polynomial in eigenvalue transform, the system can be
projected onto the eigenspace corresponding to the eigenvalue $\lambda$.

The eigenstate $|\mathbf{v}_{N+1}\rangle$ of the matrix $B$ is separated from
the nearest eigenstate of $B$ by a minimum gap of $\Delta^{\ast}=\kappa^{-1}$.
We set the initial state of the quantum circuit as
$\left|\mathbf{1}\right\rangle$, and apply the QEF algorithm to project out
the eigenstate $|\mathbf{v}_{N+1}\rangle$, while filtering out all other
eigenstates. The overlap between the states $\left|\mathbf{1}\right\rangle$
and $|\mathbf{v}_{N+1}\rangle$ is $d_{1}=O\left(1\right)$. Therefore the
eigenstate $|\mathbf{v}_{N+1}\rangle$ can be obtained with success probability
$O\left(1\right)$. The number of qubits required is $\left(n+m+1\right)$.

The implementation of the algorithm requires querying the matrix $B$ that
contains the matrix $C$ which is composed of the matrix $A$ and the column
vector $|b\rangle$. The matrix $A$ can be accessed by oracles $O_{A,1}$ and
$O_{A,2}$ as
$O_{A,1}\left|j,l\right\rangle=\left|j,\nu\left(j,l\right)\right\rangle$,
$O_{A,2}\left|j,k,z\right\rangle=\left|j,k,A_{jk}\oplus z\right\rangle$,
($j,k,l,z\in[N]$), where the oracle $O_{A,1}$ accepts a row $j$ index and
calculates the column index $\nu\left(j,l\right)$ of the $l$th nonzero element
in the $j$th row of the matrix $A$, and the oracle $O_{A,2}$ accepts a row $j$
and a column $k$ index and returns the matrix element $A_{jk}$ bc . The vector
$b$ is prepared by using an oracle $O_{b}$ as
$O_{b}\left|0\right\rangle=\left|b\right\rangle$, and its elements can be
accessed by an oracle $O_{b,1}$ as
$O_{b,1}\left|k,z\right\rangle=\left|k,b_{k}\oplus z\right\rangle$,
($k,z\in[N]$), where the oracle $O_{b,1}$ accepts the input $k$ and returns
the $k$th element $b_{k}$ of the vector $\left|b\right\rangle$. The matrix $C$
contains the matrix $A$ and the vector $\left|b\right\rangle$ as the
$\left(N+1\right)$-th column vector. If the matrix $C$ is $s$-sparse, then it
can be accessed by using oracles $O_{C,1}$ and $O_{C,2}$ similar to the
oracles for accessing the matrix $A$ above as
$O_{C,1}\left|j,l\right\rangle=\left|j,f\left(j,l\right)\right\rangle$,
$O_{C,2}\left|j,k,z\right\rangle=\left|j,k,C_{jk}\oplus z\right\rangle$,
($k,l,z\in[N],j\in[N+1]$), where the oracle $O_{C,1}$ calculates
$f\left(j,l\right)$ which is the column index of the $l$th nonzero element in
the $j$th row of the matrix $C$, and the oracle $O_{C,2}$ accepts the input
$\left(j,k\right)$ and returns the matrix element $C_{jk}$ of the matrix $C$.
In this case, a ($s$, $n+2$, $0$)-block-encoding of $C$ can be constructed by
using the $O_{C,1}$ and $O_{C,2}$ childs ; berry . The complexity of the QEF
algorithm scales as
$O\left(\left(\alpha/\Delta^{\ast}\right)\log\left(1/\epsilon\right)\right)$
Lin , therefore the complexity of using the QEF algorithm for obtaining the
state $|\mathbf{v}_{N+1}\rangle$ scales as
$O\left(s\kappa\log\left(1/\epsilon\right)\right)$ by querying the oracles
$O_{C,1}$ and $O_{C,2}$.

##  IV Solving QLSP via QRT

In Ref. whf0 ; whf1 , we proposed a quantum algorithm for finding an
eigenstate with known corresponding eigenvalue of a Hamiltonian based on
quantum resonant transition. By coupling a probe qubit to a system, a resonant
transition occurs when the transition frequency of the probe qubit matches a
transition in the system, and the system is evolved to the eigenstate with
known eigenvalue. This algorithm can be applied for solving the QLSP by
obtaining the eigenstate $|\mathbf{v}_{N+1}\rangle$ of the matrix $B$. The
algorithm requires $n+1$ qubits with one probe qubit and an $n$-qubit register
$R$ representing the matrix $B$. The Hamiltonian of the algorithm is
constructed as

$H=-\frac{1}{2}\omega\sigma_{z}\otimes I_{n}+H_{R}+c\sigma_{x}\otimes I_{n},$
(32)

where

$H_{R}=\varepsilon_{0}|1\rangle\langle
1|\otimes|\mathbf{1}\rangle\langle\mathbf{1}|+|0\rangle\langle 0|\otimes
B\mathbf{,}$ (33)

and $\sigma_{x}$ and $\sigma_{z}$ are the Pauli matrices. The first term in
Eq. ($11$) is the Hamiltonian of the probe qubit, the second term contains the
Hamiltonian of the register $R$ and describes the interaction between the
probe qubit and $R$, and the third term is a perturbation. The parameter
$\varepsilon_{0}$ is used as a reference point to the eigenstate
$|\mathbf{v}_{N+1}\rangle$ of the matrix $B$ with eigenvalue $0$, and $c\ll
1$. The condition for resonant transition between states
$|1\rangle|\mathbf{1}\rangle$ and $|0\rangle|\mathbf{v}_{N+1}\rangle$ is
satisfied as $E_{N+1}-\varepsilon_{0}=\omega$, that is,
$\omega=-\varepsilon_{0}$.

The algorithm is run as follows:

$i$) Set the initial state of the $n+1$ qubits as
$|1\rangle|\mathbf{1}\rangle$, which is an eigenstate of $H_{R}$ with
eigenvalue $\varepsilon_{0}$.

$ii$) Implement the unitary operator $U(t)=\exp\left(-iHt\right)$ by setting
$\omega=1$ and $\varepsilon_{0}=-1$.

$iii$) Read out the state of the probe qubit.

As the resonant transition occurs, the system is approximately in state
$\sqrt{1-p}|1\rangle|\mathbf{1}\rangle+\sqrt{p}|0\rangle|\mathbf{v}_{N+1}\rangle$,
where $p=\sin^{2}\left(ctd_{1}\right)$ is the decay probability of the probe
qubit, and $c<\Delta^{\ast}$ and
$d_{1}=\langle\mathbf{1}|\mathbf{v}_{N+1}\rangle$. By performing a measurement
on the probe qubit, if the probe decays to its ground state $|0\rangle$, it
indicates that a resonant transition occurs and the system evolves to the
state $|0\rangle|\mathbf{v}_{N+1}\rangle$; otherwise if the probe qubit stays
in state $|1\rangle$, it means that the register $R$ remains in state
$|\mathbf{1}\rangle$, then we repeat steps $ii$)-$iii$) until the probe qubit
decays to its ground state $|0\rangle$. The number of times the procedures
need to be repeated is proportional to $1/p$.

Errors are introduced by excitations from the initial state
$|1\rangle|\mathbf{1}\rangle$ to the states $|0\rangle|\varphi_{j}\rangle$
($j=N+2,\cdots,2N+1$). Similar to the method in Ref. wyx , we estimate errors
introduced when running the algorithm with the evolution time
$t=\pi/(2cd_{1})$ by setting $\varepsilon_{0}=-1$ and $\omega=1$ as follows.
Let

$H_{0}=-\frac{1}{2}\sigma_{z}\otimes I_{n}-|1\rangle\langle
1|\otimes|\mathbf{1}\rangle\langle\mathbf{1}|+|0\rangle\langle 0|\otimes B,$
(34)

then the algorithm Hamiltonian $H$ can be written as
$H=H_{0}+c\sigma_{x}\otimes I_{n}$. The Hamiltonian $H_{0}$ is the unperturbed
term and has eigenstates

$H_{0}|1\rangle|\mathbf{1}\rangle=\\!-\frac{1}{2}|1\rangle|\mathbf{1}\rangle,$
(35)

and

$H_{0}|0\rangle|\varphi_{j}\rangle=\left(-\frac{1}{2}\\!+\\!E_{j}\right)|0\rangle|\varphi_{j}\rangle,$
(36)

where $j=1$, $\ldots$, $2N+1$. The system is initialized in state
$|1\rangle|\mathbf{1}\rangle$.

For a probe qubit coupled to a two-level system described by the Hamiltonian
in Eq. ($11$), the maximum transition probability from the ground state to the
excited state of the two-level system becomes higher as the transition
frequency between the two-level system gets closer to the frequency of the
probe qubit. Therefore the upper bound of the transition probability from the
initial state to the states other than $|\varphi_{N+1}\rangle$ can be
estimated by assuming all the other states are degenerate at the eigenstate
$|\varphi_{N+2}\rangle$ and without considering competition of the transition
from the initial state to the target state $|0\rangle|\varphi_{N+1}\rangle$.
Then the upper bound of the error of the algorithm is the transition
probability from the initial state to the state
$|0\rangle|\varphi_{N+2}\rangle$ described by the following Hamiltonian in
basis of
$\left\\{|1\rangle|\mathbf{1}\rangle,|0\rangle|\varphi_{N+2}\rangle\right\\}$:

${H_{err}=}\left(\begin{array}[]{cc}-\frac{1}{2}&c\sqrt{1-d_{1}^{2}}\\\
c\sqrt{1-{d}_{1}^{2}}&-\frac{1}{2}+E_{N+2}\end{array}\right).$ (37)

When the coefficient $c$ is much less than the energy difference between the
ground state and the first excited state of $H_{0}$, the excitation to the
excited states can be described by the Rabi’s formula (cohen, , p414). The
transition probability from the initial state $|1\rangle|\mathbf{1}\rangle$ to
the state $|0\rangle|\varphi_{N+2}\rangle$ is

$\frac{4c^{2}\left(1-d_{1}^{2}\right)}{4c^{2}\left(1-d_{1}^{2}\right)+E_{N+2}^{2}}\sin^{2}\left[\frac{t}{2}\sqrt{4c^{2}\left(1-d_{1}^{2}\right)+E_{N+2}^{2}}\right].$
(38)

The upper bound of the transition probability from the initial state to the
states $|0\rangle|\varphi_{j}\rangle$ ($j=N+2,\cdots,2N+1$) can be estimated
as

$\sum_{j=N+2}^{2N+1}p_{j}\leq\frac{4c^{2}\left(1-d_{1}^{2}\right)}{4c^{2}\left(1-d_{1}^{2}\right)+E_{N+2}^{2}}<\frac{4c^{2}\left(1-d_{1}^{2}\right)}{E_{N+2}^{2}}=\frac{4c^{2}\left(1-d_{1}^{2}\right)}{\sigma_{N}^{2}},$
(39)

which can be controlled as a small number by setting the coefficient $c$ to be
small, since $c<\sigma_{N}$.

The evolution time $t$ of the QRT algorithm is in the order of $1/cd_{1}$ such
that the success probability of the algorithm is $O\left(1\right)$. It scales
as $O\left(\kappa\right)$ since $d_{1}$ is in $O\left(1\right)$, and
$\Delta^{\ast}$ thus $c$ is in the order of $O\left(\kappa^{-1}\right)$ in the
worst case. The complexity of the QRT algorithm is determined by Hamiltonian
simulation of the algorithm $U(t)=\exp\left(-iHt\right)$. By applying a
$(\alpha,m,0)$-block-encoding of $H$ in a unitary matrix, the Hamiltonian
simulation of $H$ is $U(t)=e^{-i\left(H/\alpha\right)\alpha t}$, and $\alpha$
scales as $O\left(1\right)$. The optimal approach for Hamiltonian simulation
is by applying the QSP algorithm, for which the query complexity scales as
$\Theta\left(\alpha
t+\frac{\log\left(1/\epsilon\right)}{\log\left(e+\frac{\log\left(1/\epsilon\right)}{\alpha
t}\right)}\right)$ by accessing the Hamiltonian $H$ grandunif . The number of
times the Hamiltonian $H$ is queried scales as
$O\left[\kappa+\log\left(1/\epsilon\right)/\log\log\left(1/\epsilon\right)\right]$
since $t$ scales as $O\left(\kappa\right)$ and $\alpha$ scales as
$O\left(1\right)$. The query complexity of the algorithm by querying the
oracles that access the Hamiltonian matrix $H$ scales as
$O\left[s\kappa+\log\left(1/\epsilon\right)/\log\log\left(1/\epsilon\right)\right]$
QSP  for an $s$-sparse matrix $C$. The total number of qubits required in this
algorithm is $\left(n+m+2\right)$.

##  V Discussion

In this work, the problem of solving the linear system of equations is
transformed to finding the eigenstate of a Hermitian matrix with eigenvalue
$0$. The QEF and QRT algorithms are used for obtaining this state by setting
an initial state whose overlap with the desired eigenstate is
$O\left(1\right)$. Because of this property, our algorithm have a better
complexity scaling compare with other algorithms where an initial state needs
to be prepared through an AQC procedure. The QEF algorithm is applied directly
on the initial state to project out the solution state without using the AQC
or QZE, thus achieves a better query complexity of
$O\left[s\kappa\log\left(1/\epsilon\right)\right]$ than that of the algorithm
in Lin . The QRT algorithm evolves the initial state to the solution state of
the QLSP with query complexity of
$O\left[s\kappa+\log\left(1/\epsilon\right)/\log\log\left(1/\epsilon\right)\right]$.
The complexity of both algorithms have linear scaling in $\kappa$, which is
optimal since the dependence on $\kappa$ cannot be made sublinear HHL . The
complex procedures such as phase estimation, amplitude amplification, AQC or
QZE used in previous algorithms are not needed, only time-independent
Hamiltonian simulation is used, for which the QSP algorithm provides an
optimal solution, thus our algorithms are easier to implement experimentally.
In the algorithms, we obtain the state $|\mathbf{v}_{N+1}\rangle$ that
contains the solution vector of LSE as
$\mathbf{v}_{N+1}=d_{0}\left(\begin{array}[]{cc}\mathbf{x^{\text{{T}}}/\left\|\mathbf{x}\right\|,}&0\end{array}\right)^{\text{{T}}}+d_{1}\left(\begin{array}[]{cc}\mathbf{0^{\text{{T}}},}&1\end{array}\right)^{\text{{T}}}$.
The contribution of the state
$\left(\begin{array}[]{cc}\mathbf{0^{\text{{T}}},}&1\end{array}\right)^{\text{{T}}}$
in applications, such as calculating expectation value of some operators, can
be removed through some adjustment.

The linear solver is a basic engine in engineering and scientific computing,
and has wide applications in many areas. It paves a way for quantum machine
learning, and acts as an important ingredient in linear regression, Bayesian
inference, least-squares fitting, least squares support vector machine. The
numerical solvers for partial differential equations and ordinary differential
equations are also built on it. After numerical discretization, such as the
finite element method, finite difference method, or finite volume method, one
usually needs to solve a sparse linear system. The quantum algorithms we
presented here can be used as a subroutine in solving these problems.

###### Acknowledgements.

This work was supported by the National Key Research and Development Program
of China (No. 2021YFA1000600), the Natural Science Fundamental Research
Program of Shaanxi Province of China under grants 2022JM-021 and the
Fundamental Research Funds for the Central Universities (Grant No.
11913291000022).

## References

  * (1) G. H. Golub, C. F. Van Loan, Matrix Computation, (4-th edition, The Johns Hopkins University Press, 2013).
  * (2) J. W. Demmel, Applied Numerical Linear Algebra, SIAM, (1997).
  * (3) D. Coppersmith and S. Winograd, J. Symb. Comput. 9, 251 (1990).
  * (4) F. Le Gall, in Proceedings of the 39th International Symposium on Symbolic and Algebraic Computation (ACM, New York, 2014), p. 296.
  * (5) A. W. Harrow, A. Hassidim, and S. Lloyd, Quantum Algorithm for Linear Systems of Equations. Phys. Rev. Lett. 103, 150502 (2009).
  * (6) A. M. Childs, R. Kothari, and R. D. Somma, Quantum algorithm for systems of linear equations with exponentially improved dependence on precision. SIAM J. Comput., 46:1920 (2017).
  * (7) Y. Subaşı, R. D. Somma, D. Orsucci, Quantum Algorithms for Systems of Linear Equations Inspired by Adiabatic Quantum Computing. Phys. Rev. Lett. 122, 060504 (2019).
  * (8) L. Lin and Y. Tong, Optimal polynomial based quantum eigenstate filtering with application to solving quantum linear systems, Quantum 4, 361 (2020).
  * (9) G. H. Low and I. L. Chuang, Optimal Hamiltonian Simulation by Quantum Signal Processing, Phys. Rev. Lett. 118, 010501 (2017).
  * (10) H. Wang, Quantum algorithm for obtaining the eigenstates of a physical system, Phys. Rev. A 93, 052334 (2016).
  * (11) Z. Li, X. Liu, H. Wang, S. Ashhab, J. Cui, H. Chen, X. Peng, and J. Du, Quantum Simulation of Resonant Transitions for Solving the Eigenproblem of an Effective Water Hamiltonian, Phys. Rev. Lett. 122, 090504 (2019).
  * (12) C. Cohen-Tannoudji, B. Diu, and F. Laloë, _Quantum Mechanics_ Vol. 1, (Wiley-Interscience Publication 1977).
  * (13) G. H. Low and I. L. Chuang, Hamiltonian simulation by qubitization, Quantum 3, 163 (2019).
  * (14) A. Gilyén, Y. Su, G. H. Low, and N. Wiebe, Quantum singular value transformation and beyond: exponential improvements for quantum matrix arithmetics, Proceedings of the 51st Annual ACM SIGACT Symposium on Theory of Computing (2019).
  * (15) J. M. Martyn, Z. M. Rossi, A. K. Tan and I. L. Chuang, Grand unification of quantum algorithms, PRX Quantum, 2, 040203 (2021).
  * (16) G. H. Low, T. J. Yoder, and I. L. Chuang, Methodology of Resonant Equiangular Composite Quantum Gates, Phys. Rev. X 6, 041067 (2016).
  * (17) J. Haah, Product decomposition of periodic functions in quantum signal processing Quantum 3, 190 (2019).
  * (18) R. Chao, D. Ding, A. Gilyén, C. Huang, and M. Szegedy, Finding angles for quantum signal processing with machine precision, arXiv:2003.02831.
  * (19) F. Arute, et al., Quantum supremacy using a programmable superconducting processor, Nature 574, 505 (2019).
  * (20) Y. Dong, X. Meng, K. B. Whaley, and L. Lin, Efficient phase factor evaluation in quantum signal processing (2020), arXiv:2002.11649.
  * (21) D. W. Berry and A. M. Childs, Black-box Hamiltonian simulation and unitary implementation, Quantum Information & Computation 12, 29 (2012).
  * (22) D. W. Berry, A. M. Childs, and R. Kothari, Hamiltonian simulation with nearly optimal dependence on all parameters. In 2015 IEEE 56th Annual Symposium on Foundations of Computer Science, pages 792–809. IEEE, 2015.
  * (23) H. Wang, S. Yu and H. Xiang, Efficient quantum algorithm for solving structured problems via multistep quantum computation, Phys. Rev. Res. 5, L012004 (2023).


††thanks<EMAIL_ADDRESS>
# Entangling Excitons with Microcavity Photons

Xuan Zuo  Interdisciplinary Center of Quantum Information, State Key
Laboratory of Modern Optical Instrumentation, and Zhejiang Province Key
Laboratory of Quantum Technology and Device, School of Physics, Zhejiang
University, Hangzhou 310027, China  Zhi-Yuan Fan  Interdisciplinary Center of
Quantum Information, State Key Laboratory of Modern Optical Instrumentation,
and Zhejiang Province Key Laboratory of Quantum Technology and Device, School
of Physics, Zhejiang University, Hangzhou 310027, China  Hang Qian
Interdisciplinary Center of Quantum Information, State Key Laboratory of
Modern Optical Instrumentation, and Zhejiang Province Key Laboratory of
Quantum Technology and Device, School of Physics, Zhejiang University,
Hangzhou 310027, China  Jie Li  Interdisciplinary Center of Quantum
Information, State Key Laboratory of Modern Optical Instrumentation, and
Zhejiang Province Key Laboratory of Quantum Technology and Device, School of
Physics, Zhejiang University, Hangzhou 310027, China

###### Abstract

We provide a systemic theory to entangle excitons with microcavity photons.
This is realized by adopting an exciton-optomechanics system and introducing a
nonlinear dispersive interaction with a mechanical oscillator. We show that
when either the exciton and cavity modes in the weak-coupling regime, or the
two exciton-polariton modes in the strong-coupling regime, are respectively
resonant with the optomechanical Stokes and anti-Stokes sidebands,
entanglement between excitons and cavity photons, or between two exciton
polaritons, can be established. The entanglement is in the steady state and
can potentially be achievable at room temperature. In both cases, genuine
tripartite entanglement is shown to be present.

##  I Introduction

An exciton is an electrically neutral quasiparticle that is formed by the
binding of an electron and a hole via the Coulomb interaction. It exists
mainly in condensed matter, e.g., insulators and semiconductors, and can
interact with electromagnetic fields through the exciton-photon dipole
interaction Keeling07 ; Deng10 . It was theoretically predicted that the
strong interaction between excitons and photons can lead to the generation of
exciton polaritons Hopfield58 , which was first experimentally observed in the
semiconductor quantum-well (QW) microcavity Weisbuch92 ; Houdre94 , benefiting
from the significantly improved coupling strength between excitons and photons
due to the strongly confined light in the microcavity. Since then, extensive
research has been made on achieving the exciton-photon strong coupling, e.g.,
by adopting various microcavity structures, including the micropillar cavity,
the photonic crystal slabs, and the whispering gallery microcavity Khitrova06
; Sun08 , and employing different materials, such as organic Lidzey98 ; Kena-
Cohen08 , wide-bandgap Sun08 , and perovskite semiconductors Su21 . Novel
exciton-polariton devices have been designed by exploiting the coherence and
nonlinearity properties of the polaritons, including spin memory Cerna13 ,
polariton light emitting diodes Tsintzos08 ; Chakraborty19 , and polariton
transistors Ballarini13 ; Zasedatelev19 .

As half-matter, half-light bosons, exciton polaritons offer a unique avenue
for exploring the interface between quantum optics, spontaneous coherence, and
quantum condensation Keeling07 . They can be manipulated and probed via the
light component and generate rich nonlinear interactions through the matter
component Keeling07 . Because the effective mass is much lower than that of
typical atomic systems, the critical temperature for achieving the Bose-
Einstein condensation (BEC) of the exciton polaritons can be several orders of
magnitude higher than the atomic one, which is potentially attainable even at
room temperature Deng10 . Despite the highly dissipative two-dimensional
system with weak particle-particle interactions, which do not conform to the
ideal BEC, the exciton-polariton BEC has been successfully achieved in the
experiments BEC1 ; BEC2 , with the observation of relevant macroscopic quantum
phenomena, including Bogoliubov excitations Bogoliubov , quantized vortices
vortices , and superfluidity fluid . The coherent nature of the polariton
condensates can be exploited to create a high-efficiency low-threshold laser
without inversion Imamoglu96 , which was first demonstrated in the experiments
Malpuech02 ; Deng03 . After that, a variety of configurations have been
adopted to realize the room-temperature lasing Christopoulos07 ; Bajoni08 ;
Kena-Cohen10 , and a more efficient polariton laser has been demonstrated by
using an electrical pump method Schneider13 ; Bhattacharya13 ; Bhattacharya14
. In addition, quantum entanglement of the exciton polaritons has been studied
by exploiting the parametric process Ciuti04 ; Liew18 ; DS  and spin squeezing
Feng , and the photon-polariton entanglement has been experimentally observed
by swapping a photon for a polariton in two-photon entangled states Cuevas .
Besides, squeezed light has been produced utilizing the strong nonlinearity of
the polariton-polariton interaction Boulier14 .

Here, we provide a complete theory to entangle excitons with microcavity
photons by coupling the latter to a mechanical oscillator via a nonlinear
radiation-pressure interaction. The system then becomes a tripartite bosonic
system, namely, exciton-optomechanics (EOM) Santos ; EOM1 ; EOM2 ; EOM3 ;
Fainstein20 ; Santos21 ; Santos23 . Due to the linear excitation-exchange
(beam-splitter-type) interaction between excitons and photons, the two modes
are intrinsically not entangled. The mechanical oscillator couples to the
optical cavity mode via a dispersive manner, which can enable an effective
optomechanical parametric down-conversion (PDC) or state-swap interaction,
corresponding to the optomechanical Stokes or anti-Stokes scattering. We
consider a complete set of two situations depending on whether excitons and
microcavity photons are strongly coupled to form exciton polaritons. We show
that when either the exciton and cavity modes (in the weak-coupling regime) or
the two exciton-polariton modes (in the strong-coupling regime) are
respectively resonant with the optomechanical Stokes and anti-Stokes
sidebands, the entanglement between excitons and cavity photons, or between
two polariton modes, can be established. The entanglement is in the steady
state and can be achieved even at room temperature for a not very high
mechanical quality factor. We further show that both the excitons and cavity
photons (or the two polaritons) are entangled with the mechanical oscillator,
exhibiting genuine tripartite entanglement in the EOM system.

The paper is organized as follows. In Sec. II, we describe the EOM system,
provide its Hamiltonian and the corresponding Langevin equations, and show how
the steady-state solutions of the system can be achieved. We then present in
Sec. III the results of the exciton-photon entanglement and the tripartite
entanglement in the weak-coupling case. In Sec. IV, we reformulate the theory
for the exciton-polariton-mechanics system in the strong-coupling regime, and
show the results of the entanglement of two polaritons and the tripartite
entanglement in Sec. V. Finally, we make a discussion on the two entanglement
protocols in the weak- and strong-coupling regimes and summarize the findings
in Sec. VI.

Figure 1: Sketch of the EOM system. A QW is placed within a semiconductor
microcavity that is formed by two DBRs, which enables the interaction between
the exciton mode ($x$) and the optical cavity mode ($c$). The DBRs are movable
(the gray area), and the mechanical motion ($b$) couples to the cavity photons
via radiation pressure.

##  II The exciton-optomechanics system

The EOM system consists of a QW and a semiconductor microcavity formed by two
movable distributed Bragg reflectors (DBRs), as depicted in Fig. 1. A DBR is
made of layers of alternating high and low refraction indices, and each layer
has an optical thickness of $\lambda/4$, with $\lambda$ being the optical
wavelength. Thus, light reflected from each interface destructively interfere
and create a stop band for transmission. The DBR then acts as a high-
reflectivity mirror when the wavelength of the incident light is within the
stop band. The DBRs are movable and the mechanical displacement couples to the
cavity photons via radiation pressure, thus forming an optomechanical system
Perrin13 . A QW is a thin layer of semiconductor, which is sandwiched between
two barrier layers with a much larger band gap. The QW is placed within the
microcavity, which enables a linear excitation-exchange (beam-splitter type)
interaction between excitons and cavity photons. The Hamiltonian of the EOM
system reads

$\displaystyle\begin{split}H/\hbar&=\\!\omega_{x}x^{\dagger}x+\omega_{c}c^{\dagger}c+\omega_{b}b^{\dagger}b+G_{0}c^{\dagger}c\left(b+b^{\dagger}\right)\\\
&+g\left(x^{\dagger}c+xc^{\dagger}\right)+i\Omega\left(c^{\dagger}e^{-i\omega_{0}t}-ce^{i\omega_{0}t}\right),\end{split}$
(1)

where $x$, $c$, and $b$ ($x^{\dagger}$, $c^{\dagger}$, and $b^{\dagger}$) are
the annihilation (creation) operators of the excitons, cavity photons, and
phonons, respectively, satisfying the commutation relation
$[j,j^{{\dagger}}]=1$ ($j=x,\,c,\,b$), and $\omega_{x}$, $\omega_{c}$, and
$\omega_{b}$ are their resonance frequencies. Here, $g$ denotes the exciton-
photon coupling strength, which is variable depending on the position of the
QW placed within the microcavity and can be very strong. When the coupling
strength exceeds the exciton and cavity decay rates $\kappa_{x}$ and
$\kappa_{c}$, i.e., $g>\kappa_{x},\kappa_{c}$, the system enters the strong-
coupling regime leading to the exciton polaritons Weisbuch92 ; Houdre94 . The
single-photon optomechanical coupling strength $G_{0}$ is typically weak MA14
, but the effective optomechanical coupling can be significantly enhanced by
driving the microcavity with an intense laser field. The last term corresponds
to the driving Hamiltonian, where $\Omega=\sqrt{2P\kappa_{c}/\hbar\omega_{0}}$
signifies the coupling strength between the cavity and the drive field with
frequency $\omega_{0}$ and power $P$.

By incorporating the dissipation and input noise of each mode, we obtain the
following quantum Langevin equations (QLEs) in the frame rotating at the drive
frequency $\omega_{0}$:

$\displaystyle\begin{split}\dot{x}=&-(i\Delta_{x}+\kappa_{x})x-igc+\sqrt[]{2\kappa_{x}}x^{in},\\\
\dot{c}=&-(i\Delta_{c}+\kappa_{c})c-igx-
iG_{0}c(b+b^{\dagger})+\Omega+\sqrt[]{2\kappa_{c}}c^{in},\\\
\dot{b}=&-(i\omega_{b}+\kappa_{b})b-iG_{0}c^{\dagger}c+\sqrt[]{2\kappa_{b}}b^{in},\end{split}$
(2)

where $\Delta_{x}=\omega_{x}-\omega_{0}$, $\Delta_{c}=\omega_{c}-\omega_{0}$,
$\kappa_{b}$ is the mechanical damping rate, and $j^{in}(t)$ ($j=x,\,c,\,b$)
are the input noise operators of the three modes, which are zero-mean and
characterized by the correlation functions Zoller : $\langle
j^{in}(t)j^{in\dagger}(t^{\prime})\rangle=[N_{j}(\omega_{j})+1]\delta(t-t^{\prime})$,
$\langle
j^{in\dagger}(t)j^{in}(t^{\prime})\rangle=N_{j}(\omega_{j})\delta(t-t^{\prime})$,
with $N_{j}(\omega_{j})=[\exp[(\hbar\omega_{j}/k_{B}T)]-1]^{-1}$ being the
equilibrium mean thermal excitation number of the mode $j$, and $T$ as the
bath temperature.

Owing to the strong driving of the cavity and the exciton-photon excitation-
exchange interaction, the cavity and exciton modes exhibit large amplitudes
$|\langle c\rangle|,\,|\langle x\rangle|\gg 1$. This allows us to linearize
the system dynamics around the steady-state values by expressing each mode
operator $j$ as the sum of its classical average $\langle j\rangle$ and
quantum fluctuation operator $\delta j$, i.e., $j=\langle j\rangle+\delta j$,
and neglecting small second-order fluctuation terms. As a result, the QLEs are
separated into two sets of equations, one for the classical averages and the
other for the quantum fluctuations. The linearized QLEs of the quantum
fluctuations are obtained as

$\displaystyle\begin{split}\dot{\delta
x}=&-\big{(}i{\Delta}_{x}+\kappa_{x}\big{)}\delta x-ig{\delta
c}+\sqrt[]{2\kappa_{x}}x^{in},\\\ \dot{\delta
c}=&-\big{(}i\tilde{\Delta}_{c}+\kappa_{c}\big{)}\delta c-ig{\delta
x}-G_{cb}\big{(}\delta b+\delta
b^{\dagger}\big{)}+\sqrt[]{2\kappa_{c}}c^{in},\\\ \dot{\delta
b}=&-\big{(}i{\omega}_{b}+\kappa_{b}\big{)}\delta b-\left(G_{cb}\delta
c^{\dagger}-G_{cb}^{*}\delta c\right)+\sqrt[]{2\kappa_{b}}b^{in},\end{split}$
(3)

where $\tilde{\Delta}_{c}={\Delta}_{c}+2G_{0}\langle b\rangle$ is the
effective cavity-drive detuning including the frequency shift due to the
optomechanical interaction, and $G_{cb}=iG_{0}\langle c\rangle$ is the
effective optomechanical coupling strength. The expressions of the classical
averages are given by

$\displaystyle\begin{split}\langle x\rangle=&\frac{-i\Omega
g}{g^{2}+(i\tilde{\Delta}_{c}+\kappa_{c})(i{\Delta}_{x}+\kappa_{x})},\\\
\langle
c\rangle=&\frac{\Omega(i{\Delta}_{x}+\kappa_{x})}{g^{2}+(i\tilde{\Delta}_{c}+\kappa_{c})(i{\Delta}_{x}+\kappa_{x})},\\\
\langle b\rangle=&-\frac{G_{0}}{\omega_{b}}|\langle c\rangle|^{2}.\end{split}$
(4)

The QLEs (3) can be expressed in a compact matrix form with the quadrature
fluctuation operators $\delta X_{j}=(\delta j+\delta
j^{\dagger})/\\!\sqrt{2}$, and $\delta Y_{j}=i(\delta j^{\dagger}-\delta
j)/\\!\sqrt{2}$, i.e.,

$\displaystyle\dot{u}(t)={\cal A}\,u(t)+n(t),$ (5)

where $u(t)=\left[\delta X_{x}(t),\delta Y_{x}(t),\delta X_{c}(t),\delta
Y_{c}(t),\delta X_{b}(t),\delta Y_{b}(t)\right]^{\rm T}$ represents the vector
of quantum fluctuations,
$n(t)=\left[\sqrt[]{2\kappa_{x}}X_{x}^{in},\sqrt[]{2\kappa_{x}}Y_{x}^{in},\sqrt[]{2\kappa_{c}}X_{c}^{in},\sqrt[]{2\kappa_{c}}Y_{c}^{in},\sqrt[]{2\kappa_{b}}X_{b}^{in},\sqrt[]{2\kappa_{b}}Y_{b}^{in}\right]^{\rm
T}$ is the vector of input noises, where $X_{j}^{in}$ and $Y_{j}^{in}$ are
defined similarly as $\delta X_{j}$ and $\delta Y_{j}$ but with the noise
operators $j^{in}$ and $j^{in{\dagger}}$, and the drift matrix ${\cal A}$ is
given by

$\displaystyle\cal A=\begin{pmatrix}-\kappa_{x}&{\Delta}_{x}&0&g&0&0\\\
-{\Delta}_{x}&-\kappa_{x}&-g&0&0&0\\\
0&g&-\kappa_{c}&\tilde{\Delta}_{c}&-2{\rm Re}\,G_{cb}&0\\\
-g&0&-\tilde{\Delta}_{c}&-\kappa_{c}&-2{\rm Im}\,G_{cb}&0\\\
0&0&0&0&-\kappa_{b}&\omega_{b}\\\ 0&0&-2{\rm Im}\,G_{cb}&2{\rm
Re}\,G_{cb}&-\omega_{b}&-\kappa_{b}\\\ \end{pmatrix}.$ (6)

Due to the linearized dynamics and the inherent Gaussian properties of the
quantum noises, the steady state of the quadrature fluctuations is a three-
mode Gaussian state, which is fully characterized by a $6\times 6$ covariance
matrix (CM) $V$, with its entries $V_{ij}=\frac{1}{2}\langle
u_{i}(t)u_{j}(t^{\prime})+u_{j}(t^{\prime})u_{i}(t)\rangle$ $(i,j=1,2,...,6)$.
The steady-state CM can be obtained by directly solving the Lyapunov equation
Parks ; Vitali07

$\displaystyle\begin{split}{\cal A}V+V{\cal A}^{\rm T}=-D,\end{split}$ (7)

where
$D=\mathrm{Diag}\big{[}\kappa_{x}(2N_{x}+1),\kappa_{x}(2N_{x}+1),\kappa_{c}(2N_{c}+1),\kappa_{c}(2N_{c}+1),\kappa_{b}(2N_{b}+1),\kappa_{b}(2N_{b}+1)\big{]}$
is the diffusion matrix, with its entries defined via
$D_{ij}\,\delta(t-t^{\prime})=\langle
n_{i}(t)n_{j}(t^{\prime})+n_{j}(t^{\prime})n_{i}(t)\rangle/2$. Once the CM $V$
is achieved, we adopt the logarithmic negativity $E_{N}$ Ent  to quantify the
quantum entanglement between any two modes of the system and the minimum
residual contangle $R^{\rm min}_{\tau}$ Cont  to calculate the tripartite
entanglement of the system.

Figure 2: (a) The mechanical motion with frequency $\omega_{b}$ scatters the
driving photons at frequency $\omega_{0}$ onto two sidebands at
$\omega_{0}\pm\omega_{b}$. When the exciton and cavity modes are resonant with
the Stokes (red) and anti-Stokes (blue) sidebands, respectively, genuine
tripartite entanglement of cavity photons, excitons, and phonons is present.
(b) Stationary exciton-photon (solid), photon-phonon (dashed), and exciton-
phonon (dot-dashed) entanglement and tripartite entanglement (inset) versus
the exciton-drive detuning $\Delta_{x}$. We take
$\tilde{\Delta}_{c}=0.9\omega_{b}$, and the other parameters are provided in
the text.

##  III Entanglement between Excitons and Microcavity Photons

In this section, we consider the weak-coupling case where the interaction
between excitons and microcavity photons does not lead to exciton polaritons.
Specifically, we consider a relatively weak exciton-photon coupling
$g<\kappa_{c},\kappa_{x}$, and the exciton and cavity modes are non-resonant
but of close frequencies, cf. Fig. 2(a).

The prerequisite for obtaining entanglement in the system is to cool the
lower-frequency mechanical mode close to its quantum ground state and
eliminate the detrimental thermal noise. To this end, we adopt a red-detuned
laser field to drive the microcavity (with
$\tilde{\Delta}_{c}\approx\omega_{b}\gg\kappa_{c}$) to activate the
optomechanical anti-Stokes scattering, which is responsible for cooling the
mechanical motion. For mechanical cooling, the drive field should be
relatively weak, such that the rotating-wave approximation (under the weak-
coupling condition $G_{cb}\ll\omega_{b}$) can be taken to obtain the effective
optomechanical beam-splitter interaction
$\propto\delta{c}^{\dagger}\delta{b}+\delta{c}\delta{b}^{\dagger}$ DV08 .
However, the optomechanical beam-splitter interaction, as well as the linear
exciton-photon interaction, does not produce any entanglement. To create
entanglement, we increase the drive power, such that the weak-coupling
condition is broken and the counter-rotating-wave (CRW) terms
$\propto\delta{c}^{\dagger}\delta{b}^{\dagger}+\delta{c}\delta{b}$ in the
linearized optomechanical interaction cannot be neglected, which correspond to
the PDC interaction yielding the optomechanical entanglement Vitali07 . The
entanglement is further distributed to the exciton-phonon and exciton-photon
subsystems when the exciton mode resonates with the mechanical Stokes
sideband, i.e., $-\Delta_{x}\approx\omega_{b}\gg\kappa_{x}$ (Fig. 2(a)) Li18 ;
Zuo . In this situation, i.e.,
$\tilde{\Delta}_{c}\approx-\Delta_{x}\approx\omega_{b}$, all bipartite
subsystems are entangled, as shown in Fig 2(b), and the EOM system shares
genuine tripartite entanglement, as witnessed by a nonzero $R^{\rm
min}_{\tau}$ (inset of Fig 2(b)). The complementary relation of the dashed
curve and the solid and dot-dashed curves in Fig 2(b) indicates that the
exciton-phonon and exciton-photon entanglement are transferred from the
photon-phonon entanglement. Similar entanglement distribution phenomena in
multipartite systems have been observed in cavity magnomechanics Li18  and
optomagnomechanics Li23 .

The exciton-photon entanglement can be understood straightforwardly as
follows. The mechanical motion scatters the driving photons at frequency
$\omega_{0}$ onto the two mechanical sidebands at frequencies
$\omega_{0}\pm\omega_{b}$. When the exciton and cavity modes are resonant with
the two sidebands, i.e.,
$\tilde{\Delta}_{c}\approx-\Delta_{x}\approx\omega_{b}$, the optomechanical
Stokes and anti-Stokes scatterings are simultaneously activated and enhanced
(here Stokes photons and excitons are of the excitation-exchange interaction).
The Stokes scattering corresponds to the PDC interaction, leading the Stokes
photons to be entangled with the mechanical oscillator, whereas the anti-
Stokes scattering results in the state-swap (beam-splitter) interaction
between the anti-Stokes photons and the oscillator. Therefore, the two
sidebands become entangled via the mediation of the mechanical motion. Since
the exciton and cavity modes resonate with the two sidebands, respectively,
the exciton and cavity modes thus get entangled.

In getting Fig 2(b), we employ the following feasible parameters Santos ; EOM1
; Perrin13 : $\omega_{b}/2\pi=20$ GHz, $\kappa_{b}/2\pi=1$ MHz,
$\kappa_{c}/2\pi=1$ GHz, $\kappa_{x}/2\pi=10^{2}$ MHz, $G_{0}/2\pi=10$ MHz,
$\Omega/2\pi=6$ THz (corresponding to the drive power $P\,\,{\approx}\,\,26$
mW for $\omega_{0}/2\pi\,\,{\approx}\,\,345$ THz), and at temperature $T=1$ K.
We take a moderate coupling $g/2\pi=0.9$ GHz, and the exciton and cavity
frequencies are nearly resonant and on both sides of the drive frequency.
Under these parameters, the mechanical motion is cooled to its quantum ground
state with the mean phonon number of $0.03$ at the optimal condition
$\tilde{\Delta}_{c}\approx-\Delta_{x}\approx\omega_{b}$. In Fig 3, we study
the impact of various dissipation rates of the system and the bath temperature
on the exciton-photon entanglement. Clearly, the entanglement is robust
against all dissipation rates and is still present for $\kappa_{c}$,
$\kappa_{x}$, $\kappa_{b}$ being up to $\sim$10 GHz, 1 GHz, and 10 MHz,
respectively, based on the parameters of Fig 2(b). Due to the relatively high
mechanical frequency in the typical EOM system, the entanglement survives at a
temperature up to $\sim 60$ K.

Figure 3: Stationary exciton-photon entanglement versus (a) dissipation rates
$\kappa_{x}$ and $\kappa_{c}$; (b) bath temperature $T$ and mechanical damping
rate $\kappa_{b}$. We take $\Delta_{x}=-1.1\omega_{b}$, and the other
parameters are the same as in Fig 2(b).

##  IV The Polariton-Mechanics system

When the exciton-photon system enters the strong-coupling regime,
$g>\kappa_{c},\kappa_{x}$, excitons and cavity photons are hybridized forming
two exciton polaritons, which are termed as the upper polariton (UP) and the
lower polariton (LP), respectively. The theory presented in Secs. II and III
is based on the original exciton and cavity modes and thus becomes inefficient
for the study of the entanglement between two polaritons. In this case, it is
more convenient to express the theory in terms of two polaritons.
Consequently, we rewrite the EOM Hamiltonian (1) with the polariton operators

$\displaystyle\begin{split}H/\hbar&=\\!\omega_{u}U^{\dagger}U+\omega_{l}L^{\dagger}L+\omega_{b}b^{\dagger}b\\\
&+G_{0}\left(b+b^{\dagger}\right)\Big{(}U^{\dagger}U\sin^{2}\theta+U^{\dagger}L\sin\theta\cos\theta\\\
&\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,+L^{\dagger}U\cos\theta\sin\theta+L^{\dagger}L\cos^{2}\theta\Big{)}\\\
&+i\Omega\Big{(}U^{\dagger}\sin\theta e^{-i\omega_{0}t}+L^{\dagger}\cos\theta
e^{-i\omega_{0}t}\\\ &\,\,\,\,\,\,\,\,\,\,\,\,\,\,-U\sin\theta
e^{i\omega_{0}t}-L\cos\theta e^{i\omega_{0}t}\Big{)}\,,\end{split}$ (8)

where $U$ and $L$ ($U^{\dagger}$ and $L^{\dagger}$) denote the annihilation
(creation) operators of the UP and LP, respectively, which are the
hybridization of the exciton and cavity modes via the unitary Hopfield
transformation Hopfield58 : $U=x\cos\theta+c\sin\theta$ and
$L=-x\sin\theta+c\cos\theta$, with the corresponding eigenfrequencies

$\displaystyle\begin{split}\omega_{u}&=\frac{1}{2}\Big{[}\omega_{x}+\omega_{c}+\sqrt{(\omega_{x}-\omega_{c})^{2}+4g^{2}}\Big{]},\\\
\omega_{l}&=\frac{1}{2}\Big{[}\omega_{x}+\omega_{c}-\sqrt{(\omega_{x}-\omega_{c})^{2}+4g^{2}}\Big{]}.\end{split}$
(9)

The proportions of excitons and cavity photons in the two polaritons are
characterized by the Hopfield coefficients $\sin^{2}\theta$ and
$\cos^{2}\theta$, where
$\theta=\frac{1}{2}\arctan\frac{2g}{\omega_{x}-\omega_{c}}\in[0,\frac{\pi}{2}]$,
implying that the proportions can be varied by altering the coupling $g$
and/or the exciton-photon detuning $\omega_{x}-\omega_{c}$. Therefore, the
mechanical mode couples to both the polaritons because both of them contain
the photon component, as shown in the Hamiltonian (8). When the cavity photons
and excitons are resonant, yielding $\theta=\frac{\pi}{4}$, both the UP and LP
are exactly half-photon, half-exciton with the minimum frequency splitting
$\omega_{u}-\omega_{l}=2g$. Since both photons and excitons are bosons, the
exciton polaritons, as their linear superposition, are also bosons, and
satisfy the bosonic commutation relation $[k,k^{\dagger}]=1$ ($k=U,L$).

Likewise, we obtain the following QLEs, in the frame rotating at the drive
frequency, concerning two polariton modes and a mechanical mode:

$\displaystyle\begin{split}\dot{U}=&-i\Delta_{u}U-iG_{0}\left(b+b^{\dagger}\right)\left(U\sin^{2}\theta+L\sin\theta\cos\theta\right)\\\
&-\kappa_{u}U-\delta\kappa L+\Omega\sin\theta+\sqrt[]{2\kappa_{u}}U^{in},\\\
\dot{L}=&-i\Delta_{l}L-iG_{0}\left(b+b^{\dagger}\right)\left(L\cos^{2}\theta+U\cos\theta\sin\theta\right)\\\
&-\kappa_{l}L-\delta\kappa U+\Omega\cos\theta+\sqrt[]{2\kappa_{l}}L^{in},\\\
\dot{b}=&-i\omega_{b}b-iG_{0}\Big{(}U^{\dagger}U\sin^{2}\theta+U^{\dagger}L\sin\theta\cos\theta\\\
&+L^{\dagger}U\cos\theta\sin\theta+L^{\dagger}L\cos^{2}\theta\Big{)}-\kappa_{b}b+\sqrt[]{2\kappa_{b}}b^{in},\end{split}$
(10)

where $\Delta_{u}=\omega_{u}-\omega_{0}$ and
$\Delta_{l}=\omega_{l}-\omega_{0}$ denote the polariton-drive detunings;
$\kappa_{u}=\kappa_{x}\cos^{2}\theta+\kappa_{c}\sin^{2}\theta$ and
$\kappa_{l}=\kappa_{x}\sin^{2}\theta+\kappa_{c}\cos^{2}\theta$ are the
dissipation rates of the two polariton modes; and
$\delta\kappa\equiv(\kappa_{c}-\kappa_{x})\sin\theta\cos\theta$ signifies the
dissipative coupling between the polaritons stemming from the unbalanced decay
rates of the original modes, i.e., $\kappa_{x}\neq\kappa_{c}$.
$U^{in}\equiv(\sqrt[]{2\kappa_{x}}\cos\theta
x^{in}+\sqrt[]{2\kappa_{c}}\sin\theta c^{in})/\sqrt[]{2\kappa_{u}}$ and
$L^{in}\equiv(-\sqrt[]{2\kappa_{x}}\sin\theta
x^{in}+\sqrt[]{2\kappa_{c}}\cos\theta c^{in})/\sqrt[]{2\kappa_{l}}$ represent
the input noises entering the UP and LP, which are the combination of the
input noises $x^{in}$ and $c^{in}$ of the exciton and cavity modes.

Similarly, the linearization of the system dynamics around the steady-state
values leads to the following QLEs for the quantum fluctuations of the system:

$\displaystyle\begin{split}\dot{\delta
U}{=}&\,{-}\big{(}i\tilde{\Delta}_{u}+\kappa_{u}\big{)}\delta
U{-}\,\big{(}iG_{b}+\delta\kappa\big{)}\delta L{-}\,G_{u,b}\big{(}\delta
b+\delta b^{\dagger}\big{)}{+}\,\sqrt[]{2\kappa_{u}}U^{in},\\\ \dot{\delta
L}{=}&\,{-}\big{(}i\tilde{\Delta}_{l}+\kappa_{l}\big{)}\delta
L{-}\,\big{(}iG_{b}+\delta\kappa\big{)}\delta U{-}\,G_{l,b}\big{(}\delta
b+\delta b^{\dagger}\big{)}{+}\,\sqrt[]{2\kappa_{l}}L^{in},\\\ \dot{\delta
b}\,{=}&\,{-}\big{(}i\omega_{b}+\kappa_{b}\big{)}\delta b\\\
&\,{-}\Big{(}G_{u,b}\delta U^{\dagger}+G_{l,b}\delta
L^{\dagger}-G_{u,b}^{*}\delta U-G_{l,b}^{*}\delta
L\Big{)}+\sqrt[]{2\kappa_{b}}b^{in},\end{split}$ (11)

where $\tilde{\Delta}_{u}=\Delta_{u}+2G_{0}\langle b\rangle\sin^{2}\theta$ and
$\tilde{\Delta}_{l}=\Delta_{l}+2G_{0}\langle b\rangle\cos^{2}\theta$ are the
effective polariton-drive detunings, including the frequency shift caused by
the optomechanical interaction;
$G_{u,b}=\left(G_{u}\sin\theta+G_{l}\cos\theta\right)\sin\theta$ and
$G_{l,b}=\left(G_{u}\sin\theta+G_{l}\cos\theta\right)\cos\theta$ are the
coupling strengths between the two polaritons and the mechanical mode,
respectively, where $G_{u}=iG_{0}\langle U\rangle$ and $G_{l}=iG_{0}\langle
L\rangle$; and $G_{b}=G_{0}\langle b\rangle\sin 2\theta$ denotes the coupling
between the two polaritons via the mediation of the mechanical mode. The
steady-state averages are given by

$\displaystyle\begin{split}\langle
U\rangle=&\frac{\Omega\left[\delta\kappa\cos\theta-i\sin\theta\left(\tilde{\Delta}_{l}-2G_{0}\langle
b\rangle\cos^{2}\theta-i\kappa_{l}\right)\right]}{\big{(}\tilde{\Delta}_{l}-i\kappa_{l}\big{)}\big{(}\tilde{\Delta}_{u}-i\kappa_{u}\big{)}+\delta\kappa^{2}-G_{b}\,\big{(}G_{b}-2i\delta\kappa\big{)}},\\\
\langle
L\rangle=&\frac{\Omega\left[\delta\kappa\sin\theta-i\cos\theta\left(\tilde{\Delta}_{u}-2G_{0}\langle
b\rangle\sin^{2}\theta-i\kappa_{u}\right)\right]}{\big{(}\tilde{\Delta}_{l}-i\kappa_{l}\big{)}\big{(}\tilde{\Delta}_{u}-i\kappa_{u}\big{)}+\delta\kappa^{2}-G_{b}\,\big{(}G_{b}-2i\delta\kappa\big{)}},\\\
\langle b\rangle=&-\frac{G_{0}}{\omega_{b}}\big{|}\langle
U\rangle\sin\theta+\langle L\rangle\cos\theta\big{|}^{2}.\end{split}$ (12)

Note that since the single-photon optomechanical coupling $G_{0}$ is typically
small MA14  and $G_{0}\langle b\rangle\propto G_{0}^{2}$, one can safely
neglect the weak coupling terms $G_{b}\left(\delta U^{\dagger}\delta L+\delta
U\delta L^{\dagger}\right)$ in the QLEs (11) and the last term of the
denominator of $\langle U\rangle$ ($\langle L\rangle$) in Eq. (12), and assume
the effective detunings to be $|\tilde{\Delta}_{u,l}|\simeq|\Delta_{u,l}|$, as
the optimal detunings for entanglement correspond to
$|\Delta_{u,l}|\simeq\omega_{b}$, which will be shown later.

In the same way, the QLEs (11) can be rewritten in terms of the quadrature
fluctuations $\left(\delta X_{u},\delta Y_{u},\delta X_{l},\delta Y_{l},\delta
X_{b},\delta Y_{b}\right)$, which can be cast in the matrix form similarly as
in Eq. (5). The steady-state CM $V^{\prime}$, with its entries defined as
$V^{\prime}_{ij}=\frac{1}{2}\langle
u^{\prime}_{i}(t)u^{\prime}_{j}(t^{\prime})+u^{\prime}_{j}(t^{\prime})u^{\prime}_{i}(t)\rangle$,
where $u^{\prime}(t)=\left[\delta X_{u}(t),\delta Y_{u}(t),\delta
X_{l}(t),\delta Y_{l}(t),\delta X_{b}(t),\delta Y_{b}(t)\right]^{\rm T}$, can
be achieved by solving the Lyapunov equation

$\displaystyle\begin{split}{\cal A^{\prime}}V^{\prime}+V^{\prime}{\cal
A^{\prime}}^{\rm T}=-D^{\prime},\end{split}$ (13)

where the diffusion matrix
$D^{\prime}=\mathrm{Diag}\big{[}\kappa_{u}(2N_{u}{+}\,1),\kappa_{u}(2N_{u}{+}\,1),$
$\kappa_{l}(2N_{l}+1),\kappa_{l}(2N_{l}+1),\kappa_{b}(2N_{b}+1),\kappa_{b}(2N_{b}+1)\big{]}+\frac{1}{2}\tan
2\theta$
$\left[-\kappa_{u}(2N_{u}{+}\,1)+\kappa_{l}(2N_{l}{+}\,1)\right]\bm{\sigma}^{x}\,\otimes\,\rm{I}_{2\times
2}\,\oplus\,\bm{0}_{2\times 2}$, with $\bm{\sigma}^{x}$ being the $x$-Pauli
matrix, and
$N_{u}=\frac{1}{2}\big{\\{}[\kappa_{x}\cos^{2}\theta(2N_{x}+1)+\kappa_{c}\sin^{2}\theta(2N_{c}+1)]/\kappa_{u}-1\big{\\}}$
and
$N_{l}=\frac{1}{2}\big{\\{}[\kappa_{x}\sin^{2}\theta(2N_{x}+1)+\kappa_{c}\cos^{2}\theta(2N_{c}+1)]/\kappa_{l}-1\big{\\}}$
as the mean thermal excitation numbers of the UP and LP, respectively. The
drift matrix ${\cal A^{\prime}}$ is given by

$\displaystyle\cal
A^{\prime}=\begin{pmatrix}-\kappa_{u}&{\Delta}_{u}&-\delta\kappa&0&-2{\rm
Re}\,G_{u,b}&0\\\ -{\Delta}_{u}&-\kappa_{u}&0&-\delta\kappa&-2{\rm
Im}\,G_{u,b}&0\\\ -\delta\kappa&0&-\kappa_{l}&{\Delta}_{l}&-2{\rm
Re}\,G_{l,b}&0\\\ 0&-\delta\kappa&-{\Delta}_{l}&-\kappa_{l}&-2{\rm
Im}\,G_{l,b}&0\\\ 0&0&0&0&-\kappa_{b}&\omega_{b}\\\ {-2{\rm
Im}\,G_{u,b}}&2{\rm Re}\,G_{u,b}&{-2{\rm Im}\,G_{l,b}}&2{\rm
Re}\,G_{l,b}&-\omega_{b}&-\kappa_{b}\\\ \end{pmatrix}.$ (14)

With the CM $V^{\prime}$ in hand, we can then calculate the bipartite
(tripartite) entanglement in the system using the logarithmic negativity (the
minimum residual contangle) introduced in Sec. II.

##  V Entanglement between two Exciton Polaritons

Figure 4: (a) Strongly-coupled excitons and cavity photons form two exciton
polaritons at frequencies $\omega_{u}$ and $\omega_{l}$. When the UP resonates
with the anti-Stokes (blue) sideband and the LP resonates with the Stokes
(red) sideband, the system exhibits genuine tripartite entanglement. (b)
Stationary UP-LP (solid), LP-phonon (dot-dashed), and UP-phonon (dashed)
entanglement and tripartite (inset) entanglement versus detuning
$\Delta_{u}=-\Delta_{l}$. See text for the other parameters.

As analyzed in Sec. III, the two optomechanical sidebands at
$\omega_{0}\pm\omega_{b}$ are entangled due to the mediation of the mechanical
oscillator involved in both the Stokes and anti-Stokes scatterings. It is
quite natural to conjecture that the entanglement between the two polaritons
can be achieved when they are respectively resonant with the two sidebands,
i.e., $\Delta_{u}=-\Delta_{l}=\omega_{b}$ (Fig 4(a)). Since both the
polaritons contain the photon component and their exciton component has no
coupling with the mechanical mode, the interaction between each polariton and
the mechanics is fully determined by the interaction between its photon
component and the mechanics, i.e., the dispersive optomechanical interaction
Zuo . The UP (LP) being resonant with the anti-Stokes (Stokes) sideband
greatly enhances the strength of the anti-Stokes (Stokes) scattering and thus
the optomechanical cooling (PDC) interaction. A unique feature and advantage
of the polariton system is that the strength of the cooling (PDC) interaction
associated with the UP (LP) is adjustable by changing the weight of the photon
component in the polaritons via altering $\theta$. As also discussed in Sec.
III, the emergence of entanglement in the system requires the combination of
both the cooling and PDC interactions for eliminating thermal noise and
creating entanglement, respectively. The fact that the cavity photons enter
the two polaritons with the weights of $\sin^{2}\theta$ and $\cos^{2}\theta$,
respectively, implies that there is a trade-off between the cooling and PDC
interactions for achieving the maximal entanglement, which indicates an
optimal value of $\theta$.

To reach such an optimal $\theta$ in a real experiment, we fix the exciton
frequency and the exciton-photon coupling, e.g., at $\omega_{x}/2\pi=345$ THz
and $g/2\pi=13$ GHz ($g<\omega_{b}$), but set the cavity frequency
$\omega_{c}$ as a variable, which can be realized by adopting the
configuration where the microcavity layer thickness is tapered by growth and
thus the cavity resonance can be continuously tuned across the sample Deng10 ;
Weisbuch92 . This leads to a continuous varying of $\theta$, since
$\theta=\frac{1}{2}\arctan\frac{2g}{\omega_{x}-\omega_{c}}$. It should be
noted that as the cavity frequency changes, the two polariton-drive detunings
$\Delta_{u}$ and $|\Delta_{l}|$ change, but they are generally not equal (cf.,
Eq. (9)). To maintain the optimal condition
$\Delta_{u}=-\Delta_{l}\approx\omega_{b}$ as the cavity frequency varies, we
set the drive frequency to be $\omega_{c}$-dependent, i.e.,
$\omega_{0}=\frac{\omega_{x}+\omega_{c}}{2}$.

In Fig. 4(b), we present the results of the bipartite and tripartite
entanglement of the system. As expected, the entanglement between the two
polaritons is maximized, $E_{N}\approx 0.28$, at the detunings
$\Delta_{u}=-\Delta_{l}\approx\omega_{b}$, and the LP-phonon entanglement is
also optimized at this point, because the LP-phonon entanglement is a direct
result of the Stokes scattering, of which the strength is maximal when the LP
resonates with the Stokes sideband. However, at this point the UP-phonon
entanglement is absent. This is because the effective interaction between the
UP and the mechanical mode is solely the cooling (beam-splitter) interaction,
which does not yield entanglement. One can increase the UP-mechanics coupling
strength by raising the weight of the photon component in the UP, such that,
as analyzed in Sec. III, the weak-coupling condition is broken and the CRW
terms start to play the role in creating the UP-phonon entanglement. The
increase of $\Delta_{u}$ in Fig. 4(b) corresponds to a continuous rise of
$\omega_{c}$, which leads to a continuous increase of $\theta$, starting from
$\theta=\frac{\pi}{4}$ corresponding to $\omega_{c}=\omega_{x}$. This results
in a growing weight of the photon component in the UP and thus an increasing
UP-mechanics coupling strength. When $\Delta_{u}>\approx 1.2\omega_{b}$, the
coupling strength becomes sufficiently strong, leading to the emergence of the
UP-phonon entanglement. A much larger detuning $\Delta_{u}$ causes the
polaritons to deviate from the mechanical sidebands, thus diminishing the
entanglement. The system also exhibits genuine tripartite entanglement around
$\Delta_{u}=-\Delta_{l}\approx\omega_{b}$, as shown in the inset of Fig 4(b).
In getting Fig. 4(b), we keep the coupling strength $G_{l}$ fixed at
$|G_{l}|/2\pi=0.6$ GHz by adjusting the drive power as the detuning
$\Delta_{u}$ varies (cf. Eq. (12)). The corresponding $\Omega/2\pi=3.5$ THz
and drive power $P\approx 8.6$ mW at the optimal detunings
$\Delta_{u}=-\Delta_{l}=\omega_{b}$. The other parameters are the same as
those in Fig 2.

Figure 5: Stationary entanglement between two polaritons versus (a)
$\kappa_{x}$ and $\kappa_{c}$; (b) bath temperature $T$ and $\kappa_{b}$. We
take the optimal detunings $\Delta_{u}=-\Delta_{l}=\omega_{b}$. The other
parameters are the same as in Fig 4(b).

The entanglement of two polaritons is much stronger than the exciton-photon
entanglement generated in Sec. III even with a smaller drive power. The
entanglement is stationary and robust with respect to the dissipation rates of
the three modes, as seen in Fig. 5. The entanglement can be achieved even at
room temperature if using a sub-MHz mechanical damping rate. For example, the
entanglement is still present for the temperature up to $300$ K with
$\kappa_{b}/2\pi=0.5$ MHz (Fig. 5(b)), corresponding to a mechanical $Q$
factor of $4\times 10^{4}$, which is high but reachable Perrin13 ; Santos .

##  VI Discussion and Conclusion

In view of the presented two entanglement protocols in the weak- and strong-
coupling regimes, they share similarities but also have their own
characteristics. In both situations, the generation of entanglement requires
the combination of the cooling and PDC interactions. For the cooling
interaction, in both cases it is achieved by activating the optomechanical
anti-Stokes scattering, realized by keeping the anti-Stokes sideband resonant
with the cavity mode in the weak-coupling regime (with the UP in the strong-
coupling regime). However, their entanglement mechanisms are not exactly the
same. In entangling excitons and cavity photons, the power of the red-detuned
drive field must be sufficiently strong to activate the CRW terms to have the
PDC interaction (Sec. III); whereas in the strong-coupling case of entangling
two polaritons, the PDC interaction can be directly provided by the
optomechanical Stokes scattering associated with the LP (Sec. V), which does
not require a strong drive field. The fact that the cavity photons enter the
two polaritons to simultaneously participate in both the optomechanical Stokes
and anti-Stokes scatterings makes the protocol in the strong-coupling case
more efficient, reflected by the fact that the entanglement between the two
polaritons is much stronger than the exciton-photon entanglement even under a
weaker drive field.

The generated bipartite and tripartite entanglement can be verified by
measuring the corresponding CMs Vitali07 ; Li18 . For the exciton-photon
entanglement, the cavity field quadratures can be measured directly by
homodyning the cavity output field. The excitonic quadratures can be detected
by sending a weak probe field that realizes a state-swap interaction with the
excitons, and by homodyning the cavity output field of the probe field. Note
that in typical exciton-microcavity systems, the cavity decay rate is much
larger, such that when the laser drive is switched off and all cavity photons
dissipate, the excitonic state remains practically unchanged, at which time
the probe field is sent. Similarly, to measure the entanglement of two
polaritons, one can successively send two weak probe fields that are
respectively resonant with the two polaritons, and perform the homodyne
detection of the cavity output fields. To verify the tripartite entanglement,
the mechanical state must be accessed. This can be realized by sending a weak
red-detuned light into the cavity, which activates an effective optomechanical
state-swap interaction and maps the mechanical state to the cavity field, and
by homodying the cavity output field. Due to the much longer mechanical
coherence time, the probe light can be sent after all cavity photons and
excitons die out.

In conclusion, we present a systemic theory for entangling excitons and
microcavity photons, or two exciton polaritons when they are strongly coupled.
The idea is to introduce a dispersively coupled mechanical mode into the
exciton-photon system, which brings in the optomechanical cooling and PDC
interactions, responsible for eliminating thermal noise and creating
entanglement in the system, respectively. By appropriately adjusting the
strengths of the two interactions, stationary exciton-photon or polariton-
polariton entanglement can be achieved. Impressively, room-temperature
polariton entanglement can potentially be obtained by improving relevant
experimental parameters of the EOM system.

Although many classical phenomena have been studied in the EOM system EOM1 ;
EOM2 ; Fainstein20 ; Santos21 ; Santos23 , quantum effects have been rarely
explored Santos ; EOM3 . This work represents the first entanglement study, to
our knowledge, in the field of EOM. The work also provides theoretical
guidance for the experimental realization of the exciton-photon entanglement
and the polariton entanglement, which may find applications in quantum
information processing with exciton polaritons, e.g., entangled polaritons can
lead to the emission of frequency-entangled photon pairs Ciuti04 .

## ACKNOWLEDGMENTS

This work has been supported by National Key Research and Development Program
of China (Grant no. 2022YFA1405200) and National Natural Science Foundation of
China (Grant no. 92265202).

## References

  * (1) H. Deng, H. Haug, and Y. Yamamoto, Exciton-polariton Bose-Einstein condensation, Rev. Mod. Phys. 82, 1489 (2010).
  * (2) J Keeling, F M Marchetti, M H Szymańska, and P B Littlewood, Collective coherence in planar semiconductor microcavities, Semicond. Sci. Technol. 22, R1 (2007).
  * (3) J. J. Hopfield, Theory of the contribution of excitons to the complex dielectric constant of crystals, Phys. Rev. 112, 1555 (1958).
  * (4) C. Weisbuch, M. Nishioka, A. Ishikawa, and Y. Arakawa, Observation of the coupled exciton-photon mode splitting in a semiconductor quantum microcavity, Phys. Rev. Lett. 69, 3314 (1992).
  * (5) R. Houdré, C. Weisbuch, R. P. Stanley, U. Oesterle, P. Pellandini, and M. Ilegems, Measurement of cavity-polariton dispersion curve from angle-resolved photoluminescence experiments, Phys. Rev. Lett. 73, 2043 (1994).
  * (6) G. Khitrova, H. M. Gibbs, M. Kira, S. W. Koch, and A. Scherer, Vacuum Rabi splitting in semiconductors, Nat. Phys. 2, 81 (2006).
  * (7) L. Sun, Z. Chen, Q. Ren, K. Yu, L. Bai, W. Zhou, H. Xiong, Z. Q. Zhu, and X. Shen, Direct observation of whispering gallery mode polaritons and their dispersion in a ZnO tapered microcavity, Phys. Rev. Lett. 100, 156403 (2008).
  * (8) D. G. Lidzey, D. D. C. Bradley, M. S. Skolnick, T. Virgili, S. Walker, and D. M. Whittaker, Strong exciton-photon coupling in an organic semiconductor microcavity, Nature 395, 53 (1998).
  * (9) S Kéna-Cohen, M. Davanço, and S. R. Forrest, Strong exciton-photon coupling in an organic single crystal microcavity, Phys. Rev. Lett. 101, 116401 (2008).
  * (10) R. Su, A. Fieramosca, Q. Zhang, H. S. Nguyen, E. Deleporte, Z. Chen, D. Sanvitto, T. C. H. Liew, and Q. Xiong, Perovskite semiconductors for room-temperature exciton-polaritonics, Nat. Mater. 20, 1315 (2021).
  * (11) R. Cerna, Y. Léger, T. K. Paraïso, M. Wouters, F. Morier-Genoud, M. T. Portella-Oberli, and B. Deveaud, Ultrafast tristable spin memory of a coherent polariton gas, Nat. Commun. 4, 2008 (2013).
  * (12) S. I. Tsintzos, N. T. Pelekanos, G. Konstantinidis, Z. Hatzopoulos, and P. G. Savvidis, A GaAs polariton light-emitting diode operating near room temperature, Nature 453, 372 (2008).
  * (13) J. Gu, B. Chakraborty, M. Khatoniar, and V. M. Menon, A room-temperature polariton light-emitting diode based on monolayer WS2, Nat. Nanotechnol. 14, 1024 (2019).
  * (14) D. Ballarini, M. D. Giorgi, E. Cancellieri, R. Houdré, E. Giacobino, R. Cingolani, A. Bramati, G. Gigli, and D. Sanvitto, All-optical polariton transistor, Nat. Commun. 4, 1778 (2013).
  * (15) A. V. Zasedatelev, A. V. Baranikov, D. Urbonas, F. Scafirimuto, U. Scherf, T. Stöferle, R. F. Mahrt, and P. G. Lagoudakis, A room-temperature organic polariton transistor, Nat. Photonics 13, 378 (2019).
  * (16) J. Kasprzak, M. Richard, S. Kundermann, A. Baas, P. Jeambrun, J. M. J. Keeling, F. M. Marchetti, M. H. Szymańska, R. André, J. L. Staehli, V. Savona, P. B. Littlewood, B. Deveaud, and Le Si Dang, Bose–Einstein condensation of exciton polaritons, Nature 443, 409 (2006).
  * (17) R. Balili, V. Hartwell, D. Snoke, L. Pfeiffer, and K. West, Bose-Einstein condensation of microcavity polaritons in a trap, Science 316, 1007 (2007).
  * (18) S. Utsunomiya, L. Tian, G. Roumpos, C. W. Lai, N. Kumada, T. Fujisawa, M. Kuwata-Gonokami, A. Löffler, S. Höfling, A. Forchel, and Y. Yamamoto, Observation of Bogoliubov excitations in exciton-polariton condensates, Nat. Phys. 4, 700 (2008).
  * (19) K. G. Lagoudakis, M. Wouters, M. Richard, A. Baas, I. Carusotto, R. André, Le Si Dang, and B. Deveaud-Plédran, Quantized vortices in an exciton–polariton condensate, Nat. Phys. 4, 706 (2008).
  * (20) A. Amo, D. Sanvitto, F. P. Laussy, D. Ballarini, E. del Valle, M. D. Martin, A. Lemaître, J. Bloch, D. N. Krizhanovskii, M. S. Skolnick, C. Tejedor, and L. Viña, Collective fluid dynamics of a polariton condensate in a semiconductor microcavity, Nature 457, 291 (2009).
  * (21) A. Imamoglu, R. J. Ram, S. Pau, and Y. Yamamoto, Nonequilibrium condensates and lasers without inversion: Exciton-polariton lasers, Phys. Rev. A 53, 4250 (1996).
  * (22) G. Malpuech, A. D. Carlo, A. Kavokin, J. J. Baumberg, M. Zamfirescu, and P. Lugli, Room-temperature polariton lasers based on GaN microcavities, Appl. Phys. Lett. 81, 412 (2002).
  * (23) H. Deng, G. Weihs, D. Snoke, J. Bloch, and Y. Yamamoto, Polariton lasing vs. photon lasing in a semiconductor microcavity, Proc. Natl. Acad. Sci. 100, 15318 (2003).
  * (24) S. Christopoulos, G. Baldassarri Höger von Högersthal, A. J. D. Grundy, P. G. Lagoudakis, A. V. Kavokin, J. J. Baumberg, G. Christmann, R. Butté, E. Feltin, J.-F. Carlin, and N. Grandjean, Room-temperature polariton lasing in semiconductor microcavities, Phys. Rev. Lett. 98, 126405 (2007).
  * (25) D. Bajoni, P. Senellart, E. Wertz, I. Sagnes, A. Miard, A. Lemaître, and J. Bloch, Polariton laser using single micropillar GaAs-GaAlAs semiconductor cavities, Phys. Rev. Lett. 100, 047401 (2008).
  * (26) S. Kéna-Cohen and S. R. Forrest, Room-temperature polariton lasing in an organic single-crystal microcavity, Nat. Photonics 4, 371 (2010).
  * (27) C. Schneider, A. Rahimi-Iman, N. Y. Kim, J. Fischer, I. G. Savenko, M. Amthor, M. Lermer, A. Wolf, L. Worschech, V. D. Kulakovskii, I. A. Shelykh, M. Kamp, S. Reitzenstein, A. Forchel, Y. Yamamoto, and S. Höfling, An electrically pumped polariton laser, Nature 497, 348 (2013).
  * (28) P. Bhattacharya, B. Xiao, A. Das, S. Bhowmick, and J. Heo, Solid state electrically injected exciton-polariton laser, Phys. Rev. Lett. 110, 206403 (2013).
  * (29) P. Bhattacharya, T. Frost, S. Deshpande, M. Z. Baten, A. Hazari, and A. Das, Room temperature electrically injected polariton laser, Phys. Rev. Lett. 112, 236802 (2014).
  * (30) C. Ciuti, Branch-entangled polariton pairs in planar microcavities and photonic wires, Phys. Rev. B 69, 245304 (2004).
  * (31) T. C. H. Liew and Y. G. Rubo, Quantum exciton-polariton networks through inverse four-wave mixing, Phys. Rev. B 97, 041302(R) (2018).
  * (32) D. Stefanatos and E. Paspalakis, Efficient entanglement generation between exciton-polaritons using shortcuts to adiabaticity, Opt. Lett. 43, 3313 (2018).
  * (33) J. Feng, H. Li, Z. Sun, and T. Byrnes, Entanglement generation and detection in split exciton-polariton condensates, Phys. Rev. A 108, 053301 (2023).
  * (34) Á. Cuevas et al., First observation of the quantized exciton-polariton field and effect of interactions on a single polariton, Sci. Adv. 4, eaao6814 (2018).
  * (35) T. Boulier, M. Bamba, A. Amo, C. Adrados, A. Lemaitre, E. Galopin, I. Sagnes, J. Bloch, C. Ciuti, E. Giacobino, and A. Bramati, Polariton-generated intensity squeezing in semiconductor micropillars, Nat. Commun. 5, 3260 (2014).
  * (36) P. V. Santos, and A. Fainstein, Polaromechanics: polaritonics meets optomechanics, Opt. Mater. Express 13, 1974 (2023).
  * (37) O. Kyriienko, T. C. H. Liew, and I. A. Shelykh, Optomechanics with cavity polaritons: Dissipative coupling and unconventional bistability, Phys. Rev. Lett. 112, 076402 (2014).
  * (38) B. Jusserand, A. N. Poddubny, A. V. Poshakinskiy, A. Fainstein, and A. Lemaitre, Polariton resonances for ultrastrong coupling cavity optomechanics in GaAs/AlAs multiple quantum wells, Phys. Rev. Lett. 115, 267402 (2015).
  * (39) D. L. Chafatinos, A. S. Kuznetsov, S. Anguiano, A. E. Bruchhausen, A. A. Reynoso, K. Biermann, P. V. Santos, and A. Fainstein, Polariton-driven phonon laser, Nat. Commun. 11, 4552 (2020).
  * (40) A. S. Kuznetsov, D. H. O. Machado, K. Biermann, and P. V. Santos, Electrically Driven Microcavity Exciton-Polariton Optomechanics at 20 GHz, Phys. Rev. X 11, 021020 (2021).
  * (41) N. Carlon Zambon, Z. Denis, R. De Oliveira, S. Ravets, C. Ciuti, I. Favero, and J. Bloch, Enhanced cavity optomechanics with quantum-well exciton polaritons, Phys. Rev. Lett. 129, 093603 (2022).
  * (42) A. S. Kuznetsov, K. Biermann, A. Reynoso, A. Fainstein, and P. V. Santos, Microcavity phonoritons—a coherent optical-to-microwave interface, Nat. Commun. 14, 5470 (2023).
  * (43) A. Fainstein, N. D. Lanzillotti-Kimura, B. Jusserand, and B. Perrin, Strong optical-mechanical coupling in a vertical GaAs/AlAs microcavity for subterahertz phonons and near-infrared light, Phys. Rev. Lett. 110, 037403 (2013).
  * (44) M. Aspelmeyer, T. J. Kippenberg, and F. Marquardt, Cavity optomechanics, Rev. Mod. Phys. 86, 1391 (2014).
  * (45) C. W. Gardiner and P. Zoller, Quantum Noise (Springer, Berlin, 2000).
  * (46) P. C. Parks and V. Hahn, Stability Theory (Prentice Hall, New York, 1993).
  * (47) D. Vitali, S. Gigan, A. Ferreira, H. R. Bohm, P. Tombesi, A. Guerreiro, V. Vedral, A. Zeilinger, and M. Aspelmeyer, Optomechanical entanglement between a movable mirror and a cavity field, Phys. Rev. Lett. 98, 030405 (2007).
  * (48) G. Vidal and R. F. Werner, Computable measure of entanglement, Phys. Rev. A 65, 032314 (2002); G. Adesso, A. Serafini, and F. Illuminati, Extremal entanglement and mixedness in continuous variable systems, Phys. Rev. A 70, 022318 (2004); M. B. Plenio, Logarithmic negativity: A full entanglement monotone that is not convex, Phys. Rev. Lett. 95, 090503 (2005).
  * (49) G. Adesso and F. Illuminati, Entanglement in continuous-variable systems: recent advances and current perspectives, J. Phys. A 40, 7821 (2007).
  * (50) C. Genes, A. Mari, P. Tombesi, and D. Vitali, Robust entanglement of a micromechanical resonator with output optical fields, Phys. Rev. A 78, 032316 (2008).
  * (51) J. Li, S.-Y. Zhu, and G. S. Agarwal, Magnon-photon-phonon entanglement in cavity magnomechanics, Phys. Rev. Lett. 121, 203601 (2018).
  * (52) X. Zuo, Z.-Y. Fan, H. Qian, M.-S. Ding, H. Tan, H. Xiong, and J. Li, Cavity magnomechanics: from classical to quantum, arXiv:2310.19237.
  * (53) Z.-Y. Fan, L. Qiu, S. Gröblacher, and J. Li, Microwave-optics Entanglement via Cavity Optomagnomechanics, Laser Photonics Rev. 2200866 (2023); arXiv:2208.10703.