Generating mechanical and optomechanical entanglement via pulsed interaction and measurement

For a pulsed optomechanical interaction [44], all the optical timescales involved are much less than the mechanical period, and the pulse is accommodated by operating in the unresolved-sideband regime, where the cavity decay rate κ far exceeds the mechanical angular frequency ω _m . We consider an optomechanical system comprising a cavity light mode and a mechanical mode interacting via radiation pressure. In a frame rotating at the cavity frequency, the optomechanical interaction is described by the Hamiltonian

$$\begin{equation}H=\hslash {\omega }_{\mathtt{\text{m}}}{b}^{{\dagger}}b-\hslash {g}_{0}{a}^{{\dagger}}a\left(b+{b}^{{\dagger}}\right),\end{equation} \tag{ 1 }$$

where g₀ is the intrinsic optomechanical coupling rate [66]. Here, a and b are the optical and mechanical annihilation operators, respectively, and we define the dimensionless optical quadratures as ${X}_{\mathtt{\text{l}}}=\left(a+{a}^{{\dagger}}\right)/\sqrt{2}$ and ${P}_{\mathtt{\text{l}}}=-\mathrm{i}\left(a-{a}^{{\dagger}}\right)/\sqrt{2}$, which satisfy the canonical commutation relation [X _l , P _l ] = i. The mechanical quadratures X _m and P _m are defined in the same way using the mechanical annihilation and creation operators.

The optomechanical Hamiltonian is linearized by transforming to a displaced frame such that a → α + a—where α is the mean intracavity amplitude of the pulse—and neglecting terms quadratic in a and a^† that describe the small quantum fluctuations. This linearization leads to

$$\begin{equation}{H}_{\mathtt{\text{lin}}}=\hslash {\omega }_{\mathtt{\text{m}}}{b}^{{\dagger}}b-\hslash {g}_{0}{\alpha }^{2}\left(b+{b}^{{\dagger}}\right)-\hslash {g}_{0}\alpha \left(a+{a}^{{\dagger}}\right)\left(b+{b}^{{\dagger}}\right),\end{equation} \tag{ 2 }$$

where we have taken $\alpha \in \mathbb{R}$ without loss of generality. In the unresolved-sideband regime, we describe a pulsed-optomechanical interaction with the unitary operator ${U}_{\mathtt{\text{lin}}}={\mathrm{e}}^{\mathrm{i}\lambda {X}_{\mathtt{\text{m}}}}\;{\mathrm{e}}^{\mathrm{i}\chi {X}_{\mathtt{\text{l}}}{X}_{\mathtt{\text{m}}}}$. The first exponential term in U _lin represents a deterministic momentum transfer to the mechanics due to the mean photon number, where λ ∝ α² g₀/κ quantifies this transfer in units of zero-point momentum. The second exponential in U _lin is a momentum transfer to the mechanical mode dependent upon the amplitude quadrature of the light, where the optomechanical interaction strength is given by χ ∝ αg₀/κ. The constants of proportionality in front of λ and χ are of order unity and depend on the optical pulse shape. From here on, we will use ${U}_{\mathtt{\text{om}}}={\mathrm{e}}^{\mathrm{i}\chi {X}_{\mathtt{\text{l}}}{X}_{\mathtt{\text{m}}}}$ to describe the optomechanical interaction, as the deterministic part ${\mathrm{e}}^{\mathrm{i}\lambda {X}_{\mathtt{\text{m}}}}$ may either be compensated for by applying a suitable displacement or subtracted from the measurement results in postprocessing.

As the linearized optomechanical Hamiltonian is bilinear in annihilation and creation operators, the unitary it generates will transform Gaussian states into other Gaussian states [67]. Such states can be fully characterized by the first and second moments of their quadrature vector $\mathbf{X}={\left({\mathbf{X}}_{\mathtt{\text{l}}},{\mathbf{X}}_{\mathtt{\text{m}}}\right)}^{\text{T}}={\left({X}_{\mathtt{\text{l}}},{P}_{\mathtt{\text{l}}},{X}_{\mathtt{\text{m}}},{P}_{\mathtt{\text{m}}}\right)}^{\text{T}}$. Here, the quadrature vectors for each mode are given by ${\mathbf{X}}_{\mathtt{\text{l}}}={\left({X}_{\mathtt{\text{l}}},{P}_{\mathtt{\text{l}}}\right)}^{\text{T}}$ and ${\mathbf{X}}_{\mathtt{\text{m}}}={\left({X}_{\mathtt{\text{m}}},{P}_{\mathtt{\text{m}}}\right)}^{\text{T}}$, respectively, and the canonical commutation relations can be written as [X, X^T] = iΩ, where the symplectic form Ω corresponding to our chosen ordering of quadrature operators, is in general, given by

$$\begin{equation}{\Omega}=\underset{i=1}{\overset{n}{\bigoplus}}\left(\begin{pmatrix}{cc}\hfill 0\hfill & \hfill 1\hfill \\ \hfill -1\hfill & \hfill 0\hfill \end{pmatrix}\right),\end{equation} \tag{ 3 }$$

where n = 2 for the bipartite optomechanical system. The first moments are given by the expectation value of the quadrature vector, ⟨X⟩, while the second moments correspond to elements of the covariance matrix σ, which are given by

$$\begin{equation}{\sigma }_{i,j}=\frac{1}{2}\langle {X}_{i}{X}_{j}+{X}_{j}{X}_{i}\rangle -\langle {X}_{i}\rangle \langle {X}_{j}\rangle .\end{equation} \tag{ 4 }$$

In this work, we utilize the block matrix form of the covariance matrix:

$$\begin{equation}\sigma =\left(\begin{pmatrix}{cc}A\hfill & \hfill C\hfill \\ \hfill {C}^{\text{T}}\hfill & \hfill B\hfill \end{pmatrix}\right).\end{equation} \tag{ 5 }$$

For a system comprising one optical and one mechanical mode, the A, B, and C blocks are 2 × 2 matrices, where A and B correspond to the optical and mechanical modes, respectively, and the off-diagonal block C corresponds to optomechanical correlations.

The set of Gaussian operations we employ in this work may be divided into unitary and non-unitary operations. A general Gaussian unitary on n modes can be written with a symplectic matrix $S\in Sp\left(2n,\mathbb{R}\right)$ acting on the quadrature vector X, where the real symplectic group $Sp\left(2n,\mathbb{R}\right)$ is defined as the group of 2n × 2n real symmetric matrices that preserve the commutation relations between the quadrature operators, namely SΩS^T = Ω. If the action of some Gaussian unitary U _g on a state ρ in Hilbert space is given by ${U}_{\text{G}}\rho {U}_{\text{G}}^{{\dagger}}$, then in quantum phase space the mapping is described by X → S X, where S is the symplectic matrix that has the same action on the quadrature vector as does the Gaussian unitary in the Heisenberg picture. For later convenience, we use the circle notation for quantum operations on density matrices, i.e. M◦ρ = MρM^†.

On the other hand, non-unitary Gaussian operations—optical loss, mechanical decoherence and optical homodyne measurements—are described differently. We model the decoherence of each mode using a phase-insensitive Gaussian channel, which corresponds to each mode interacting with a Gaussian environment on a beamsplitter [68, 69]. In Hilbert space we represent the action of these phase-insensitive Gaussian decoherence channels on state ρ as ${\mathcal{E}}_{\eta }\left(\rho \right)$ and ${\mathcal{E}}_{\gamma }\left(\rho \right)$, where η and γ are the efficiency of the optical channel and the damping rate of the mechanical mode, respectively. As we work in the pulsed regime of optomechanics, the mechanical decoherence may be neglected over the course of the optomechanical interactions and optical measurements.

Gaussian measurements may be described using the Krauss operator representation, i.e. ρ ∝ K _g ◦ρ, where K _g is the Krauss operator corresponding to the Gaussian measurement. Also, by considering the first-moments vector ⟨X⟩ and the covariance matrix of equation (5), a Gaussian measurement on the optical mode can be described through the mapping [69]

$$\begin{equation}\langle {\mathbf{X}}_{\mathtt{\text{m}}}\rangle \to \langle {\mathbf{X}}_{\mathtt{\text{m}}}\rangle +{C}^{\text{T}}{\left(A+{\sigma }_{\text{meas}}\right)}^{-1}\left(\langle {\mathbf{X}}_{\text{meas}}\rangle -\langle {\mathbf{X}}_{\mathtt{\text{l}}}\rangle \right),\end{equation} \tag{ 6 }$$

$$\begin{equation}B\to B-{C}^{\text{T}}{\left(A+{\sigma }_{\text{meas}}\right)}^{-1}C,\end{equation} \tag{ 7 }$$

where ⟨X_meas⟩ corresponds to the measurement outcome and σ_meas describes the Gaussian measurement [70].

Throughout this paper we describe quantum states as either density matrices in Hilbert space or as covariance matrices in phase space. We use these descriptions interchangeably and the corresponding symplectic transformations, together with unitary and non-unitary operations, can be found in appendix A.

The symplectic formalism outlined above may then be used to efficiently compute the generation of optical–mechanical entanglement and also be used to describe the effects of optical-quadrature measurements. Moreover, following a pulsed optomechanical interaction with an optical measurement allows one to reduce thermal noise along an axis of mechanical phase space, which we refer to as mechanical precooling.

The generation of an entangled optomechanical state is achieved via interaction of light and mechanics through U _om . However, the optical mode will then be subject to optical loss described by the channel ${\mathcal{E}}_{\eta }$, and we assume that the optical environment is well described by the vacuum state. After this optical loss channel, the mechanical state evolves freely through a phase-space angle θ, described by the unitary ${U}_{\mathtt{\text{rot}}}\left(\theta \right)={\mathrm{e}}^{\mathrm{i}\theta {b}^{{\dagger}}b}$. During this time, the mechanical mode interacts with a thermal environment with mean occupation $\overline{N}$ at a rate γ, described by ${\mathcal{E}}_{\gamma }$. Including this decoherence, the map that describes the generation of an optical–mechanical entangled state from initial state ρ is therefore

$$\begin{equation}\rho \to {\mathcal{E}}_{\gamma }\left({U}_{\mathtt{\text{rot}}}\left(\theta \right){\circ}{\mathcal{E}}_{\eta }\left({U}_{\mathtt{\text{om}}}{\circ}\rho \right)\right).\end{equation} \tag{ 8 }$$

Note that the phase insensitive channel commutes with the free evolution of the mechanical mode, i.e. ${\mathcal{E}}_{\gamma }\left({U}_{\mathtt{\text{rot}}}{\circ}\rho \right)={U}_{\mathtt{\text{rot}}}{\circ}{\mathcal{E}}_{\gamma }\left(\rho \right)$.

The mechanical oscillator may be cooled prior to the generation of optical–mechanical entanglement by homodyning the output light following an optomechanical interaction. More specifically, this precooling stage is described by

$$\begin{equation}\rho \to {\mathcal{E}}_{\gamma }\left({U}_{\mathtt{\text{rot}}}\left(\pi /2\right)\vert {P}_{\mathtt{\text{l}}}\rangle \langle {P}_{\mathtt{\text{l}}}\vert {\circ}{\mathcal{E}}_{\eta }\left({U}_{\mathtt{\text{om}}}{\circ}\rho \right)\right).\end{equation} \tag{ 9 }$$

Mechanical precooling, equation (9), begins with a pulsed optomechanical interaction, followed by optical decoherence. Subsequently, a phase-quadrature measurement, with projector |P _l ⟩⟨P _l |, is made on the light which projects the mechanical mode into a squeezed thermal state with a lower effective thermal occupation [44]. The mechanical state is then allowed to evolve freely over a quarter of a period prior to the generation of optical–mechanical entanglement, in which time the mechanical mode undergoes some decoherence. In this way, equation (9) prepares a mechanical state with reduced thermal noise along the momentum quadrature of phase space. The entanglement and precooling stages are illustrated in the circuit diagram of figure 1(a).

Zoom In Zoom Out Reset image size

Our proposed experimental setup for the preparation of an entangled state of light and mechanics is shown in figure 1(b). A pulse of light in a coherent state |α⟩ is injected into an optomechanical cavity, which then interacts via the entangling unitary U _om with a mechanical oscillator in an initial thermal state, with mean thermal occupation $\overline{n}$. This interaction, followed by optical loss and mechanical decoherence is described by equation (8) and generates an entangled state of light and mechanics.

Prior to the entangling stage, the precooling stage, described by equation (9), creates an initial mechanical state which can generate more entanglement for a given interaction strength χ. This increase in performance is due to the reduction in thermal noise in the direction of the phase-space translations generated by the optomechanical unitary, and hence a higher amount of entanglement can be prepared between light and mechanics. Multiple precooling stages can be employed to further increase the degree of entanglement generated. However, within the range of parameters we consider in this work, we find that additional stages are unnecessary and lead to a negligible increase in the amount of entanglement generated. The analytic results of the precooling and entangling stages are given in appendix A.

Once an entangled state of light and mechanics has been generated, the full phase-space structure of the two-mode Gaussian state can be reconstructed. In this reconstruction process many runs of an experiment will be needed to build up adequate statistics for the different quadrature combinations. Experimentally, it will be important to track the first moments vector to accurately construct probability histograms of each observed quadrature. More specifically, in order to reconstruct the statistics of each quadrature the first moment must be recorded and subtracted from the measurement results, removing the effects of the random homodyne measurement outcomes from the data. Tracking the first moments vector may be achieved by considering the action of each Gaussian unitary and non-unitary process in the precooling and entangling stages of the protocol, as well as recording the measurement outcomes of the homodyne detection events. As mentioned in section 2, when constructing the probability histograms of the mechanical quadratures, the projection of the deterministic momentum transfer along the mechanical quadrature under consideration may also be subtracted from the measurement results to avoid the necessity of a feedback force.

Turning our attention now to the reconstruction of the covariances, we refer to the reconstructed covariance matrix as σ _ver , which in block matrix form is

$$\begin{equation}{\sigma }_{\mathtt{\text{ver}}}=\left(\begin{pmatrix}{cc}\hfill {A}_{\mathtt{\text{ver}}}\hfill & \hfill {C}_{\mathtt{\text{ver}}}\hfill \\ \hfill {C}_{\mathtt{\text{ver}}}^{\text{T}}\hfill & \hfill {B}_{\mathtt{\text{ver}}}\hfill \end{pmatrix}\right),\end{equation} \tag{ 10 }$$

and will approach σ in the absence of decoherence processes. In our analysis, we assume that the interaction strength χ and the total optical efficiency η are precisely known, which allows σ _ver to be constructed accurately. So that the reader may easily follow the symplectic transformations, we present formulas for σ _ver without the inclusion of optical loss in only the verification stage—optical losses are present in the precooling and entangling stages. However, including optical loss in the formulas for σ _ver is straightforward and does not affect the results we present. In appendix B, we discuss how the elements of σ _ver are obtained in more detail.

After the optomechanical interaction in the entangling stage of the protocol, the optical pulse exits the cavity and either a homodyne measurement is made on a rotated quadrature P _l (φ) = P _l  cos φ − X _l  sin φ, or the light is sent down a delay line. If the entangling pulse of light is directly homodyned, this allows the A _ver block of σ _ver to be determined and the required set of homodyne measurement angles may be noted from the expression

$$\begin{equation}{A}_{\mathtt{\text{ver}}}=\left(\begin{pmatrix}{cc}\mathrm{Var}\left({X}_{\mathtt{\text{l}}}\right)\hfill & \hfill \frac{1}{2}\left[\mathrm{Var}\left({P}_{\mathtt{\text{l}}}\left(\frac{3\pi }{4}\right)\right)-\mathrm{Var}\left({P}_{\mathtt{\text{l}}}\left(\frac{\pi }{4}\right)\right)\right]\hfill \\ \hfill \frac{1}{2}\left[\mathrm{Var}\left({P}_{\mathtt{\text{l}}}\left(\frac{3\pi }{4}\right)\right)-\mathrm{Var}\left({P}_{\mathtt{\text{l}}}\left(\frac{\pi }{4}\right)\right)\right]\hfill & \hfill \mathrm{Var}\left({P}_{\mathtt{\text{l}}}\right)\hfill \end{pmatrix}\right).\end{equation} \tag{ 11 }$$

After the entangling pulse has been homodyned, a verification pulse in a pure coherent state is then introduced into the cavity after a controllable time delay t = θ/ω _m in which time the mechanics evolves freely and undergoes some decoherence. This decoherence of the mechanical mode between the interaction with the entangling and verification pulses leads to imperfect reconstruction of the covariance matrix σ. After the optomechanical interaction, the quadrature vector of the verification pulse is ${\mathbf{X}}_{\mathtt{\text{v}}}\left(\theta \right)={\left({X}_{\mathtt{\text{v}}}\left(\theta \right),{P}_{\mathtt{\text{v}}}\left(\theta \right)\right)}^{\text{T}}$. The X _v (θ) quadrature contains no information about the mechanical mode, however the phase quadrature of the verification pulse reads P _v (θ) = P _in + χX _m (θ), where the initial optical phase quadrature P _in introduces noise in the reconstruction process and X _m (θ) = X _m  cos θ + P _m  sin θ. Note that the argument θ of P _v (θ) refers to the phase-space evolution of the mechanics and not a rotation in the optical subspace. To access the B _ver block of the covariance matrix, optical homodyne measurements are made directly on the P _v (θ) quadrature at various time delays and the input noise Var(P _in ) can be subtracted to increase the accuracy of the reconstruction process. Importantly, to verify the B _ver block the homodyne measurement outcomes on the entangling pulse of light must be ignored so the light is traced out in the reconstruction process. The mechanical block is then given by

$$\begin{equation}{B}_{\mathtt{\text{ver}}}=\frac{1}{{\chi }^{2}}\left(\begin{pmatrix}{cc}\mathrm{Var}\left({P}_{\mathtt{\text{v}}}\left(0\right)\right)-\mathrm{Var}\left({P}_{\mathtt{\text{in}}}\right)\hfill & \hfill \frac{1}{2}\left[\mathrm{Var}\left({P}_{\mathtt{\text{v}}}\left(\frac{\pi }{4}\right)\right)-\mathrm{Var}\left({P}_{\mathtt{\text{v}}}\left(\frac{3\pi }{4}\right)\right)\right]\hfill \\ \hfill \frac{1}{2}\left[\mathrm{Var}\left({P}_{\mathtt{\text{v}}}\left(\frac{\pi }{4}\right)\right)-\mathrm{Var}\left({P}_{\mathtt{\text{v}}}\left(\frac{3\pi }{4}\right)\right)\right]\hfill & \hfill \mathrm{Var}\left({P}_{\mathtt{\text{v}}}\left(\frac{\pi }{2}\right)\right)-\mathrm{Var}\left({P}_{\mathtt{\text{in}}}\right)\hfill \end{pmatrix}\right).\end{equation} \tag{ 12 }$$

To access the C _ver -block elements of the covariance matrix, the entangling light is not directly homodyned, but instead it is sent down a path towards a controllable delay line, while the verification pulse is sent down another path. These two pulses interact with each other on a 50:50 beamsplitter and each individual output is directly homodyned. Redirection of the entangling and verification pulses is achieved by using an optical switch as shown in figure 1(b). The timescale over which this optical switch operates must be less than a quarter of the mechanical cycle as this corresponds to the time between the precooling pulse and the entangling pulse of light. This means there is no need for high-speed switching or small-scale integrated operation, one can use a larger switch that operates with very low optical loss, such as a Pockels cell.

The beamsplitter operation, which mixes the entangling and verification pulses, is described by the 4 × 4 symplectic matrix S _bs (α, β) given in appendix A, where α and β are the real beamsplitter parameters. To access every element of C _ver , two different beamsplitter configurations are required: ${S}_{\mathtt{\text{bs}}}\left(\frac{\pi }{4},0\right)$ and ${S}_{\mathtt{\text{bs}}}\left(\frac{\pi }{4},\frac{\pi }{2}\right)$. These two configurations may be implemented with a single 50:50 beamsplitter and a controllable phase shifter. We define the quadrature vectors after the ${S}_{\mathtt{\text{bs}}}\left(\frac{\pi }{4},0\right)$ and ${S}_{\mathtt{\text{bs}}}\left(\frac{\pi }{4},\frac{\pi }{2}\right)$ operations as ${\left({\mathbf{X}}_{\mathtt{\text{v}}}^{\prime }\left(\theta \right),{\mathbf{X}}_{\mathtt{\text{l}}}^{\prime }\left(\theta \right)\right)}^{\text{T}}$ and ${\left({\mathbf{X}}_{\mathtt{\text{v}}}^{{\prime\prime}}\left(\theta \right),{\mathbf{X}}_{\mathtt{\text{l}}}^{{\prime\prime}}\left(\theta \right)\right)}^{\text{T}}$, respectively. Then the C _ver block are determined by

$$\begin{equation}{C}_{\mathtt{\text{ver}}}=\frac{1}{2\chi }\left(\begin{pmatrix}{cc}\hfill \mathrm{Var}\left({P}_{\mathtt{\text{v}}}^{{\prime\prime}}\left(2\pi \right)\right)-\mathrm{Var}\left({X}_{\mathtt{\text{l}}}^{{\prime\prime}}\left(2\pi \right)\right)\hfill & \hfill \mathrm{Var}\left({P}_{\mathtt{\text{v}}}^{{\prime\prime}}\left(\frac{5\pi }{2}\right)\right)-\mathrm{Var}\left({X}_{\mathtt{\text{l}}}^{{\prime\prime}}\left(\frac{5\pi }{2}\right)\right)\hfill \\ \hfill \mathrm{Var}\left({P}_{\mathtt{\text{v}}}^{\prime }\left(2\pi \right)\right)-\mathrm{Var}\left({P}_{\mathtt{\text{l}}}^{\prime }\left(2\pi \right)\right)\hfill & \hfill \mathrm{Var}\left({P}_{\mathtt{\text{v}}}^{\prime }\left(\frac{5\pi }{2}\right)\right)-\mathrm{Var}\left({P}_{\mathtt{\text{l}}}^{\prime }\left(\frac{5\pi }{2}\right)\right)\hfill \end{pmatrix}\right).\end{equation} \tag{ 13 }$$

Interference between the entangling and verification pulses is necessary to obtain the optomechanical covariances. Whereas, if there was no interference one would only be able to access the A _ver and B _ver blocks and not the full covariance matrix. Note, that we have included an addition of 2π in the arguments of the C _ver elements, which corresponds to a full rotation of the mechanical mode in phase space. The necessity of this additional 2π is explained in section 3.2. Furthermore, when the additional factor of 2π is included, the time between the entangling and verification pulses is always greater than the time between the precooling and entangling pulses.

In order to quantify entanglement, we employ the logarithmic negativity, which is a well-suited entanglement monotone for low-dimensional discrete variable systems [71–74] and for bipartite Gaussian states [75–77]. This monotone quantifies the extent to which the positivity of the density operator is violated after a transpose operation is applied to a single mode. For bipartite Gaussian systems, partial transposition leads to a covariance matrix with eigenvalues given by

$$\begin{equation}{\tilde {\nu }}_{{\pm}}=\sqrt{\frac{\tilde {{\Delta}}{\pm}\sqrt{{\tilde {{\Delta}}}^{2}-4\;\mathrm{det}\;\sigma }}{2}}\end{equation} \tag{ 14 }$$

where $\tilde {{\Delta}}=\mathrm{det}\;A+\mathrm{det}\;B-2\;\mathrm{det}\;C$. A necessary and sufficient condition for bipartite entanglement in Gaussian continuous variable systems is given by ${\tilde {\nu }}_{-}{< }1/2$. The logarithmic negativity of the two-mode state with covariance matrix σ is then calculated as

$$\begin{equation}{E}_{\mathcal{N}}\left(\sigma \right)=\mathrm{max}\left\{0,-{\mathrm{log}}_{2}\left(2{\tilde {\nu }}_{-}\right)\right\}.\end{equation} \tag{ 15 }$$

Having now outlined the protocol for optical–mechanical entanglement preparation and verification, we now consider two additional subtleties the protocol presents. Firstly, we investigate the non-monotonic behaviour in the logarithmic negativity with increasing interaction strength in the presence of both optical loss and mechanical decoherence, and how the protocol can be made more resilient to these decoherence effects by using optical squeezers. Secondly, we consider the time order in which elements of the covariance matrix must be measured to conservatively estimate the entanglement generated in the protocol. To facilitate this quantitative discussion, in table 1 we list values for our system parameters based on the sliced-silicon nanobeam architecture comprising the optomechanical system that appears in references [51–53]. In the following, we also discuss the parameters we change from these experiments, namely the thermal occupation $\overline{n}$ and the optical efficiency η.

Table 1.  List of proposed experimental parameters for the optical–mechanical and mechanical entanglement protocols. The parameter values are based on the sliced-silicon nanobeam system that appears in references [51–53]. Here, we assume that the mechanical oscillator and its environment are in equilibrium, hence $\overline{n}=\overline{N}=500$—this corresponds to a temperature of approximately 0.1 K. We propose an improved optical efficiency of η = 0.855. This optical efficiency may be decomposed into a part that accounts for output cavity coupling losses from the cavity, η _cav , and a part that accounts for homodyne detection efficiency η _det . In this work, we consider a range of values for the optomechanical interaction strength χ.

System parameter	Value
ω m /2π	4 MHz
γ/2π	100 Hz
$\overline{n}=\overline{N}$	500
κ/2π	20 GHz
g0/2π	30 MHz
η cav	0.9
η det	0.95

In reference [51], a reduction in the thermal occupation from $\overline{n}=22\;200$ to an effective occupation of 3 400 was achieved using a series of pulsed-optomechanical interactions at 3.2 K. Therefore, multiple precooling steps can be used to increase the temperature at which the protocol can operate at. However, to simplify the discussions, we do not consider multiple precooling steps and we assume a bath temperature of approximately 0.1 K, hence $\overline{n}=\overline{N}=500$. Such temperatures are readily achieved with commercially available dilution refrigerators. Furthermore, in the table we propose an improved optical efficiency of η = 0.855. Optical efficiencies of 1.3% have been achieved in experiments, and without implementing further cooling techniques a reduction in the uncertainty below the zero-point fluctuations may be achieved with η > 8% [51, 53]. In section 5, we further discuss the requirements on optical efficiency that our protocols present. In particular, the optical–mechanical entanglement protocol can generate entanglement at very low optical efficiencies, less than 1%, while the protocols for generating mechanical entanglement require total optical efficiencies greater than 10%. As noted in reference [53], improvements to the optical design and employing new waveguide coupling techniques provide an encouraging route for further experimental progress in pulsed optomechanics.

After the optomechanical interaction between the entangling pulse and the precooled mechanical oscillator, the reduced state of the optical mode will satisfy Var(P _l ) > Var(X _l ). This is due to the phase-space displacements generated by U _om along the P _l axis of optical phase space, which are proportional to the mechanical position quadrature operator. The optical mode is then subject to the phase-insensitive decoherence channel ${\mathcal{E}}_{\eta }$, which introduces vacuum noise and degrades the optomechanical correlations. However, to reduce these deleterious effects, an optical squeezer may be applied immediately after the optomechanical interaction. Optical squeezing is a useful way to protect single-mode quantum states from interactions with the environment, for example squeezers may be used to improve state fidelity in teleportation [78] and the performance of quantum memory devices [79]. Furthermore, the utility of such operations in protecting quantum coherence has been investigated theoretically [80–82] and verified experimentally [83]. Unlike these recent applications of optical squeezing, in this work we utilize optical squeezing to protect the entanglement of a bipartite state specifically. Equation (8), which describes the generation of optical–mechanical entanglement, is therefore modified to

$$\begin{equation}\rho \to {\mathcal{E}}_{\gamma }\left({U}_{\mathtt{\text{rot}}}\left(\theta \right){\circ}{\mathcal{E}}_{{\eta }_{\mathtt{\text{det}}}}\left({U}_{\mathtt{\text{sq}}}{\circ}{\mathcal{E}}_{{\eta }_{\mathtt{\text{cav}}}}\left({U}_{\mathtt{\text{om}}}{\circ}\rho \right)\right)\right).\end{equation} \tag{ 16 }$$

Here, the squeezing unitary is ${U}_{\mathtt{\text{sq}}}=\mathrm{exp}\left\{\frac{r}{2}\left[{\left({a}^{{\dagger}}\right)}^{2}-{a}^{2}\right]\right\}$, with real squeezing parameter r. Moreover, the channel ${\mathcal{E}}_{{\eta }_{\mathtt{\text{cav}}}}$ represents losses from the optical cavity and ${\mathcal{E}}_{{\eta }_{\mathtt{\text{det}}}}$ is the channel describing detection inefficiencies. We assume that the optical squeezer can be turned off during the precooling and verification stages of the protocol. However, if for a particular experimental implementation of our scheme, operating the squeezer in this way is not feasible, then the description of the precooling and verification stages is straightforwardly modified by using the symplectic transformation S _sq (r), which appears in appendix A.

In our entanglement scheme, we are interested in how squeezing operations allow for the protection of entanglement—rather than quantum superposition as studied in references [80–83]. For a Gaussian bosonic channel, it is well known that the output state with maximum purity and minimum von-Neumann entropy corresponds to a coherent state input [84, 85]. Likewise, when a squeezing operation is applied to an input state of fixed given purity, the input state that maximizes the output purity and minimizes the output entropy, with respect to the real squeezing parameter r, is given by a symmetric mode in phase space, i.e. a thermal state, see appendix C. This may suggest that the optical squeezing parameter ${r}_{\mathtt{\text{sym}}}=\frac{1}{4}\;\mathrm{ln}\left[1-{\eta }_{\mathtt{\text{cav}}}+{\eta }_{\mathtt{\text{cav}}}\left(1+2{V}_{\text{x}}{\chi }^{2}\right)\right]$, which symmetrizes the optical mode such that Var(P _l ) = Var(X _l ), is a good candidate for the optimal squeezing parameter to also maximize the output logarithmic negativity. Here, V_x is the position variance of the precooled mechanical state, which is derived in appendix A.

However, the entangling unitary U _om generates optomechanical correlations which are distributed asymmetrically between the optical phase and amplitude quadratures. Namely, only the optical phase quadrature contains information about the mechanical position, which may be lost to the optical environment leading to a reduction in the logarithmic negativity. Furthermore, the optical environment will be less sensitive to displacements along the phase quadrature for an optical state of a given purity with Var(P _l ) > Var(X _l ). Therefore, a state with Var(P _l ) > Var(X _l ) at the input of the optical loss channel ${\mathcal{E}}_{{\eta }_{\mathtt{\text{det}}}}$, reduces the sensitivity of the optical environment to the mechanical position quadrature and the optomechanical correlations.

A balance between a symmetric optical phase-space distribution, which minimises the output entropy of the optical decoherence channel, and an optical phase-space distribution with Var(P _l ) > Var(X _l ), which reduces the sensitivity of the optical environment to optomechanical correlations, leads to the optimal squeezing parameter r _opt that maximizes the logarithmic negativity. In the upper plot of figure 2, r _opt is plotted as a function of the interaction strength χ. This plot shows the logarithmic negativity of the optical–mechanical entangled state as a function of r and χ. Here, the entangled state has been subject to both optical and mechanical decoherence processes since the time of entanglement generation. For comparison, we also present r _sym as a function of χ, and we note that r _sym > r _opt for all χ—demonstrating the balance between the two effects described above. The experimental implementation of this squeezing operation, which we introduce to protect entanglement in pulsed optomechanics, is made feasible by experimental advances in squeezing pulses of light [86, 87] and non-classical continuous travelling waves of light [88].

Zoom In Zoom Out Reset image size

Importantly, as is evident from the lower plot of figure 2, non-monotonic behaviour in the logarithmic negativity is observed when both optical loss and mechanical decoherence are present after the generation of entanglement. The existence of this non-monotonic behaviour is because—at large values of χ—the increase in the magnitude of the optomechanical correlations with interaction strength scales in the same way as the decrease in the logarithmic negativity due to the decoherence of either the optical or mechanical mode alone. Hence, when both optical and mechanical decoherence processes are present, the logarithmic negativity first rises to a maximum before tending towards zero as χ is increased and decoherence processes outweigh the effect of the entangling optomechanical unitary. When the optical squeezing operation can be applied before the decoherence channels, the non-monotonic behaviour in logarithmic negativity can be completely eliminated. In this case, the entangled state is made more resilient to loss of information to the environment before it enters the decoherence channels and so higher values of ${E}_{\mathcal{N}}\left(\sigma \right)$ can be reached than in figure 2. We refer the reader to appendix D for a discussion of this case.

The lower plot in figure 2 also demonstrates how accurately the optical–mechanical entanglement is measured in the verification process. As it is necessary to allow the mechanical oscillator to undergo free mechanical evolution to obtain information about the rotated mechanical quadratures X _m (θ), we encounter the problem that different elements of the covariance matrix must be measured at different times, having suffered different amounts of mechanical decoherence. Moreover, care must be taken when calculating the logarithmic negativity from the reconstructed covariance matrix, in particular one must ensure that ${E}_{\mathcal{N}}\left(\sigma \right){\geqslant}{E}_{\mathcal{N}}\left({\sigma }_{\mathtt{\text{ver}}}\right)$, meaning that the logarithmic negativity is not overestimated. To ensure this condition, we demand that the elements of C _ver are accessed at later times than those of A _ver and B _ver , hence the addition of 2π in the argument of the elements in equation (13). Therefore, the C _ver elements experience a larger amount of decoherence than the A _ver and B _ver elements, meaning that the optical–mechanical correlations cannot be falsely enhanced in the verification process. We refer to this approach as conservative in time, and the conservative nature of this approach is demonstrated in the lower plot of figure 2. Moreover, we consider the case in which state reconstruction is performed without assuming that the noise statistics of the verification pulse Var(P _in ) is known and can be subtracted. Similarly, we refer to this approach as conservative in time and noise.

If the mechanical damping rate γ and thermal occupation $\overline{n}$ are well known, the inverse of the mechanical decoherence may be applied to each element of σ _ver to compensate for decoherence, as illustrated in the lower plot of figure 2 by the vertical arrow. The application of the inverse map to the elements of σ _ver that depend on two different times is more involved and is described in appendix B. Also, if decoherence can be accurately compensated for in this way, the need to be conservative with respect to time is no longer needed and so additional factors of 2π can be removed from the arguments of the C _ver -elements.

The interferometric protocol to entangle two mechanical modes is shown in the circuit diagram of figure 3, where, for brevity, the details of the precooling stage have been omitted. The entanglement stage of the protocol starts with a pulse of light in a coherent state impinging on one input port of a Mach–Zehnder interferometer. The other input port is in the vacuum state, such that after the first 50:50 beamsplitter interaction, described by U _bs , two coherent pulses of light propagate in the upper and lower arms towards the two mechanical modes. We label the quadratures of the optical modes ${\mathbf{X}}_{\mathtt{\text{l}}1}={\left({X}_{\mathtt{\text{l}}1},{P}_{\mathtt{\text{l}}1}\right)}^{\text{T}}$ and ${\mathbf{X}}_{\mathtt{\text{l}}2}={\left({X}_{\mathtt{\text{l}}2},{P}_{\mathtt{\text{l}}2}\right)}^{\text{T}}$, and the mechanical modes with which these light modes interact as X _m1 and X _m2, respectively. Each light mode interacts with its respective mechanical mode through the unitary U _om , hence the total optomechanical interaction is described by ${U}_{\mathtt{\text{om}}}\otimes {U}_{\mathtt{\text{om}}}={U}_{\mathtt{\text{om}}}^{\otimes 2}$. After these interactions, the decoherence channel ${\mathcal{E}}_{{\eta }_{\mathtt{\text{cav}}}}$ describing cavity losses acts on the two-mode optical subspace, and optional unitary squeezing operations ${U}_{\mathtt{\text{sq}}}^{\otimes 2}$ may be applied to each optical mode to increase the final amount of entanglement generated. Following this, the optical modes experience the loss channel ${\mathcal{E}}_{{\eta }_{\mathtt{\text{det}}}}$ and then pass through a 50:50 beamsplitter, which erases which-path information about the optomechanical interactions. In our model, the ordering of the final loss channel and beamsplitter is arbitrary as they both lead to the same final mechanical state. To turn the optical–mechanical entanglement to mechanical entanglement, optical homodyne measurements are made on the rotated quadratures X _l1(φ) and X _l2(ψ), described by the projector |X _l1(φ), X _l2(ψ)⟩⟨X _l1(φ), X _l2(ψ)|. After the mechanical entanglement has been generated, the first and second mechanical modes evolve through phase-space angles θ and ϕ, respectively, which is described by U _rot (θ, ϕ) = U _rot (θ) ⊗ U _rot (ϕ). Different periods of mechanical free evolution can be achieved using controllable delay lines, as is shown in figure 3(b) and discussed in further detail below. During this free evolution the mechanical oscillators decohere via the channel ${\mathcal{E}}_{\gamma }$. Including mechanical decoherence, the interferometric entanglement protocol is therefore given by

$$\begin{equation}\rho \to {\mathcal{E}}_{\gamma }\left({U}_{\mathtt{\text{rot}}}\left(\theta ,\phi \right)\vert {X}_{\mathtt{\text{l}}1}\left(\varphi \right),{X}_{\mathtt{\text{l}}2}\left(\psi \right)\rangle \langle {X}_{\mathtt{\text{l}}1}\left(\varphi \right),{X}_{\mathtt{\text{l}}2}\left(\psi \right)\vert {U}_{\mathtt{\text{bs}}}{\circ}\right.\left.{\mathcal{E}}_{{\eta }_{\mathtt{\text{det}}}}\left({U}_{\mathtt{\text{sq}}}^{\otimes 2}{\circ}{\mathcal{E}}_{{\eta }_{\mathtt{\text{cav}}}}\left({U}_{\mathtt{\text{om}}}^{\otimes 2}{U}_{\mathtt{\text{bs}}}{\circ}\rho \right)\right)\right).\end{equation} \tag{ 17 }$$

Zoom In Zoom Out Reset image size

From equation (17), one finds that at the input ports of the homodyne detectors the optical phase quadratures P _l1 and P _l2 contain information about the mechanical EPR quadratures X _m1 ± X _m2, whilst no knowledge of mechanical position can be gleaned from the the corresponding optical amplitude quadratures, X _l1 and X _l2. In what follows, we therefore choose that the homodyne angles are set to one of two equivalent configurations given by $\left(\varphi ,\psi \right)=\left\{\left(0,\pi /2\right),\left(\pi /2,0\right)\right\}$. These configurations correspond to Bayesian inference of a mechanical EPR quadrature, which projects the mechanics into an entangled state. Conversely, if both phase quadratures are measured, then no mechanical entanglement can be created, as this would allow information about each mechanical quadrature to be obtained separately. Interestingly, there exists other optimal homodyne angles than those stated above, which also maximize the logarithmic negativity. These angles correspond to a balance between the Bayesian inference of an EPR quadrature and an effective unitary operation, and this is discussed further in appendix E.

The state prepared by the precooling and entangling stages, as illustrated in figure 3, can be reconstructed by sending a subsequent verification pulse into the input port of the interferometer. The verification pulses in the upper and lower arms then interact with the mechanical oscillators after time delays t₁ = θ/ω _m and t₂ = ϕ/ω _m , respectively. The phase quadrature of the verification pulse in the upper arm of the Mach–Zehnder interferometer evolves from P _in1 to P _v1(θ) = P _in1 + χX _m1(θ), during the optomechanical interaction. While the verification pulse in the lower arm evolves from P _in2 to P _v2(ϕ) = P _in2 + χX _m2(ϕ). As in the optical–mechanical entanglement scheme, for clarity, we present results for the verification procedure without the inclusion of optical loss and refer the reader to appendix B for a more complete discussion, including optical loss.

We represent the covariance matrix obtained in the verification procedure as σ _ver , which has the same block-matrix form as equation (10). However, now both the A _ver and B _ver blocks correspond to a covariance matrix of a reduced mechanical state. These matrix blocks are obtained by performing phase homodyne measurements on the verification pulses after they interact with the mechanical modes. To allow for this direct detection, optical switches are placed after the optomechanical interactions in figure 3(b) to redirect the verification pulses away from the second 50:50 beamsplitter of the Mach–Zehnder interferometer. The elements of the A _ver and B _ver blocks, which depend on the mechanical phase-space angles θ and ϕ, respectively, are therefore given by equivalent expressions to equation (12)—with P _v (θ) replaced by P _v1(θ) or P _v2(ϕ).

On the other hand, the correlations between the two mechanical modes, defined by the C _ver -block elements, are obtained by mixing the two verification pulses on a beamsplitter after the optomechanical interactions. Similarly to the verification of optical–mechanical entanglement, it is the interference of these two optical pulses that allows for the full reconstruction of the covariance matrix. After this operation, the resulting phase quadrature operators at the beamsplitter outputs are given by ${P}_{\mathtt{\text{v}}{\pm}}\left(\theta ,\phi \right)=\left({P}_{\mathtt{\text{v}}1}\left(\theta \right){\pm}{P}_{\mathtt{\text{v}}2}\left(\phi \right)\right)/\sqrt{2}$. Homodyne measurements are then performed on these quadratures in order to construct C _ver . Introducing δ(θ, ϕ) = Var(P _v+(θ, ϕ)) − Var(P _v−(θ, ϕ)) allows one to write C _ver as

$$\begin{equation}{C}_{\mathtt{\text{ver}}}=\frac{1}{2{\chi }^{2}}\left(\begin{pmatrix}{cc}\delta \left(2\pi ,2\pi \right)\hfill & \hfill \delta \left(2\pi ,\frac{5\pi }{2}\right)\hfill \\ \hfill \delta \left(\frac{5\pi }{2},2\pi \right)\hfill & \hfill \delta \left(\frac{5\pi }{2},\frac{5\pi }{2}\right)\hfill \end{pmatrix}\right).\end{equation} \tag{ 18 }$$

The off-diagonal elements of C _ver , correspond to different periods of free evolution for each mechanical mode. To introduce this difference—and to ensure the recombination of the verification pulses on the final beamsplitter—controllable delay lines may be inserted before and after the optomechanical interactions in the interferometer of figure 3(b). We have again included an additional factor of 2π in the arguments of the C _ver -block such that the logarithmic negativity is conservatively estimated as discussed in section 3.2.

Figure 4 shows the logarithmic negativity of the entangled state of two mechanical oscillators, and demonstrates how optical squeezing may be used to increase the amount of entanglement generated. The squeezing changes the phase-space distributions of the light modes to ensure the protocol remains robust to decoherence processes in the large χ limit. The reduced state of each mechanical mode is identical—see appendix A—meaning the Gaussian state is symmetric, which allows the entanglement of formation to be computed analytically [89, 90]. But as this entanglement measure is also a monotonically decreasing function of ${\tilde {\nu }}_{-}$, it is an equivalent measure to ${E}_{\mathcal{N}}$. The lower plot of figure 4 also reaffirms the conclusions of section 3.2. Namely, that the origin of the non-monotonic behaviour in the logarithmic negativity is due to the presence of both optical loss and mechanical decoherence. Conservative approaches for state verification are also shown in the lower plot of figure 4. As in the case of the optical–mechanical entanglement protocol, the statistics obtained via the conservative approaches may be compensated for by using an inverse decoherence map.

Zoom In Zoom Out Reset image size

The interferometric scheme we have introduced above offers an alternative strategy to the mechanical entanglement scheme presented in reference [33]—which utilizes Stokes scattering in the resolved-sideband regime. Namely, in this work, the potential to generate entanglement using short optical pulses in the unresolved sideband regime is explored, which allows rapid preparation and verification of entanglement, and we introduce and detail a full-state-characterization procedure using these short optical pulses.

The non-interferometric scheme for entanglement preparation and verification using pulsed optomechanics is shown in the circuit diagram of figure 5. As before, a precooling stage is implemented to increase the amount of mechanical entanglement generated. The entangling stage is then initiated by a pulse of coherent light entering an optomechanical cavity and interacting with the first mechanical mode—identified by its quadrature vector X _m1. The light is then subject to optical losses, each described by ${\mathcal{E}}_{{\eta }_{\mathtt{\text{cav}}}}$, corresponding to the output-cavity-coupling and input-cavity-coupling efficiencies of the first and second cavities, respectively. Here, we assume equal outcoupling and incoupling efficiencies, which are given by η _cav . An optional optical squeezer may be placed between these two optical loss channels to provide resilience to interactions with the optical environment.

Zoom In Zoom Out Reset image size

For the first optomechanical interaction, it was sufficient to absorb the input-cavity-coupling losses into the definition of the interaction strength χ, but when the light interacts with the second mechanical mode—described by X _m2—it is important to separate out this effect. This is because the light incident on the second cavity contains information about X _m1, and so input coupling losses lead to loss of information about the mechanical position and a reduction in the final amount of entanglement generated. The attenuation of the light between the two optomechanical interactions leads to a reduction in the optomechanical interaction strength: χ → η _cav χ. Hence the second optomechanical interaction is modified from ${U}_{\mathtt{\text{om}}}={\mathrm{e}}^{\mathrm{i}\chi {X}_{\mathtt{\text{l}}}{X}_{\mathtt{\text{m}}}}$ to ${U}_{\mathtt{\text{om}}}^{\prime }={\mathrm{e}}^{\mathrm{i}{\eta }_{\mathtt{\text{cav}}}\chi {X}_{\mathtt{\text{l}}}{X}_{\mathtt{\text{m}}}}$.

After the optomechanical interaction with the second mechanical oscillator, the light passes through the optical loss channels ${\mathcal{E}}_{{\eta }_{\mathtt{\text{cav}}}}$ and ${\mathcal{E}}_{{\eta }_{\mathtt{\text{det}}}}$, with an optional squeezing operation in between. Finally, a homodyne measurement is made on the phase quadrature of light which contains information about an EPR-like quadrature variable of X _m1 and X _m2. Hence, this measurement projects the mechanical subspace into an entangled state via Bayesian inference of a joint-mechanical quadrature—see appendix E for a more complete discussion. This non-interferometric entanglement protocol, followed by free mechanical evolution U _rot (θ, ϕ) and mechanical decoherence ${\mathcal{E}}_{\gamma }$, is represented by equation (19).

$$\begin{equation}\rho \to {\mathcal{E}}_{\gamma }\left({U}_{\mathtt{\text{rot}}}\left(\theta ,\phi \right)\vert {P}_{\mathtt{\text{l}}}\rangle \langle {P}_{\mathtt{\text{l}}}\vert {\circ}{\mathcal{E}}_{{\eta }_{\mathtt{\text{det}}}}\left({U}_{\mathtt{\text{sq}}}{\circ}\right.\right.\left.\left.{\mathcal{E}}_{{\eta }_{\mathtt{\text{cav}}}}\left({U}_{\mathtt{\text{om}}}^{\prime }{\circ}{\mathcal{E}}_{{\eta }_{\mathtt{\text{cav}}}}\left({U}_{\mathtt{\text{sq}}}{\circ}{\mathcal{E}}_{{\eta }_{\mathtt{\text{cav}}}}\left({U}_{\mathtt{\text{om}}}{\circ}\rho \right)\right)\right)\right)\right).\end{equation} \tag{ 19 }$$

Focussing on the verification procedure, the A _ver and B _ver blocks of the verified covariance matrix σ _ver are obtained by addressing each mechanical oscillator individually with a pulse of light, which is achieved by the use of optical switches as shown in figure 5. The A _ver and B _ver blocks are therefore given by an equivalent expression to equation (12). A single verification pulse, with an initial phase quadrature P _in , is then used to construct C _ver . By allowing both mechanical modes to evolve over a phase-space angle θ and implementing a controllable delay line between the two mechanical modes, such that the second mechanical oscillator evolves over a total phase-space angle ϕ, the phase quadrature entering the homodyne detector is given by P _v (θ, ϕ) = P _in + χ(X _m1(θ) + X _m2(ϕ)). Due to the non-interferometric arrangement, we have the condition ϕ ⩾ θ, which leads to further decay of the C _ver -elements—on top of the decay due to conservative in time approach. By defining (θ, ϕ) = Var(P _v (θ, ϕ)) − Var(P _v (θ, ϕ + π)), we find that the elements describing the mechanical correlations are given by

$$\begin{equation}{C}_{\mathtt{\text{ver}}}=\frac{1}{4{\chi }^{2}}\left(\begin{pmatrix}{cc}{\epsilon}\left(2\pi ,2\pi \right)\hfill & \hfill {\epsilon}\left(2\pi ,\frac{5\pi }{2}\right)\hfill \\ \hfill -{\epsilon}\left(\frac{5\pi }{2},3\pi \right)\hfill & \hfill {\epsilon}\left(\frac{5\pi }{2},\frac{5\pi }{2}\right)\hfill \end{pmatrix}\right).\end{equation} \tag{ 20 }$$

Plots analysing the behaviour of the logarithmic negativity obtained using the non-interferometric scheme are shown in figure 6, and we see that the protocol is relatively insensitive to the squeezing operations as compared to the optical–mechanical and interferometric entanglement schemes. This difference in sensitivity occurs because there is an additional optomechanical interaction and loss channel which further modifies the same optical phase-space distribution towards the ideal, slightly asymmetric, state with Var(P _l ) > Var(X _l ). This additional modification to the phase-space distribution means a lower squeezing parameter r is needed to protect the state from optical loss and mechanical decoherence. The lower plot of figure 6 shows the logarithmic negativity which may be measured in the verification process, compared to the logarithmic negativity of the unverified state at different times since entanglement generation.

Zoom In Zoom Out Reset image size

In the Heisenberg picture, the action of a Gaussian unitary U _g on the 2n-dimensional quadrature vector X is described by the equivalent transformation S X. Here, S belongs to the real symplectic group—$S\in Sp\left(2n,\mathbb{R}\right)$. Therefore, under the action of a Gaussian unitary operation, the first moments and the covariance matrix of a state transform as follows:

$$\begin{equation}\langle \mathbf{X}\rangle \to S\langle \mathbf{X}\rangle \end{equation} \tag{ A.1 }$$

$$\begin{equation}\sigma \to S\sigma {S}^{\text{T}}.\end{equation} \tag{ A.2 }$$

A necessary and sufficient condition for the covariance matrix σ to represent a Gaussian state is given by the uncertainty relation σ + iΩ/2 ⩾ 0 [91]. As the symplectic matrix S preserves the symplectic form SΩS^T = Ω, the symplectic transform preserves the Gaussian character of the quantum state.

A.1.1. Single-mode transformation

Consider a single mode described by the quadrature vector X = (X,P)^T, with $X=\left(a+{a}^{{\dagger}}\right)/\sqrt{2}$ and $P=-\mathrm{i}\left(a-{a}^{{\dagger}}\right)/\sqrt{2}$. Below, we list the symplectic single-mode transformations used in this work.

Rotation—Evolution over time t = θ/ω under the free Hamiltonian H = ℏωa^† a, corresponds to a rotation by an angle θ in phase space, and has the symplectic matrix

$$\begin{equation}{S}_{\mathtt{\text{rot}}}\left(\theta \right)=\left(\begin{pmatrix}{cc}\mathrm{cos} \theta \hfill & \hfill \mathrm{sin} \theta \hfill \\ \hfill -\mathrm{sin} \theta \hfill & \hfill \mathrm{cos} \theta \hfill \end{pmatrix}\right).\end{equation} \tag{ A.3 }$$

Squeezer—A squeezing operation is described by the unitary operator ${U}_{\mathtt{\text{sq}}}=\mathrm{exp}\left\{\frac{1}{2}\left[r{\left({a}^{{\dagger}}\right)}^{2}-{r}^{{\ast}}{a}^{2}\right]\right\}$. For $r\in \mathbb{R}$, the corresponding symplectic matrix is given by

$$\begin{equation}{S}_{\mathtt{\text{sq}}}\left(r\right)=\left(\begin{pmatrix}{cc}\hfill {\mathrm{e}}^{r}\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill {\mathrm{e}}^{-r}\hfill \end{pmatrix}\right).\end{equation} \tag{ A.4 }$$

A.1.2. Two-mode transformations

A two-mode quadrature vector may be written as $\mathbf{X}={\left({X}_{1},{P}_{1},{X}_{2},{P}_{2}\right)}^{\text{T}}$. Here, ${X}_{i}=\left({a}_{i}+{a}_{i}^{{\dagger}}\right)/\sqrt{2}$, ${P}_{i}=-\mathrm{i}\left({a}_{i}-{a}_{i}^{{\dagger}}\right)/\sqrt{2}$, and i = 1, 2.

Beamsplitter—The unitary operation ${U}_{\mathtt{\text{bs}}}=\mathrm{exp}\left\{\alpha \;\mathrm{cos}\;\beta \left({a}_{1}^{{\dagger}}{a}_{2}-{a}_{1}{a}_{2}^{{\dagger}}\right)+\mathrm{i}\alpha \;\mathrm{sin}\;\beta \left({a}_{1}^{{\dagger}}{a}_{2}+{a}_{1}{a}_{2}^{{\dagger}}\right)\right\}$ acts as a general beamsplitter operation and corresponds to the symplectic transformation

$$\begin{equation}{S}_{\mathtt{\text{bs}}}\left(\alpha ,\beta \right)=\left(\begin{pmatrix}{cccc}\mathrm{cos}\;\alpha \hfill & 0\hfill & \hfill \mathrm{sin}\;\alpha \;\mathrm{cos}\;\beta \hfill & \hfill -\mathrm{sin}\;\alpha \;\mathrm{sin}\;\beta \hfill \\ 0\hfill & \hfill \mathrm{cos}\;\alpha \hfill & \hfill \mathrm{sin}\;\alpha \;\mathrm{sin}\;\beta \hfill & \hfill \mathrm{sin}\;\alpha \;\mathrm{cos}\;\beta \hfill \\ \hfill -\mathrm{sin}\;\alpha \;\mathrm{cos}\;\beta \hfill & \hfill -\mathrm{sin}\;\alpha \;\mathrm{sin}\;\beta \hfill & \hfill \mathrm{cos}\;\alpha \hfill & 0\hfill \\ \hfill \mathrm{sin}\;\alpha \;\mathrm{sin}\;\beta \hfill & \hfill -\mathrm{sin}\;\alpha \;\mathrm{cos}\;\beta \hfill & 0\hfill & \hfill \mathrm{cos}\;\alpha \hfill \end{pmatrix}\right).\end{equation} \tag{ A.5 }$$

Linearized-pulsed optomechanical interaction—The optomechanical interaction described by ${\mathrm{e}}^{\mathrm{i}\chi {X}_{1}{X}_{2}}$ corresponds to the symplectic transformation

$$\begin{equation}{S}_{\mathtt{\text{om}}}=\left(\begin{pmatrix}{cccc}\hfill 1\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 1\hfill & \hfill \chi \hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill 1\hfill & \hfill 0\hfill \\ \hfill \chi \hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 1\hfill \end{pmatrix}\right).\end{equation} \tag{ A.6 }$$

A.1.3. Decoherence

A phase-insensitive Gaussian decoherence channel acts on the quantum state ρ according to ${\mathcal{E}}_{\mathtt{\text{g}}}\left(\rho \right)$. The channel ${\mathcal{E}}_{\mathtt{\text{g}}}\left(\rho \right)$ is the result of carrying out a trace operation on a phase-insensitive bath after a system–bath unitary interaction. At the level of the first moments vector ⟨X⟩ and covariance matrix σ, the corresponding phase-space mapping is of the form [69]

$$\begin{equation}\langle \mathbf{X}\rangle \to {G}^{1/2}\langle \mathbf{X}\rangle \end{equation} \tag{ A.7 }$$

$$\begin{equation}\sigma \to {G}^{1/2}\sigma {G}^{1/2}+\left(\mathbb{1}-G\right){\sigma }_{\mathtt{\text{env}}}.\end{equation} \tag{ A.8 }$$

This map allows one to model both optical and mechanical decoherence in the entanglement protocols. The description of the optical loss channel ${\mathcal{E}}_{\eta }$ may be achieved by setting $G=\eta {\mathbb{1}}_{2}$ and σ _env to the vacuum optical state. Here, η is the intensity efficiency of the loss channel. Describing the mechanical decoherence channel ${\mathcal{E}}_{\gamma }$ is accomplished by letting $G={\mathrm{e}}^{-\gamma t}{\mathbb{1}}_{2}$ and σ _env be a thermal state with mean occupation $\overline{N}$. Here, the decoherence channel acts on the mechanical modes for a time t and γ is the intrinsic mechanical decay rate.

A.1.4. Measurements

In Hilbert space, we may describe a Gaussian measurement with the Krauss operator K _g , such that after the measurement the state ρ transforms to a state proportional to K _g ◦ρ. By considering the first-moments vector ⟨X⟩ and the covariance matrix of equation (5), a Gaussian measurement on the optical mode can be described through the mapping

$$\begin{equation}\langle {\mathbf{X}}_{\mathtt{\text{m}}}\rangle \to \langle {\mathbf{X}}_{\mathtt{\text{m}}}\rangle +{C}^{\text{T}}{\left(A+{\sigma }_{\text{meas}}\right)}^{-1}\left(\langle {\mathbf{X}}_{\text{meas}}\rangle -\langle {\mathbf{X}}_{\mathtt{\text{l}}}\rangle \right),\end{equation} \tag{ A.9 }$$

$$\begin{equation}B\to B-{C}^{\text{T}}{\left(A+{\sigma }_{\text{meas}}\right)}^{-1}C.\end{equation} \tag{ A.10 }$$

Here, ⟨X_meas⟩ corresponds to the measurement outcome and σ_meas describes the Gaussian measurement [70]. In the limit of a projection onto an infinitely squeezed Gaussian state, which constitutes an optical homodyne measurement, the term ${\left(A+{\sigma }_{\text{meas}}\right)}^{-1}$ becomes ${\left({{\Pi}}_{\varphi }A{{\Pi}}_{\varphi }\right)}^{\text{MP}}$. Here, MP refers to the Moore–Penrose pseudo-inverse and Π_φ is the projection operator onto the X _l (φ) quadrature

$$\begin{equation}{{\Pi}}_{\varphi }=\left(\begin{matrix}{cc}{\mathrm{cos}}^{2}\;\varphi \hfill & \hfill \mathrm{cos}\;\varphi \;\mathrm{sin}\;\varphi \hfill \\ \hfill \mathrm{cos}\;\varphi \;\mathrm{sin}\;\varphi \hfill & \hfill {\mathrm{sin}}^{2}\;\varphi \hfill \end{matrix}\right).\end{equation} \tag{ A.11 }$$

The Krauss operator corresponding to Π_φ is given by the Hilbert-space projector |X _l (φ)⟩⟨X _l (φ)|, with ${\int }_{-\infty }^{+\infty }\mathrm{d}{X}_{\mathtt{\text{l}}}\left(\varphi \right)\vert {X}_{\mathtt{\text{l}}}\left(\varphi \right)\rangle \langle {X}_{\mathtt{\text{l}}}\left(\varphi \right)\vert =\mathbb{1}$.

Equation (9) describes the precooling stage of the mechanical oscillators. The symplectic formalism outlined above may be used to neatly compute the analytic result of this map. Namely, if the initial state of the optical–mechanical system is given by the covariance matrix

$$\begin{equation}{\sigma }_{\mathtt{\text{in}}}=\frac{1}{2}\left(\begin{pmatrix}{cccc}\hfill 1\hfill & \hfill 0\hfill & 0\hfill & 0\hfill \\ \hfill 0\hfill & \hfill 1\hfill & 0\hfill & 0\hfill \\ 0\hfill & 0\hfill & \hfill 1+2\overline{n}\hfill & 0\hfill \\ 0\hfill & 0\hfill & 0\hfill & \hfill 1+2\bar{n}\hfill \end{pmatrix}\right),\end{equation} \tag{ A.12 }$$

where $\overline{n}$ is the initial mechanical occupation, then the precooled state of the mechanics is given by

$$\begin{equation}{\sigma }_{\mathtt{\text{cool}}}=\left(\begin{pmatrix}{cc}\hfill {V}_{\text{x}}\hfill & \hfill 0\hfill \\ 0\hfill & \hfill {V}_{\text{p}}\hfill \end{pmatrix}\right).\end{equation} \tag{ A.13 }$$

Here,

$$\begin{equation}2{V}_{\text{x}}=1+2\overline{N}+2\;{\mathrm{e}}^{-\gamma t}\left(\overline{n}-\overline{N}\right)+{\mathrm{e}}^{-\gamma t}{\chi }^{2}\end{equation} \tag{ A.14 }$$

$$\begin{equation}2{V}_{\text{p}}=\left(1+2\overline{N}\right)\left(1-{\mathrm{e}}^{-\gamma t}\right)+{\mathrm{e}}^{-\gamma t}\frac{1+2\overline{n}}{1+\eta {\chi }^{2}\left(1+2\overline{n}\right)}.\end{equation} \tag{ A.15 }$$

The preparation of an entangled state of light and a precooled mechanical oscillator, followed by optical and mechanical decoherence, is described by equation (16). Taking $\frac{1}{2}{\mathbb{1}}_{2}\oplus {\sigma }_{\mathtt{\text{cool}}}$ as the initial covariance matrix of light and mechanics, then the entangled state has the covariance matrix

$$\begin{equation}{\sigma }_{\mathtt{\text{lm}}}=\left(\begin{pmatrix}{ccc}\hfill {A}_{\mathtt{\text{lm}}}\hfill & \hfill \hfill & \hfill {C}_{\mathtt{\text{lm}}}\hfill \\ \hfill {C}_{\mathtt{\text{lm}}}^{\text{T}}\hfill & \hfill \hfill & \hfill {B}_{\mathtt{\text{lm}}}\hfill \end{pmatrix}\right),\end{equation} \tag{ A.16 }$$

with

$$\begin{equation}{A}_{\mathtt{\text{lm}}}=\left(\begin{pmatrix}{cc}\hfill {a}_{11}\hfill & \hfill 0\hfill \\ 0\hfill & \hfill {a}_{22}\hfill \end{pmatrix}\right),\end{equation} \tag{ A.17 }$$

$$\begin{equation}{B}_{\mathtt{\text{lm}}}=\left(\begin{pmatrix}{cc}\hfill {b}_{11}\hfill & \hfill 0\hfill \\ 0\hfill & \hfill {b}_{22}\hfill \end{pmatrix}\right),\end{equation} \tag{ A.18 }$$

$$\begin{equation}{C}_{\mathtt{\text{lm}}}=\left(\begin{pmatrix}{cc}0\hfill & \hfill {c}_{12}\hfill \\ \hfill {c}_{21}\hfill & \hfill 0\hfill \end{pmatrix}\right),\end{equation} \tag{ A.19 }$$

where

$$\begin{equation}2{a}_{11}=1-{\eta }_{\mathtt{\text{det}}}+{\mathrm{e}}^{2r}{\eta }_{\mathtt{\text{det}}}\end{equation} \tag{ A.20 }$$

$$\begin{equation}2{a}_{22}=1-{\eta }_{\mathtt{\text{det}}}+{\mathrm{e}}^{-2r}{\eta }_{\mathtt{\text{det}}}\left(1+2{V}_{\text{x}}{\eta }_{\mathtt{\text{cav}}}{\chi }^{2}\right)\end{equation} \tag{ A.21 }$$

$$\begin{equation}2{b}_{11}=\left(1+2\overline{N}\right)\left(1-{\mathrm{e}}^{-\gamma t}\right)+2{V}_{\text{x}}\;{\mathrm{e}}^{-\gamma t}\end{equation} \tag{ A.22 }$$

$$\begin{equation}2{b}_{22}=\left(1+2\overline{N}\right)\left(1-{\mathrm{e}}^{-\gamma t}\right)+\left(2{V}_{\text{p}}+{\chi }^{2}\right){\mathrm{e}}^{-\gamma t}\end{equation} \tag{ A.23 }$$

$$\begin{equation}2{c}_{12}={\mathrm{e}}^{r}\;{\mathrm{e}}^{-\gamma t/2}\sqrt{\eta }\chi \end{equation} \tag{ A.24 }$$

$$\begin{equation}2{c}_{21}=2{V}_{\text{x}}\;{\mathrm{e}}^{-r}\;{\mathrm{e}}^{-\gamma t/2}\sqrt{\eta }\chi .\end{equation} \tag{ A.25 }$$

For brevity, σ _lm is presented in a frame rotating at the mechanical frequency, and we will also present the mechanical entanglement protocols in the same rotating frame. However, as the decoherence channel ${\mathcal{E}}_{\gamma }$ and the unitary operation described by U _rot commute, the form of σ _lm in a non-rotating frame may be obtained by computing ${\mathbb{1}}_{2}\oplus {S}_{\mathtt{\text{rot}}}\left(\theta \right){\sigma }_{\mathtt{\text{lm}}}{\mathbb{1}}_{2}\oplus {S}_{\mathtt{\text{rot}}}^{\text{T}}\left(\theta \right)$.

The protocol that generates an entangled state of two mechanical oscillators using the interferometric scheme is described by equation (17). If the initial state of each mechanical oscillator is described by the covariance matrix σ _cool , the covariance matrix of each light mode interacting with the mechanical modes is $\frac{1}{2}{\mathbb{1}}_{2}$, and the homodyne angles are chosen to be (φ, ψ) = (0, π/2), the final bipartite mechanical covariance matrix is given by

$$\begin{equation}{\sigma }_{\mathtt{\text{int}}}=\left(\begin{pmatrix}{ccc}\hfill {A}_{\mathtt{\text{int}}}\hfill & \hfill \hfill & \hfill {C}_{\mathtt{\text{int}}}\hfill \\ \hfill {C}_{\mathtt{\text{int}}}^{\text{T}}\hfill & \hfill \hfill & \hfill {B}_{\mathtt{\text{int}}}\hfill \end{pmatrix}\right),\end{equation} \tag{ A.26 }$$

$$\begin{equation}{A}_{\mathtt{\text{int}}}={B}_{\mathtt{\text{int}}}=\left(\begin{pmatrix}{cc}\hfill {a}_{11}\hfill & 0\\ 0& \hfill {a}_{22}\hfill \end{pmatrix}\right),\end{equation} \tag{ A.27 }$$

$$\begin{equation}{C}_{\mathtt{\text{int}}}=\left(\begin{pmatrix}{cc}\hfill {c}_{11}\hfill & 0\\ 0& \hfill {c}_{22}\hfill \end{pmatrix}\right),\end{equation} \tag{ A.28 }$$

where

$$\begin{equation*}2{a}_{11}=\left(1+2\overline{N}\right)\left(1-{\mathrm{e}}^{-\gamma t}\right)+2{V}_{\text{x}}\;{\mathrm{e}}^{-\gamma t}\left(1-\frac{\eta {V}_{\text{x}}{\chi }^{2}}{{\mathrm{e}}^{2r}\left(1-{\eta }_{\mathtt{\text{det}}}\right)+{\eta }_{\mathtt{\text{det}}}\left(1+2{V}_{\text{x}}{\eta }_{\mathtt{\text{cav}}}{\chi }^{2}\right)}\right),\end{equation*}$$

$$\begin{align}\hfill 2{a}_{22}& =\left(1+2\overline{N}\right)\left(1-{\mathrm{e}}^{-\gamma t}\right)+{\mathrm{e}}^{-\gamma t}\left[2{V}_{\text{p}}+{\chi }^{2}-\frac{1}{2}\left(\frac{{\mathrm{e}}^{2r}\eta {\chi }^{2}}{1+{\eta }_{\mathtt{\text{det}}}\left({\mathrm{e}}^{2r}-1\right)}\right)\right],\hfill \\ \hfill {c}_{11}& =\frac{{\mathrm{e}}^{-\gamma t}\eta {V}_{\text{x}}^{2}{\chi }^{2}}{{\mathrm{e}}^{2r}\left(1-{\eta }_{\mathtt{\text{det}}}\right)+{\eta }_{\mathtt{\text{det}}}\left(1+2{V}_{\text{x}}{\eta }_{\mathtt{\text{cav}}}{\chi }^{2}\right)},\hfill \\ \hfill {c}_{22}& =-\frac{1}{4}\left(\frac{{\mathrm{e}}^{2r}\;{\mathrm{e}}^{-\gamma t}\eta {\chi }^{2}}{1+{\eta }_{\mathtt{\text{det}}}\left({\mathrm{e}}^{2r}-1\right)}\right).\hfill \end{align}$$

Note that as A _int = B _int , the state is symmetric and so, in addition to the logarithmic negativity, the entanglement of formation may be also computed.

Equation (19) describes the generation of an entangled state of two mechanical oscillators using the non-interferometric scheme. As in the previous subsection, if the initial state of each mechanical oscillator is described by the covariance matrix σ _cool , and the initial covariance matrix of the light mode is $\frac{1}{2}{\mathbb{1}}_{2}$, then the final bipartite mechanical covariance matrix is

$$\begin{equation}{\sigma }_{\mathtt{\text{non}}}=\left(\begin{pmatrix}{ccc}\hfill {A}_{\mathtt{\text{non}}}\hfill & \hfill \hfill & \hfill {C}_{\mathtt{\text{non}}}\hfill \\ \hfill {C}_{\mathtt{\text{non}}}^{\text{T}}\hfill & \hfill \hfill & \hfill {B}_{\mathtt{\text{non}}}\hfill \end{pmatrix}\right),\end{equation} \tag{ A.29 }$$

$$\begin{equation}{A}_{\mathtt{\text{non}}}=\left(\begin{pmatrix}{cc}\hfill {a}_{11}\hfill & \hfill 0\hfill \\ 0& \hfill {a}_{22}\hfill \end{pmatrix}\right),\end{equation} \tag{ A.30 }$$

$$\begin{equation}{B}_{\mathtt{\text{non}}}=\left(\begin{pmatrix}{cc}\hfill {b}_{11}\hfill & \hfill 0\hfill \\ 0& \hfill {b}_{22}\hfill \end{pmatrix}\right),\end{equation} \tag{ A.31 }$$

$$\begin{equation}{C}_{\mathtt{\text{non}}}=\left(\begin{pmatrix}{cc}\hfill {c}_{11}\hfill & \hfill 0\hfill \\ 0& \hfill {c}_{22}\hfill \end{pmatrix}\right),\end{equation} \tag{ A.32 }$$

where

$$\begin{align*}\hfill 2{a}_{11}& =\left(1+2\overline{N}\right)\left(1-{\mathrm{e}}^{-\gamma t}\right)+{\mathrm{e}}^{-\gamma t}{V}_{\text{x}}{J}_{a}\hfill \\ \hfill 2{a}_{22}& =\left(1+2\overline{N}\right)\left(1-{\mathrm{e}}^{-\gamma t}\right)+{\mathrm{e}}^{-\gamma t}\left(2{V}_{\text{p}}+{\chi }^{2}\right)\hfill \\ \hfill 2{b}_{11}& =\left(1+2\overline{N}\right)\left(1-{\mathrm{e}}^{-\gamma t}\right)+{\mathrm{e}}^{-\gamma t}{V}_{\text{x}}{J}_{b}\hfill \\ \hfill 2{b}_{22}& =\left(1+2\overline{N}\right)\left(1-{\mathrm{e}}^{-\gamma t}\right)+{\mathrm{e}}^{-\gamma t}\left[2{V}_{\text{p}}+{\eta }_{\mathtt{\text{cav}}}^{2}{\chi }^{2}\left[1+{\eta }_{\mathtt{\text{cav}}}\left({\mathrm{e}}^{2r}-1\right)\right]\right]\hfill \\ \hfill {c}_{11}& =-\frac{1}{8}\;{\mathrm{e}}^{2r}{V}_{\text{x}}^{2}\;{\mathrm{e}}^{-\gamma t}{\eta }_{\mathtt{\text{cav}}}^{3}{\eta }_{\mathtt{\text{det}}}{\chi }^{2}{J}_{c}\hfill \\ \hfill {c}_{22}& =\frac{1}{2}\;{\mathrm{e}}^{r}\;{\mathrm{e}}^{-\gamma t}{\eta }_{\mathtt{\text{cav}}}^{2}{\chi }^{2}.\hfill \end{align*}$$

The functions J_a , J_b , and J_c are given by

$$\begin{align*}\hfill {J}_{a}& =\left\{{\mathrm{e}}^{8r}{\left(1-{\eta }_{\mathtt{\text{det}}}\right)}^{2}+2\;{\mathrm{e}}^{6r}{\eta }_{\mathtt{\text{det}}}\left(1-{\eta }_{\mathtt{\text{det}}}\right)\left(1-{\eta }_{\mathtt{\text{cav}}}^{2}+2{V}_{\text{x}}{\eta }_{\mathtt{\text{cav}}}^{3}{\chi }^{2}\right)\right.\hfill \\ \hfill & \left.\quad +{\mathrm{e}}^{4r}{\eta }_{\mathtt{\text{det}}}\left\{{\eta }_{\mathtt{\text{det}}}+{\eta }_{\mathtt{\text{cav}}}^{2}\left[2-\left(4-{\eta }_{\mathtt{\text{cav}}}^{2}\right){\eta }_{\mathtt{\text{det}}}+2{V}_{\text{x}}{\eta }_{\mathtt{\text{cav}}}\left(1+\left(1-2{\eta }_{\mathtt{\text{cav}}}^{2}\right){\eta }_{\mathtt{\text{det}}}\right){\chi }^{2}\right.\right.\right.\hfill \\ \hfill & \left.\left.\left. +4{V}_{\text{x}}^{2}{\eta }_{\mathtt{\text{cav}}}^{4}{\eta }_{\mathtt{\text{det}}}{\chi }^{4}\right]\right\}+{\mathrm{e}}^{2r}{\eta }_{\mathtt{\text{cav}}}^{2}{\eta }_{\mathtt{\text{det}}}^{2}\left\{2+{\eta }_{\mathtt{\text{cav}}}\left[2{V}_{\text{x}}{\chi }^{2}-2{\eta }_{\mathtt{\text{cav}}}\left(1-{V}_{\text{x}}{\eta }_{\mathtt{\text{cav}}}{\chi }^{2}\right.\right.\right.\right.\hfill \\ \hfill & \left.\left.\left. {\times}\left(1+2{V}_{\text{x}}{\eta }_{\mathtt{\text{cav}}}{\chi }^{2}\right)\right)\right]\right\}+\left.{\eta }_{\mathtt{\text{cav}}}^{4}{\eta }_{\mathtt{\text{det}}}^{2}\left[1+2{V}_{\text{x}}{\eta }_{\mathtt{\text{cav}}}{\chi }^{2}\right]\right\}/\left\{{\mathrm{e}}^{8r}{\left(1-{\eta }_{\mathtt{\text{det}}}\right)}^{2}\right.\hfill \\ \hfill & \left.\quad +2\;{\mathrm{e}}^{6r}{\eta }_{\mathtt{\text{det}}}\left(1-{\eta }_{\mathtt{\text{det}}}\right)\left(1-{\eta }_{\mathtt{\text{cav}}}^{2}+2{V}_{\text{x}}{\eta }_{\mathtt{\text{cav}}}^{3}{\chi }^{2}\right)\right.\hfill \\ \hfill & \left.\quad +{\mathrm{e}}^{4r}{\eta }_{\mathtt{\text{det}}}\left\{{\eta }_{\mathtt{\text{det}}}+{\eta }_{\mathtt{\text{cav}}}^{2}\left[2+4{V}_{\text{x}}{\eta }_{\mathtt{\text{cav}}}{\chi }^{2}-{\eta }_{\mathtt{\text{det}}}\left(4-{\eta }_{\mathtt{\text{cav}}}^{2}\left(1-4{V}_{\text{x}}{\chi }^{2}\right.\right.\right.\right.\right.\hfill \\ \hfill & \left.\left.\left.\left.\left.{\times}\left(1-{V}_{\text{x}}{\eta }_{\mathtt{\text{cav}}}{\chi }^{2}\right)\right)\right)\right]\right\}+{\mathrm{e}}^{2r}{\eta }_{\mathtt{\text{cav}}}^{2}{\eta }_{\mathtt{\text{det}}}^{2}\left[2-2{\eta }_{\mathtt{\text{cav}}}^{2}+4{V}_{\text{x}}{\eta }_{\mathtt{\text{cav}}}{\chi }^{2}+8{V}_{\text{x}}^{2}{\eta }_{\mathtt{\text{cav}}}^{4}{\chi }^{4}\right]\right.\hfill \\ \hfill & \quad +\left.{\eta }_{\mathtt{\text{cav}}}^{4}{\eta }_{\mathtt{\text{det}}}^{2}\left[1+4{V}_{\text{x}}{\eta }_{\mathtt{\text{cav}}}{\chi }^{2}\left(1+{V}_{\text{x}}{\eta }_{\mathtt{\text{cav}}}{\chi }^{2}\right)\right]\right\},\hfill \end{align*}$$

$$\begin{align*}\hfill {J}_{b}& =\left\{{\mathrm{e}}^{8r}{\left(1-{\eta }_{\mathtt{\text{det}}}\right)}^{2}+2\;{\mathrm{e}}^{6r}{\eta }_{\mathtt{\text{det}}}\left(1-{\eta }_{\mathtt{\text{det}}}\right)\left(1-{\eta }_{\mathtt{\text{cav}}}^{2}+{V}_{\text{x}}{\eta }_{\mathtt{\text{cav}}}^{3}\right)\right.\hfill \\ \hfill & \quad +{\mathrm{e}}^{4r}{\eta }_{\mathtt{\text{det}}}\left\{{\eta }_{\mathtt{\text{det}}}+{\eta }_{\mathtt{\text{cav}}}^{2}\left[2-\left(4-{\eta }_{\mathtt{\text{cav}}}^{2}\right){\eta }_{\mathtt{\text{det}}}+2{V}_{\text{x}}{\eta }_{\mathtt{\text{cav}}}\left(2-\left(1+{\eta }_{\mathtt{\text{cav}}}^{2}\right){\eta }_{\mathtt{\text{det}}}\right){\chi }^{2}\right]\right\}\hfill \\ \hfill & \quad +{\mathrm{e}}^{2r}{\eta }_{\mathtt{\text{cav}}}^{2}{\eta }_{\mathtt{\text{det}}}^{2}\left\{2+{\eta }_{\mathtt{\text{cav}}}\left[4{V}_{\text{x}}{\chi }^{2}-2{\eta }_{\mathtt{\text{cav}}}\left(1+{V}_{\text{x}}{\eta }_{\mathtt{\text{cav}}}{\chi }^{2}\left(1-2{V}_{\text{x}}{\eta }_{\mathtt{\text{cav}}}{\chi }^{2}\right)\right)\right]\right\}\hfill \\ \hfill & \quad +\left.{\eta }_{\mathtt{\text{cav}}}^{4}{\eta }_{\mathtt{\text{det}}}^{2}\left[1+4{V}_{\text{x}}{\eta }_{\mathtt{\text{cav}}}{\chi }^{2}\left(1+{V}_{\text{x}}{\eta }_{\mathtt{\text{cav}}}{\chi }^{2}\right)\right]\right\}/\left\{{\mathrm{e}}^{8r}{\left(1-{\eta }_{\mathtt{\text{det}}}\right)}^{2}\right.\hfill \\ \hfill & \quad +2{\mathrm{e}}^{6r}{\eta }_{\mathtt{\text{det}}}\left(1-{\eta }_{\mathtt{\text{det}}}\right)\left(1-{\eta }_{\mathtt{\text{cav}}}^{2}+2{V}_{\text{x}}{\eta }_{\mathtt{\text{cav}}}^{3}\right)+{\mathrm{e}}^{4r}{\eta }_{\mathtt{\text{det}}}\left\{{\eta }_{\mathtt{\text{det}}}+{\eta }_{\mathtt{\text{cav}}}^{2}\left[2+4{V}_{\text{x}}{\eta }_{\mathtt{\text{cav}}}{\chi }^{2}\right.\right.\hfill \\ \hfill & \left.\left. -{\eta }_{\mathtt{\text{det}}}\left(4-{\eta }_{\mathtt{\text{cav}}}^{2}\left(1-4{V}_{\text{x}}{\eta }_{\mathtt{\text{cav}}}{\chi }^{2}\left(1-{V}_{\text{x}}{\eta }_{\mathtt{\text{cav}}}{\chi }^{2}\right)\right)\right)\right]\right\}\hfill \\ \hfill & \quad +{\mathrm{e}}^{2r}{\eta }_{\mathtt{\text{cav}}}^{2}{\eta }_{\mathtt{\text{det}}}^{2}\left[2-2{\eta }_{\mathtt{\text{cav}}}^{2}+4{V}_{\text{x}}{\eta }_{\mathtt{\text{cav}}}{\chi }^{2}+8{V}_{\text{x}}^{2}{\eta }_{\mathtt{\text{cav}}}^{4}{\chi }^{4}\right]\hfill \\ \hfill & \quad +\left.{\eta }_{\mathtt{\text{cav}}}^{4}{\eta }_{\mathtt{\text{det}}}^{2}\left[1+4{V}_{\text{x}}{\eta }_{\mathtt{\text{cav}}}{\chi }^{2}\left(1+{V}_{\text{x}}{\eta }_{\mathtt{\text{cav}}}{\chi }^{2}\right)\right]\right\},\hfill \end{align*}$$

$$\begin{equation*}{J}_{c}=\frac{4}{{\mathrm{e}}^{5r}\left(1-{\eta }_{\mathtt{\text{det}}}\right)+{\mathrm{e}}^{3r}{\eta }_{\mathtt{\text{det}}}\left[1-{\eta }_{\mathtt{\text{cav}}}^{2}\left(1-2{V}_{\text{x}}{\eta }_{\mathtt{\text{cav}}}{\chi }^{2}\right)\right]+{\mathrm{e}}^{r}{\eta }_{\mathtt{\text{cav}}}^{2}{\eta }_{\mathtt{\text{det}}}\left(1+2{V}_{\text{x}}{\eta }_{\mathtt{\text{cav}}}{\chi }^{2}\right)}.\end{equation*}$$

For r = 0, these functions simplify to the form

$$\begin{align*}\hfill {J}_{a}& =\frac{1}{2}\left[1+\frac{1+4{V}_{\text{x}}{\eta }_{\mathtt{\text{cav}}}^{3}{\eta }_{\mathtt{\text{det}}}{\chi }^{2}}{1+8{V}_{\text{x}}{\eta }_{\mathtt{\text{cav}}}^{3}{\eta }_{\mathtt{\text{det}}}{\chi }^{2}\left(1+2{V}_{\text{x}}{\eta }_{\mathtt{\text{cav}}}^{3}{\eta }_{\mathtt{\text{det}}}{\chi }^{2}\right)}\right]\hfill \\ \hfill {J}_{b}& =\frac{1}{2}\left[1+\frac{1+4{V}_{\text{x}}{\eta }_{\mathtt{\text{cav}}}^{3}{\eta }_{\mathtt{\text{det}}}{\chi }^{2}}{1+8{V}_{\text{x}}{\eta }_{\mathtt{\text{cav}}}^{3}{\eta }_{\mathtt{\text{det}}}{\chi }^{2}\left(1+2{V}_{\text{x}}{\eta }_{\mathtt{\text{cav}}}^{3}{\eta }_{\mathtt{\text{det}}}{\chi }^{2}\right)}\right]\hfill \\ \hfill {J}_{c}& =\frac{4}{1+4{V}_{\text{x}}{\eta }_{\mathtt{\text{cav}}}^{3}{\eta }_{\mathtt{\text{det}}}{\chi }^{2}}.\hfill \end{align*}$$

The covariance matrix σ(t) of the entangled state after a time t since entanglement generation is given by

$$\begin{equation}\sigma \left(t\right)={G}^{1/2}\sigma \left(0\right){G}^{1/2}+\left(\mathbb{1}-G\right){\sigma }_{\mathtt{\text{env}}},\end{equation} \tag{ B.1 }$$

and we define the reconstructed covariance matrix σ _ver so that in the absence of loss σ _ver = σ(0). As in the previous section, we work in a rotating frame at the mechanical frequency only for computational convenience—and hence we do not consider free mechanical evolution in σ(t)—but this in principle is not necessary. We now give two examples of elements defined for σ _ver in the optical–mechanical entanglement scheme. The elements of σ _ver in the mechanical entanglement protocols are defined in precisely the same way.

B _ver,12 and C _ver,11 in the optical–mechanical entanglement scheme—From equation (12), we have that

$$\begin{equation}{B}_{\mathtt{\text{ver}},12}=\frac{1}{2{\chi }^{2}}\left[\mathrm{Var}\left({P}_{\mathtt{\text{v}}}\left(\frac{\pi }{4}\right)\right)-\mathrm{Var}\left({P}_{\mathtt{\text{v}}}\left(\frac{3\pi }{4}\right)\right)\right],\end{equation} \tag{ B.2 }$$

and by using P _v (θ) = P _in + χX _m (θ) one finds that, in the absence of decoherence,

$$\begin{align}\hfill {B}_{\mathtt{\text{ver}},12}& =\frac{1}{2{\chi }^{2}}\left[\mathrm{Var}\left({P}_{\mathtt{\text{in}}}+\frac{\chi }{\sqrt{2}}\left({X}_{\mathtt{\text{m}}}+{P}_{\mathtt{\text{m}}}\right)\right)-\mathrm{Var}\left({P}_{\mathtt{\text{in}}}+\frac{\chi }{\sqrt{2}}\left({X}_{\mathtt{\text{m}}}-{P}_{\mathtt{\text{m}}}\right)\right)\right]\hfill \\ \hfill & =\frac{1}{2}\langle \left\{{X}_{\mathtt{\text{m}}},{P}_{\mathtt{\text{m}}}\right\}\rangle -\langle {X}_{\mathtt{\text{m}}}\rangle \langle {P}_{\mathtt{\text{m}}}\rangle \hfill \\ \hfill \;& ={B}_{12}\left(0\right).\hfill \end{align} \tag{ B.3 }$$

Similarly, from equation (13) we have that

$$\begin{equation}{C}_{\mathtt{\text{ver}},11}=\frac{1}{2\chi }\left[\mathrm{Var}\left({P}_{\mathtt{\text{v}}}^{{\prime\prime}}\left(2\pi \right)\right)-\mathrm{Var}\left({X}_{\text{l}}^{{\prime\prime}}\left(2\pi \right)\right)\right].\end{equation} \tag{ B.4 }$$

P _v ^''(θ) and X _l ^''(θ) are defined implicitly in the main text and are given by ${{P}_{\mathtt{\text{v}}}}^{{\prime\prime}}\left(\theta \right)=\frac{1}{\sqrt{2}}\left({P}_{\mathtt{\text{v}}}\left(\theta \right)+{X}_{\mathtt{\text{l}}}\right)$ and ${{X}_{\mathtt{\text{l}}}}^{{\prime\prime}}\left(\theta \right)=\frac{1}{\sqrt{2}}\left({X}_{\mathtt{\text{l}}}-{P}_{\mathtt{\text{v}}}\left(\theta \right)\right)$, respectively. Therefore, we have

$$\begin{align}\hfill {C}_{\mathtt{\text{ver}},11}& =\frac{1}{4\chi }\left[\mathrm{Var}\left({P}_{\mathtt{\text{v}}}\left(2\pi \right)+{X}_{\mathtt{\text{l}}}\right)-\mathrm{Var}\left({X}_{\mathtt{\text{l}}}-{P}_{\mathtt{\text{v}}}\left(2\pi \right)\right)\right]\hfill \\ \hfill & =\frac{1}{2}\langle \left\{{X}_{\mathtt{\text{l}}},{X}_{\mathtt{\text{m}}}\right\}\rangle -\langle {X}_{\mathtt{\text{l}}}\rangle \langle {X}_{\mathtt{\text{m}}}\rangle \hfill \\ \hfill \;& ={C}_{11}\left(0\right).\hfill \end{align} \tag{ B.5 }$$

However, the reconstruction process will be imperfect due to decoherence processes and as such σ _ver ≠ σ(0). A conservative in time approach is therefore taken to ensure that entanglement is not overestimated and a conservative in time and noise approach is also taken when Var(P _in ) is unknown and cannot be subtracted.

To describe the effect of decoherence in the reconstruction process and address the fact that all the matrix elements of σ _ver are accessed at different times, we use the elements of σ(t) to calculate those of σ _ver . Describing the elements of σ _ver using those of σ(t) is slightly more complicated when the elements of σ _ver depend on two different mechanical phase-space angles θ—and therefore two different times t = θ/ω _m , as illustrated in the calculation of B _ver,12 below.

B _ver,22 , an element that depends on one time—The B _ver,22 element is given by

$$\begin{align}\hfill {B}_{\mathtt{\text{ver}},22}& =\frac{1}{{\chi }^{2}}\left[\mathrm{Var}\left({P}_{\mathtt{\text{v}}}\left(\frac{\pi }{2}\right)\right)-\mathrm{Var}\left({P}_{\text{in}}\right)\right]\hfill \\ \hfill & =\begin{cases}\quad \quad {B}_{22}\left(\frac{\pi }{2}\right)\quad \text{conservative}\;\text{in}\;\text{time}\quad \hfill \\ {B}_{22}\left(\frac{\pi }{2}\right)+\frac{1}{2{\chi }^{2}}\quad \text{conservative}\;\text{in}\;\text{time}\;\text{and}\;\text{noise}.\quad \hfill \end{cases}\hfill \end{align} \tag{ B.6 }$$

The conservative in time and noise approach corresponds to not subtracting Var(P _in ) and here we assume that vacuum noise contaminates B _ver,22.

B _ver,12 , an element that depends on two times—As above

$$\begin{align}\hfill {B}_{\mathtt{\text{ver}},12}& =\frac{1}{2{\chi }^{2}}\left[\mathrm{Var}\left({P}_{\mathtt{\text{v}}}\left(\frac{\pi }{4}\right)\right)-\mathrm{Var}\left({P}_{\mathtt{\text{v}}}\left(\frac{3\pi }{4}\right)\right)\right]\hfill \\ \hfill & =\frac{1}{4}\left[{B}_{11}\left(\frac{\pi }{4}\right)-{B}_{11}\left(\frac{3\pi }{4}\right)+{B}_{22}\left(\frac{\pi }{4}\right)-{B}_{22}\left(\frac{3\pi }{4}\right)\right.\hfill \\ \hfill & \left.\quad +2{B}_{12}\left(\frac{\pi }{4}\right)+2{B}_{12}\left(\frac{3\pi }{4}\right)\right].\hfill \end{align} \tag{ B.7 }$$

In the absence of decoherence, σ(0) ≡ σ(t) and so B _ver,12 = B₁₂(0).

If we include optical losses in the verification stage, the formulas for σ _ver in the main text are modified. To derive these new forms for σ _ver , we assume that the variance of each loss mode is equal to that of P _in and at each loss channel an optical quadrature Q transforms according to

$$\begin{equation}Q\to \sqrt{\eta }Q+\sqrt{1-\eta }{Q}_{\mathtt{\text{vac}}},\end{equation} \tag{ B.8 }$$

where Q _vac is the environmental quadrature for a given channel and Var(Q _vac ) = Var(P _in ). Below we present how the formulas for σ _ver are modified in each entanglement scheme.

Optical–mechanical entanglement scheme—By including optical losses in the verification stage, the A _ver block is modified as follows:

$$\begin{equation*}{A}_{\mathtt{\text{ver}}}\to \frac{1}{\eta }\left(\begin{pmatrix}{cc}\mathrm{Var}\left({X}_{\mathtt{\text{l}}}\right)-\left(1-\eta \right)\mathrm{Var}\left({P}_{\mathtt{\text{in}}}\right)\hfill & \hfill \frac{1}{2}\left[\mathrm{Var}\left({P}_{\mathtt{\text{l}}}\left(\frac{3\pi }{4}\right)\right)-\mathrm{Var}\left({P}_{\mathtt{\text{l}}}\left(\frac{\pi }{4}\right)\right)\right]\hfill \\ \hfill \frac{1}{2}\left[\mathrm{Var}\left({P}_{\mathtt{\text{l}}}\left(\frac{3\pi }{4}\right)\right)-\mathrm{Var}\left({P}_{\mathtt{\text{l}}}\left(\frac{\pi }{4}\right)\right)\right]\hfill & \hfill \mathrm{Var}\left({P}_{\mathtt{\text{l}}}\right)-\left(1-\eta \right)\mathrm{Var}\left({P}_{\mathtt{\text{in}}}\right)\hfill \end{pmatrix}\right).\end{equation*}$$

The B _ver and C _ver blocks rescale according to B _ver → B _ver /η, and C _ver → C _ver /η.

Interferometric mechanical entanglement scheme—In this scheme, all the elements of the reconstructed covariance matrix rescale in the same way σ _ver → σ _ver /η.

Non-interferometric mechanical entanglement scheme—Here, the A _ver and B _ver blocks are modified in the same way A _ver → A _ver /η and B _ver → B _ver /η. While C _ver is modified to C _ver → C _ver /η².

Assuming that η is a well-known quantity, none of the results in the reconstruction stage are affected.

Once σ _ver has been reconstructed, an inverse map may be applied to calculate the covariance matrix σ(0) of the bipartite state at the moment of entanglement generation. In the verification procedure, the elements of σ _ver are defined in terms of the state σ(t). From equation (A.8), for the elements of σ _ver defined at a single time, we have

$$\begin{align}\hfill {\sigma }_{\mathtt{\text{ver}},ij}& ={\sigma }_{ij}\left({t}_{ij}\right)\hfill \\ \hfill & ={\left[{G}^{1/2}\sigma \left(0\right){G}^{1/2}+\left(\mathbb{1}-G\right){\sigma }_{\mathtt{\text{env}}}\right]}_{ij}.\hfill \end{align} \tag{ B.9 }$$

Here, each element of σ _ver may be defined at different times, hence the inclusion of the indices ij to the temporal parameter t_ij . Using the fact that G and σ _env are diagonal matrices we may simplify the expression for σ _ver,ij to

$$\begin{equation}{\sigma }_{\mathtt{\text{ver}},ij}={G}_{ii}^{1/2}{\sigma }_{ij}\left(0\right){G}_{jj}^{1/2}+{\left(\mathbb{1}-G\right)}_{ii}{\sigma }_{\mathtt{\text{env}},jj}{\delta }_{ij},\end{equation} \tag{ B.10 }$$

or inversely

$$\begin{equation}{\sigma }_{ij}\left(0\right)={G}_{ii}^{-1/2}\left[{\sigma }_{\mathtt{\text{ver}},ij}-{\left(\mathbb{1}-G\right)}_{ii}{\sigma }_{\mathtt{\text{env}},jj}{\delta }_{ij}\right]{G}_{jj}^{-1/2}.\end{equation} \tag{ B.11 }$$

Elements defined at two different times, t_ij,1 and t_ij,2, may be written as

$$\begin{align}\hfill {\sigma }_{\mathtt{\text{ver}},ij}& =\sum _{k}\sum _{nm}{c}_{nm,k}{\sigma }_{nm}\left({t}_{ij,k}\right)\hfill \\ \hfill & =\sum _{k}\sum _{nm\ne ij}{c}_{nm,k}{\sigma }_{nm}\left({t}_{ij,k}\right)+\sum _{k}{c}_{ij,k}{\sigma }_{ij}\left({t}_{ij,k}\right).\hfill \end{align} \tag{ B.12 }$$

Here, c_nm,k are real coefficients and σ_ij (t_ij,k ) is given by

$$\begin{equation}{\sigma }_{ij}\left({t}_{ij,k}\right)={G}_{ii}^{1/2}\left({t}_{ij,k}\right){\sigma }_{ij}\left(0\right){G}_{jj}^{1/2}\left({t}_{ij,k}\right)+{\left(\mathbb{1}-G\right)}_{ii}\left({t}_{ij,k}\right){\sigma }_{\mathtt{\text{env}},jj}{\delta }_{ij},\end{equation} \tag{ B.13 }$$

where a temporal argument has to now be added to the matrix G to facilitate the two different times. Rearranging the expression for σ _ver,ij leads to

$$\begin{equation}{\sigma }_{ij}\left(0\right)=\frac{{\sigma }_{\mathtt{\text{ver}},ij}-\sum _{k}\sum _{nm\ne ij}{c}_{nm,k}{\sigma }_{nm}\left({t}_{ij,k}\right)-\sum _{k}{c}_{ij,k}{\left(\mathbb{1}-G\right)}_{ii}\left({t}_{ij,k}\right){\sigma }_{\mathtt{\text{env}},jj}{\delta }_{ij}}{\sum _{k}{c}_{ij,k}{G}_{ii}^{1/2}\left({t}_{ij,k}\right){G}_{jj}^{1/2}\left({t}_{ij,k}\right)}.\end{equation} \tag{ B.14 }$$

Providing the quantities σ_nm (0), with nm ≠ ij, have been calculated from the elements of σ _ver defined at one time, then the quantities σ_nm (t_ij,k ) may also be calculated. After this, the elements σ_ij (0) can be calculated from the elements of σ _ver that depend on two times. To summarise: in order to calculate σ(0), the inverse map is applied to elements of σ _ver that depend on one time and these are then used to find the elements which depend on two times. Below is a short example of how this procedure works for the B _ver block of the optical–mechanical covariance matrix.

Calculating elements of σ(0) using the inverse map—Consider again the element B _ver,12 given by

$$\begin{align}\hfill {B}_{\mathtt{\text{ver}},12}& =\frac{1}{4}\left[{B}_{11}\left(\frac{\pi }{4}\right)-{B}_{11}\left(\frac{3\pi }{4}\right)+{B}_{22}\left(\frac{\pi }{4}\right)-{B}_{22}\left(\frac{3\pi }{4}\right)\right.+\left.2{B}_{12}\left(\frac{\pi }{4}\right)+2{B}_{12}\left(\frac{3\pi }{4}\right)\right].\hfill \end{align} \tag{ B.15 }$$

The elements B _ver,11 = B₁₁(0) and ${B}_{\mathtt{\text{ver}},22}={B}_{22}\left(\frac{\pi }{2}\right)$ only depend on one time, and so we have that B₁₁(0) and B₂₂(0) are given by equation (B.11). Explicitly we have

$$\begin{equation}{B}_{11}\left(0\right)={B}_{\mathtt{\text{ver}},11},\end{equation} \tag{ B.16 }$$

$$\begin{equation}{B}_{22}\left(0\right)={G}_{22}^{-1/2}\left(\pi /2\right)\left[{B}_{\mathtt{\text{ver}},22}-{\left(\mathbb{1}-G\right)}_{22}{\sigma }_{\mathtt{\text{env}},22}\right]{G}_{22}^{-1/2}\left(\pi /2\right).\end{equation} \tag{ B.17 }$$

This allows one to work out ${B}_{11}\left(\frac{\pi }{4}\right)$, ${B}_{11}\left(\frac{3\pi }{4}\right),{B}_{22}\left(\frac{\pi }{4}\right)$, and ${B}_{22}\left(\frac{3\pi }{4}\right)$ using equation (B.13). Equation (B.14) may then be used to find an expression for B₁₂(0). Namely,

$$\begin{equation}{B}_{12}\left(0\right)=\frac{{B}_{\mathtt{\text{ver}},12}-\frac{1}{4}\left[{B}_{11}\left(\frac{\pi }{4}\right)-{B}_{11}\left(\frac{3\pi }{4}\right)+{B}_{22}\left(\frac{\pi }{4}\right)-{B}_{22}\left(\frac{3\pi }{4}\right)\right]}{\frac{1}{2}{G}_{11}^{1/2}\left(\pi /4\right){G}_{11}^{1/2}\left(\pi /4\right)+\frac{1}{2}{G}_{11}^{1/2}\left(3\pi /4\right){G}_{11}^{1/2}\left(3\pi /4\right)}.\end{equation} \tag{ B.18 }$$

The above analysis assumes that the inverse maps can be applied perfectly, whereas in reality this will not be the case. Uncertainties in the values of η and γ will limit the precision to which these inverse operations can be applied and hence limit the amount of entanglement one can verify when reconstructing σ(0) from σ _ver .

The Krauss operator corresponding to the optical–mechanical entanglement protocol is given simply by

$$\begin{equation}{\Upsilon}\left({X}_{\mathtt{\text{l}}},{X}_{\mathtt{\text{m}}}\right)={U}_{\mathtt{\text{sq}}}\left(r\right){\mathrm{e}}^{\mathrm{i}\chi {X}_{\mathtt{\text{l}}}{X}_{\mathtt{\text{m}}}}.\end{equation} \tag{ E.1 }$$

The two mode unitary ${\mathrm{e}}^{\mathrm{i}\chi {X}_{\mathtt{\text{l}}}{X}_{\mathtt{\text{m}}}}$ entangles the optical and mechanical states, while the optical squeezer protects the state from optical loss and mechanical decoherence.

ϒ(X _l , X _m ) is completely unitary, which reflects the continuous behaviour of the entanglement montone with respect to the system parameters, such as interaction strength χ and optical efficiency η.

The Krauss operator ϒ(X _m , φ) describing the precooling stage results from a single optomechanical interaction, where the input light is in the state |α⟩, followed by optical squeezing and homodyne detection at angle φ. Hence,

$$\begin{align}\hfill {\Upsilon}\left({X}_{\mathtt{\text{m}}},\varphi \right)=\langle {X}_{\mathtt{\text{l}}}\left(\varphi \right)\vert {U}_{\mathtt{\text{sq}}}\left(r\right){\mathrm{e}}^{\mathrm{i}\chi {X}_{\mathtt{\text{l}}}{X}_{\mathtt{\text{m}}}}\vert \alpha \rangle \\ \hfill ={\mathrm{e}}^{\mathrm{i}\chi {X}_{\alpha }{X}_{\mathtt{\text{m}}}/2}\langle {X}_{\mathtt{\text{l}}}\left(\varphi \right)\vert {U}_{\mathtt{\text{sq}}}\left.\left(r\right)\right\vert \left({X}_{\alpha }+\mathrm{i}\chi {X}_{\mathtt{\text{m}}}\right)/\sqrt{2}\rangle .\end{align} \tag{ E.2 }$$

Therefore,

$$\begin{align}\hfill {\Upsilon}\left({X}_{\mathtt{\text{m}}},\varphi \right)={\left({{\Gamma}}_{\mathtt{\text{r}}}/\pi \right)}^{1/4}\;\mathrm{exp}\left\{-{\Gamma}{\left({X}_{\mathtt{\text{l}}}\left(\varphi \right)-{X}_{0}\right)}^{2}/2\right\}\\ \hfill \mathrm{exp}\left\{\mathrm{i}{P}_{0}{X}_{\mathtt{\text{l}}}\left(\varphi \right)\right\}\mathrm{exp}\left\{-\mathrm{i}{X}_{0}{P}_{0}/2\right\}\mathrm{exp}\left\{\mathrm{i}\chi {X}_{\alpha }{X}_{\mathtt{\text{m}}}/2\right\}.\end{align} \tag{ E.3 }$$

Here,

$$\begin{equation}{X}_{0}={\mathrm{e}}^{r}{X}_{\alpha }\;\mathrm{cos}\;\varphi +{\mathrm{e}}^{-r}\chi {X}_{\mathtt{\text{m}}}\;\mathrm{sin}\;\varphi ,\end{equation} \tag{ E.4 }$$

$$\begin{equation}{P}_{0}={\mathrm{e}}^{-r}\chi {X}_{\mathtt{\text{m}}}\;\mathrm{cos}\;\varphi -{\mathrm{e}}^{r}{X}_{\alpha }\;\mathrm{sin}\;\varphi ,\end{equation} \tag{ E.5 }$$

$$\begin{equation}{\Gamma}\left(\varphi \right)={{\Gamma}}_{\mathtt{\text{r}}}\left(\varphi \right)+\mathrm{i}{{\Gamma}}_{\mathtt{\text{i}}}\left(\varphi \right),\end{equation} \tag{ E.6 }$$

$$\begin{equation}{{\Gamma}}_{\mathtt{\text{r}}}\left(\varphi \right)=\frac{1}{\mathrm{cosh}\;2r+\mathrm{cos}\;2\varphi \;\mathrm{sinh}\;2r},\end{equation} \tag{ E.7 }$$

$$\begin{equation}{{\Gamma}}_{\mathtt{\text{i}}}\left(\varphi \right)=\frac{2\;\mathrm{sin}\;2\varphi \;\mathrm{sinh}\;2r}{\mathrm{cosh}\;2r+\mathrm{cos}\;2\varphi \;\mathrm{sinh}\;2r}.\end{equation} \tag{ E.8 }$$

The optimal homodyne angle to precool the mechanical mode via measurement is φ = π/2. This is because the mechanical position quadrature is imprinted on the optical phase quadrature of light. In this case the Krauss operator is

$$\begin{equation}{\Upsilon}\left({X}_{\mathtt{\text{m}}},\pi /2\right)={\left({\mathrm{e}}^{2r}/\pi \right)}^{1/4}\mathrm{exp}\left\{-{\mathrm{e}}^{2r}{\left({P}_{\mathtt{\text{l}}}-{\mathrm{e}}^{-r}\chi {X}_{\mathtt{\text{m}}}\right)}^{2}/2-\mathrm{i}\;{\mathrm{e}}^{r}{X}_{\alpha }{P}_{\mathtt{\text{l}}}\right\}.\end{equation} \tag{ E.9 }$$

We write the action of the Krauss operator on the mechanical mode as ϒ◦ρ. Applying the Krauss operator m times produces the state ϒ^m ◦ρ, where

$$\begin{align}\hfill {\left[{\Upsilon}\left({X}_{\mathtt{\text{m}}},\varphi \right)\right]}^{m}={\left({{\Gamma}}_{\mathtt{\text{r}}}/\pi \right)}^{m/4}\mathrm{exp}\left\{-{\Gamma}{\left(\sqrt{m}{X}_{\mathtt{\text{l}}}\left(\varphi \right)-\sqrt{m}{X}_{0}\right)}^{2}/2\right\}\\ \hfill \mathrm{exp}\left\{\mathrm{i}m{P}_{0}{X}_{\mathtt{\text{l}}}\left(\varphi \right)\right\}\mathrm{exp}\left\{-\mathrm{i}m{X}_{0}{P}_{0}/2\right\}\mathrm{exp}\left\{\mathrm{i}m\chi {X}_{\alpha }{X}_{\mathtt{\text{m}}}/2\right\}.\end{align} \tag{ E.10 }$$

This operation is equivalent to a single application of the Krauss operator but with rescalings $\chi \to \sqrt{m}\chi$, $\alpha \to \sqrt{m}\alpha$, and ${X}_{\mathtt{\text{l}}}\left(\varphi \right)\to \sqrt{m}{X}_{\mathtt{\text{l}}}\left(\varphi \right)$. Hence, increasing the number of pulsed optomechanical interactions by a factor m is equivalent to a single pulsed interaction with m-times more photons.

The Krauss operator which acts on the initially separable state of the two mechanical oscillators, corresponding to equation (17), is

$$\begin{equation}{\Upsilon}\left({X}_{\mathtt{\text{m}}1},{X}_{\mathtt{\text{m}}2},\varphi ,\psi \right)=\langle {X}_{\mathtt{\text{l}}1}\left(\varphi \right)\vert \langle {X}_{\mathtt{\text{l}}2}\left(\psi \right)\vert {U}_{\mathtt{\text{bs}}}{U}_{\mathtt{\text{sq}}}^{\otimes 2}\;{\mathrm{e}}^{\mathrm{i}\chi {X}_{\mathtt{\text{l}}2}{X}_{\mathtt{\text{m}}2}}\;{\mathrm{e}}^{\mathrm{i}\chi {X}_{\mathtt{\text{l}}1}{X}_{\mathtt{\text{m}}1}}{U}_{\mathtt{\text{bs}}}\vert \sqrt{2}\alpha \rangle \vert 0\rangle .\end{equation} \tag{ E.11 }$$

Up to irrelevant phase terms, the Krauss operator is therefore given by

$$\begin{align}\hfill {\Upsilon}\left({X}_{\mathtt{\text{m}}1},{X}_{\mathtt{\text{m}}2},\varphi ,\psi \right)& ={\left({{\Gamma}}_{\mathtt{\text{r}}}\left(\varphi \right){{\Gamma}}_{\mathtt{\text{r}}}\left(\psi \right)/{\pi }^{2}\right)}^{1/4}\mathrm{exp}\left\{-{\Gamma}\left(\varphi \right){\left({X}_{\mathtt{\text{l}}1}\left(\varphi \right)-{X}_{1}\right)}^{2}/2\right\}\hfill \\ \hfill & \quad {\times}\mathrm{exp}\left\{\mathrm{i}{P}_{1}{X}_{\mathtt{\text{l}}1}\left(\varphi \right)\right\}\mathrm{exp}\left\{-\mathrm{i}{X}_{1}{P}_{1}/2\right\}\mathrm{exp}\left\{-{\Gamma}\left(\psi \right)\left({X}_{\mathtt{\text{l}}2}\left(\psi \right)\right.\right.\hfill \\ \hfill & \left.{\left.\quad -{X}_{2}\right)}^{2}/2\right\}\mathrm{exp}\left\{\mathrm{i}{P}_{2}{X}_{\mathtt{\text{l}}2}\left(\psi \right)\right\}\mathrm{exp}\left\{-\mathrm{i}{X}_{2}{P}_{2}/2\right\},\hfill \end{align} \tag{ E.12 }$$

where Γ is defined as above, while

$$\begin{equation}{X}_{1}={\mathrm{e}}^{-r}\chi \left({X}_{\mathtt{\text{m}}1}+{X}_{\mathtt{\text{m}}2}\right)\mathrm{sin}\;\varphi /\sqrt{2},\end{equation} \tag{ E.13 }$$

$$\begin{equation}{P}_{1}={\mathrm{e}}^{-r}\chi \left({X}_{\mathtt{\text{m}}1}+{X}_{\mathtt{\text{m}}2}\right)\mathrm{cos}\;\varphi /\sqrt{2},\end{equation} \tag{ E.14 }$$

$$\begin{equation}{X}_{2}=-\sqrt{2}\;{\mathrm{e}}^{r}{X}_{\alpha }\;\mathrm{cos}\;\psi -{\mathrm{e}}^{-r}\chi \left({X}_{\mathtt{\text{m}}1}-{X}_{\mathtt{\text{m}}2}\right)\mathrm{sin}\;\psi /\sqrt{2},\end{equation} \tag{ E.15 }$$

$$\begin{equation}{P}_{2}=\sqrt{2}\;{\mathrm{e}}^{r}{X}_{\alpha }\;\mathrm{sin}\;\psi -{\mathrm{e}}^{-r}\chi \left({X}_{\mathtt{\text{m}}1}-{X}_{\mathtt{\text{m}}2}\right)\mathrm{cos}\;\psi /\sqrt{2}.\end{equation} \tag{ E.16 }$$

In this case, the polar decomposition of the Krauss operator is given by ϒ(X _m1, X _m2, φ, ψ) = U(X _m1, X _m2, φ, ψ)P(X _m1, X _m2, φ, ψ), with positive operator

$$\begin{align}\hfill P\left({X}_{\mathtt{\text{m}}1},{X}_{\mathtt{\text{m}}2},\varphi ,\psi \right)& ={\left({{\Gamma}}_{\mathtt{\text{r}}}\left(\varphi \right){{\Gamma}}_{\mathtt{\text{r}}}\left(\psi \right)/{\pi }^{2}\right)}^{1/4}\mathrm{exp}\left\{-{{\Gamma}}_{\mathtt{\text{r}}}\left(\varphi \right){\left({X}_{\mathtt{\text{l}}1}\left(\varphi \right)-{X}_{1}\right)}^{2}/2\right\}\hfill \\ \hfill & \quad {\times}\mathrm{exp}\left\{-{{\Gamma}}_{\mathtt{\text{r}}}\left(\psi \right){\left({X}_{\mathtt{\text{l}}2}\left(\psi \right)-{X}_{2}\right)}^{2}/2\right\},\hfill \end{align} \tag{ E.17 }$$

and effective unitary operator

$$\begin{align}\hfill U\left({X}_{\mathtt{\text{m}}1},{X}_{\mathtt{\text{m}}2},\varphi ,\psi \right)& ={\left({{\Gamma}}_{\mathtt{\text{r}}}\left(\varphi \right){{\Gamma}}_{\mathtt{\text{r}}}\left(\psi \right)/{\pi }^{2}\right)}^{1/4}\mathrm{exp}\left\{-\mathrm{i}{{\Gamma}}_{\mathtt{\text{i}}}\left(\varphi \right){\left({X}_{\mathtt{\text{l}}1}\left(\varphi \right)-{X}_{1}\right)}^{2}/2\right\}\hfill \\ \hfill & \quad {\times}\mathrm{exp}\left\{\mathrm{i}{P}_{1}{X}_{\mathtt{\text{l}}1}\left(\varphi \right)\right\}\mathrm{exp}\left\{-\mathrm{i}{X}_{1}{P}_{1}/2\right\}\mathrm{exp}\left\{-\mathrm{i}{{\Gamma}}_{\mathtt{\text{i}}}\left(\psi \right)\left({X}_{\mathtt{\text{l}}2}\left(\psi \right)\right.\right.\hfill \\ \hfill & \left.\quad {\left.-{X}_{2}\right)}^{2}/2\right\}\mathrm{exp}\left\{\mathrm{i}{P}_{2}{X}_{\mathtt{\text{l}}2}\left(\psi \right)\right\}\mathrm{exp}\left\{-\mathrm{i}{X}_{2}{P}_{2}/2\right\},\hfill \end{align} \tag{ E.18 }$$

which also depends on the homodyne measurement outcomes. The parts of these operators which depend on the square of the mechanical EPR quadratures may induce entanglement as they contain terms proportional to the product X _m1 X _m2. We call these parts P _ent and U _ent , respectively, and they are given by

$$\begin{align}\hfill {P}_{\mathtt{\text{ent}}}& ={\left({{\Gamma}}_{\mathtt{\text{r}}}\left(\varphi \right){{\Gamma}}_{\mathtt{\text{r}}}\left(\psi \right)/{\pi }^{2}\right)}^{1/4}\mathrm{exp}\left\{-{{\Gamma}}_{\mathtt{\text{r}}}\left(\varphi \right){\mathrm{e}}^{-2r}{\chi }^{2}{\left({X}_{\mathtt{\text{m}}1}+{X}_{\mathtt{\text{m}}2}\right)}^{2}{\mathrm{sin}}^{2}\;\varphi /4\right\}\hfill \\ \hfill & \quad {\times}\mathrm{exp}\left\{-{{\Gamma}}_{\mathtt{\text{r}}}\left(\psi \right){\mathrm{e}}^{-2r}{\chi }^{2}{\left({X}_{\mathtt{\text{m}}1}-{X}_{\mathtt{\text{m}}2}\right)}^{2}{\mathrm{sin}}^{2}\;\psi /4\right\},\hfill \end{align} \tag{ E.19 }$$

and

$$\begin{align}\hfill {U}_{\mathtt{\text{ent}}}& ={\left({{\Gamma}}_{\mathtt{\text{r}}}\left(\varphi \right){{\Gamma}}_{\mathtt{\text{r}}}\left(\psi \right)/{\pi }^{2}\right)}^{1/4}\mathrm{exp}\left\{-\mathrm{i}{{\Gamma}}_{\mathtt{\text{i}}}\left(\varphi \right){\mathrm{e}}^{-2r}{\chi }^{2}{\left({X}_{\mathtt{\text{m}}1}+{X}_{\mathtt{\text{m}}2}\right)}^{2}{\mathrm{sin}}^{2}\;\varphi /4\right\}\hfill \\ \hfill & \quad {\times}\mathrm{exp}\left\{-\mathrm{i}{{\Gamma}}_{\mathtt{\text{i}}}\left(\psi \right){\mathrm{e}}^{-2r}{\chi }^{2}{\left({X}_{\mathtt{\text{m}}1}-{X}_{\mathtt{\text{m}}2}\right)}^{2}{\mathrm{sin}}^{2}\;\psi /4\right\}\mathrm{exp}\left\{-\mathrm{i}\;{\mathrm{e}}^{-2r}{\chi }^{2}{\left({X}_{\mathtt{\text{m}}1}+{X}_{\mathtt{\text{m}}2}\right)}^{2}\right.\hfill \\ \hfill & \quad \left.{\times}\mathrm{sin}\;\varphi \;\mathrm{cos}\;\varphi /4\right\}\mathrm{exp}\left\{-\mathrm{i}\;{\mathrm{e}}^{-2r}{\chi }^{2}{\left({X}_{\mathtt{\text{m}}1}-{X}_{\mathtt{\text{m}}2}\right)}^{2}\mathrm{sin}\;\psi \;\mathrm{cos}\;\psi /4\right\}.\hfill \end{align} \tag{ E.20 }$$

As stated in the main text, the homodyne configurations (φ, ψ) = (0, π/2) and (φ, ψ) = (π/2, 0) therefore corresponds to Bayesian inference of the mechanical EPR quadratures (X _m1 + X _m2) and (X _m1 − X _m2), respectively. To see this explicitly, consider the case when (φ, ψ) = (0, π/2), then we have that

$$\begin{align}\hfill {P}_{\mathtt{\text{ent}}}=\mathrm{exp}\left\{-{\chi }^{2}{\left({X}_{\mathtt{\text{m}}1}-{X}_{\mathtt{\text{m}}2}\right)}^{2}/4\right\},\\ \hfill \quad \quad {U}_{\mathtt{\text{ent}}}=\mathbb{1}.\end{align} \tag{ E.21 }$$

Hence, the total operation ϒ(X _m1, X _m2, θ, ϕ) corresponds to Bayesian inference up to local unitaries.

Whilst at (φ, ψ) = (π/4, 3π/4) and (φ, ψ) = (3π/4, π/4) the measurement induces a two-mode entangling unitary operation and the positive operator P _ent does not lead to any entanglement generation. From a Bayesian perspective, this operation amounts to unitarily evolving the prior towards a more entangled state. Taking (φ, ψ) = (π/4, 3π/4), we have that

$$\begin{align}\hfill {P}_{\mathtt{\text{ent}}}=\mathrm{exp}\left\{-{\mathrm{e}}^{-2r}{\chi }^{2}\left({X}_{\mathtt{\text{m}}1}^{2}+{X}_{\mathtt{\text{m}}2}^{2}\right)/4\;\mathrm{cosh}\;2r\right\},\\ \hfill {U}_{\mathtt{\text{ent}}}=\mathrm{exp}\left\{-\mathrm{i}\;{\mathrm{e}}^{-2r}{\chi }^{2}\left(1+2\mathrm{tanh}\;2r\right){X}_{\mathtt{\text{m}}1}{X}_{\mathtt{\text{m}}2}/2\right\}.\end{align} \tag{ E.22 }$$

Figure E1 shows the balance between Bayesian inference and unitary action with and without the use of optical squeezers. These contour plots for the logarithmic negativity, show that there is a continuum of optimal homodyne configurations, regardless of whether or not optical squeezers are implemented. The configurations (φ, ψ) = (0, π/2) and (φ, ψ) = (π/2, 0) chosen in the main text correspond to just two of these optimal choices.

Zoom In Zoom Out Reset image size

When optical squeezers are not used the contour plot is more heavily weighted around the set of points (φ, ψ) = {(0, π/2), (π/2, 0), (π, π/2), (π/2, π)}, which corresponds to maximizing the contribution of P _ent , as compared to the regions surrounding (φ, ψ) = {(0, π/2), (3π/4, π/4)}, which corresponds to maximizing the entangling effect of U _ent . This observation demonstrates that Bayesian inference of the mechanical EPR quadratures is the more effective way to create mechanical entanglement from the optical–mechanical entanglement. This is because, without the use of squeezers, the optical subspace is asymmetric with more information about the mechanical position quadratures encoded in the phase quadratures of the optical modes.

It is interesting to note that when the optimal choice of optical squeezers is utilized, the contour plot shows a more equal balance between the effects of P _ent and U _ent . With the best choice of homodyne angles corresponding to approximately ψ = φ − π/2 and ψ = φ + π/2. Optical squeezers bring the optical modes to a more symmetric state, and hence the contour plot is less heavily weighted to areas that correspond to a combination of a phase and an amplitude quadrature measurement.

As mentioned in the main text, no entanglement can be generated when φ = ψ as this allows one to gain knowledge of X _m1 and X _m2 separately.

The Krauss operator that acts on the initial separable mechanical state, corresponding to equation (19), is given by

$$\begin{equation}{\Upsilon}\left({X}_{\mathtt{\text{m}}1},{X}_{\mathtt{\text{m}}2},\varphi \right)=\langle {X}_{\mathtt{\text{l}}}\left(\varphi \right)\vert {U}_{\mathtt{\text{sq}}}\;{\mathrm{e}}^{\mathrm{i}\chi {X}_{\mathtt{\text{l}}}{X}_{\mathtt{\text{m}}2}}{U}_{\mathtt{\text{sq}}}\;{\mathrm{e}}^{\mathrm{i}\chi {X}_{\mathtt{\text{l}}}{X}_{\mathtt{\text{m}}1}}\vert \alpha \rangle .\end{equation} \tag{ E.23 }$$

Up to irrelevant local unitaries, the Krauss operator may be expressed as

$$\begin{align}\hfill {\Upsilon}\left({X}_{\mathtt{\text{m}}1},{X}_{\mathtt{\text{m}}2},\varphi \right)& ={\left({{\Gamma}}_{\mathtt{\text{r}}}/\pi \right)}^{1/4}\mathrm{exp}\left\{-{\Gamma}{\left({X}_{\mathtt{\text{l}}}\left(\varphi \right)-{X}_{0}\right)}^{2}/2\right\}\hfill \\ \hfill & \quad {\times}\mathrm{exp}\left\{\mathrm{i}{P}_{0}{X}_{\mathtt{\text{l}}}\left(\varphi \right)\right\}\mathrm{exp}\left\{-\mathrm{i}{X}_{0}{P}_{0}/2\right\}.\hfill \end{align} \tag{ E.24 }$$

Here,

$$\begin{equation}{X}_{0}={\mathrm{e}}^{2r}{X}_{\alpha }\;\mathrm{cos}\;\varphi +\chi \left({\mathrm{e}}^{-2r}{X}_{\mathtt{\text{m}}1}+{\mathrm{e}}^{-r}{X}_{\mathtt{\text{m}}2}\right)\mathrm{sin}\;\varphi ,\end{equation} \tag{ E.25 }$$

$$\begin{equation}{P}_{0}=\chi \left({\mathrm{e}}^{-2r}{X}_{\mathtt{\text{m}}1}+{\mathrm{e}}^{-r}{X}_{\mathtt{\text{m}}2}\right)\mathrm{cos}\;\varphi -{\mathrm{e}}^{2r}{X}_{\alpha }\;\mathrm{sin}\;\varphi ,\end{equation} \tag{ E.26 }$$

$$\begin{equation}{\Gamma}\left(\varphi \right)={{\Gamma}}_{\mathtt{\text{r}}}\left(\varphi \right)+\mathrm{i}{{\Gamma}}_{\mathtt{\text{i}}}\left(\varphi \right),\end{equation} \tag{ E.27 }$$

$$\begin{equation}{{\Gamma}}_{\mathtt{\text{r}}}\left(\varphi \right)=\frac{1}{\mathrm{cosh}\;4r+\mathrm{cos}\;2\varphi \;\mathrm{sinh}\;4r},\end{equation} \tag{ E.28 }$$

$$\begin{equation}{{\Gamma}}_{\mathtt{\text{i}}}\left(\varphi \right)=\frac{2\;\mathrm{sin}\;2\varphi \;\mathrm{sinh}\;4r}{\mathrm{cosh}\;4r+\mathrm{cos}\;2\varphi \;\mathrm{sinh}\;4r}.\end{equation} \tag{ E.29 }$$

The operator ϒ(X _m1, X _m2, φ) for the non-interferometric scheme admits a similar polar decomposition as the operator for the interferometric scheme. However, from figure E2 we see that regardless of whether or not optical squeezers are used, there only exists one angle φ that corresponds to maximum entanglement. This occurs at φ = π/2 and corresponds to Bayesian inference of the EPR quadrature. The region between φ = 1.13 and φ = 2.01 is the region where the use of optical squeezers leads to an increase in logarithmic negativity.

Zoom In Zoom Out Reset image size