Generating mechanical and optomechanical entanglement via pulsed interaction and measurement

Entanglement generation at a macroscopic scale offers an exciting avenue to develop new quantum technologies and study fundamental physics on a tabletop. Cavity quantum optomechanics provides an ideal platform to generate and exploit such phenomena owing to the precision of quantum optics combined with recent experimental advances in optomechanical devices. In this work, we propose schemes operating outside the resolved-sideband regime, to prepare and verify both optical-mechanical and mechanical-mechanical entanglement. Our schemes employ pulsed interactions with a duration much less than the mechanical period and, together with homodyne measurements, can both generate and characterize these types of entanglement. To improve the performance of our schemes, a precooling stage comprising prior pulses can be utilized to increase the amount of entanglement prepared, and local optical squeezers may be used to provide resilience against open-system dynamics. The entanglement generated by our schemes is quantified using the logarithmic negativity and is analysed with respect to the strength of the pulsed optomechanical interactions for realistic experimental scenarios including mechanical decoherence and optical loss. Two separate schemes for mechanical entanglement generation are introduced and compared: one scheme based on an optical interferometric design, and the other comprising sequential optomechanical interactions. The pulsed nature of our protocols provides more direct access to these quantum correlations in the time domain, with applications including quantum metrology and tests of quantum decoherence. By considering a parameter set based on recent experiments, the feasibility to generate significant entanglement with our schemes, even with large optical losses, is demonstrated.


Introduction
Entanglement is one of most striking features of quantum mechanics and allows for correlations between two or more objects to be much stronger than is allowed classically. Though this behaviour troubled the founders of quantum mechanics [1,2], entanglement is now very well established, and is viewed as a powerful resource for quantum technologies, such as quantum metrology [3] and communication [4], and for tests of fundamental physics [5].
Entanglement is now also being actively studied both theoretically and experimentally for motional degrees of freedom, opening a rich new avenue for quantum science. Notably, the first demonstration of entanglement between spatially distinct mechanical oscillators has been achieved using trapped ions [17], which reported the observation of non-classical correlations between the motion of two separated pairs of ions. One route to study motional entanglement for more massive systems is via cavity quantum optomechanics [18,19], which utilizes radiation-pressure and other optical forces, such as electrostriction, to generate and study non-classical motional states of mechanical oscillators from the zeptogram to kilogram scale.
Studying the various forms of optomechanical entanglement is an active current area of research and there have been several theoretical, and more recently experimental, studies exploring this avenue. While we do not aim to provide a thorough review here, we briefly describe and contrast some key studies of both optical-mechanical and mechanical-mechanical entanglement.
Focusing first on theoretical proposals, a notable early direction on optomechanical entanglement was to use photon-phonon transfer operations [20], where non-classical and entangled states of light are mapped onto the motion of mechanical oscillators via light-matter beamsplitter interactions [21]. Interactions of this type require operation in the resolved-sideband regime of optomechanics, which is realized if the mechanical frequency far exceeds the cavity decay rate. Following these proposals, detailed studies into the generation and detection of Gaussian optical-mechanical entanglement have been carried out [22,23], which focus on the linearized dynamics around the steady state. In the steady state, optomechanical systems are subject to certain stability conditions, which preclude parts of parameter space, particularly for optical drive on the blue sideband. Long-pulsed optical drives operating in the resolved-sideband regime were suggested as a means to access this part of parameter space [24], and the effect of multiple interactions, pulse shaping, and different optical detunings have also been studied for such long optical pulses [25].
In addition to optical-mechanical entanglement, developing methods to establish entanglement between a pair of mechanical oscillators is also of key interest. In particular, the steady-state linearized dynamics of two mechanical oscillators interacting with the same entangling optical field have been studied to investigate Gaussian entanglement [26]. Alternative experimental configurations have also been proposed, such as the suspension of two mechanical membranes within the same optical cavity [27], and the injection of squeezed light into both a double cavity system [28] and a ring cavity design [29]. The generation of mechanical entanglement via reservoir engineering of a single cavity mode by a multi-tone drive has also been suggested [30]. Recently, a more comprehensive analysis of steadystate Gaussian entanglement has been carried out, which investigates opticalmechanical, mechanical-mechanical, and tripartite light-mechanics entanglement [31]. Importantly, this analysis goes beyond the usual resolved-sideband regime and does not employ the rotating-wave or adiabatic approximations typically used in the literature.
Swapping optical-mechanical to mechanical entanglement by long-pulsed interactions, optical interferometry, and measurements provides another exciting route for generating entangled states of mechanical oscillators. In this setting, linear interactions and measurements on the optical field have been discussed as a way to generate Gaussian continuous-variable entanglement [32]. Furthermore, complementary approaches that implement photoncounting measurements to generate and witness non-Gaussian mechanical entanglement have also been considered [33,34]. Additionally, recent theoretical work, which proposes entanglement swapping between two optical and two mechanical modes offers a promising avenue to use mechanical entanglement to investigate spatially dependent decoherence in massive systems [35].
In recent years, field-mechanics and mechanics-mechanics entanglement experiments have been performed, demonstrating the interest in, and feasibility of, generating optomechanical entanglement. Notably, long-pulsed optomechanical interactions in the resolved-sideband regime have allowed entanglement between microwave fields and mechanical motion to be measured using optical-quadrature measurements [36]. While in the optical domain, photon counting measurements have been used to witness non-classical optomechanical correlations between optical pulses and phonon modes in bulk diamond via off-resonant Raman scattering [37], and then similarly in a nanomechanical resonator operating in the resolved-sideband regime [38]. Furthermore, mechanical entanglement has been established using an interferometric pump-probe scheme, first between two spatially-separated diamonds with photon-counting measurements on Stokes scattered photons [39], and then between mechanical resonators in the resolved-sideband regime of cavity optomechanics [40]. Moreover, continuous driving of a non-interferometric configuration with two micromechanical drum oscillators has allowed for mechanical entanglement to be established mediated by microwave fields [41].
In this work, we propose methods for preparing and verifying both opticalmechanical and mechanical-mechanical entanglement that explore new parameter regimes in the landscape of optomechanical entanglement. Our protocols utilize short optical pulses with a temporal width much less than the mechanical period, such that over the optomechanical interaction both mechanical evolution and decoherence can be ignored. Pulsed optomechanics [42] has allowed for the development of numerous non-classical state preparation and verification schemes, see e.g. Refs [43][44][45][46][47], by utilizing resonant interactions in the unresolved-sideband regime-where the cavity decay rate greatly exceeds the mechanical frequency. This regime of optomechanics also allows the optical and mechanical quadratures to become strongly correlated over a pulsed interaction, which provides a route for generating optomechanical entanglement. We consider pulsed-linearized interactions that enable the generation and full characterization of Gaussian continuous-variable entangled states, which allows for greater entanglement than is possible in low-dimensional discrete variable states. The pulsed nature of our protocols also allows for more direct access to the entanglement in the time domain and permits operation in discrete time steps. Additionally, we compare and contrast interferometric and non-interferometric configurations and analyze the utility of optical squeezing in each case. Here, we see that optical squeezing provides a route to increase the resilience of the entangled state to optical loss and mechanical decoherence.
Experimental progress in optomechanics outside the resolved-sideband regime has advanced to a point where the entanglement schemes we propose are feasible with current systems. For example, mechanical cooling via pulsed-interaction and measurement has been performed [48,49], and developments in sliced-photonic crystal structures now enable large optomechanical coupling rates [50,51]. These experimental advances allow our pulsed optical-mechanical entanglement protocol to be realized with present-day parameters, and our approach to mechanical-mechanical entanglement can be achieved even with total optical efficiencies much less than unity, i.e. of order 10%.
Further work in this direction has a wide range of applications for quantum technologies and fundamental science. To name a few, investigating optomechanical entanglement in the unresolved-sideband regime using short optical pulses provides new routes to: develop quantum sensors [18] and transducers [52], measure decoherence rates in massive systems [35], enable entanglement distillation via optomechanics [53], and investigate non-classical correlations mediated by gravity [54][55][56].
In Section 2, we describe the physical operations that are used in our protocols for optomechanical entangled state preparation and verification. In particular, we describe linearized pulsed-optomechanical interactions, thermal decoherence channels, and optical homodyne measurements. We do so by using the covariance matrix description of quantum states in phase space to efficiently calculate the effect of these Gaussian processes. In Section 3, we use these operations as building blocks to develop a protocol for the preparation and verification of entanglement between light and mechanics. In Section 4, we then build upon these ideas to introduce two protocols that exploit this optical-mechanical entanglement to generate entanglement between two mechanical oscillators.

Pulsed optomechanics and Gaussian operations
For a pulsed optomechanical interaction [42], 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 where g 0 is the intrinsic optomechanical coupling rate [57]. Here, a and b are the optical and mechanical annihilation operators, respectively, and we define the dimensionless optical quadratures as X L = (a + a † )/ √ 2 and P L = −i(a−a † )/ √ 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 where we have taken α ∈ R without loss of generality. In the unresolvedsideband regime, we describe a pulsed-optomechanical interaction with the unitary operator U lin = e iλX M e iχX L X M . The first exponential term in U lin represents a deterministic momentum transfer to the mechanics due to the mean photon number, where λ ∝ α 2 g 0 /κ 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 0 /κ. 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 om = e iχXlXm to describe the optomechanical interaction, as the deterministic part e iλXm 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 [58] . Such states can be fully characterized by the first and second moments of their quadrature vector X = (X l , X m ) T = (X l , P l , X m , P m ) T . Here, the quadrature vectors for each mode are given by X l = (X l , P l ) T and X m = (X m , P m ) 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 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 In this work, we utilize the block matrix form of the covariance matrix: 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∈Sp(2n, R) acting on the quadrature vector X, where the real symplectic group Sp(2n, R) 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 G ρU † G , then in quantum phase space the mapping is described by X → SX, 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.
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 [59,60]. In Hilbert space we represent the action of these phase-insensitive Gaussian decoherence channels on state ρ as E η (ρ) and E γ (ρ), 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 Eq. (5), a Gaussian measurement on the optical mode can be described through the mapping [60] where X meas corresponds to the measurement outcome and σ meas describes the Gaussian measurement [61]. 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.

Precooling and optomechanical entanglement
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 E η , 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 rot (θ) = e iθb † b . During this time, the mechanical mode interacts with a thermal environment with mean occupationN at a rate γ, described by E γ . Including this decoherence, the 8 map that describes the generation of an optical-mechanical entangled state from initial state ρ is therefore Note that the phase insensitive channel commutes with the free evolution of the mechanical mode, i.e.
The mechanical oscillator may be cooled prior to the generation of opticalmechanical entanglement by homodyning the output light following an optomechanical interaction. More specifically, this precooling stage is described by Mechanical precooling, Eq. (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 [42]. 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, Eq. (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 Fig. 1(a).

Entangled state preparation and verification
Our proposed experimental setup for the preparation of an entangled state of light and mechanics is shown in Fig. 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 occupationn. This interaction, followed by optical loss and mechanical decoherence is described by Eq. (8) and generates an entangled state of light and mechanics.
Prior to the entangling stage, the precooling stage, described by Eq. (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  Figure 1: Scheme for pulsed optical-mechanical entanglement preparation and verification. (a) Circuit diagram for the separate stages of the protocol-mechanical precooling, entanglement generation, and verification-with red lines corresponding to the optical modes and the black line representing the mechanical mode. The interaction between the optical pulses and the mechanical state is given by U om = e iχXlXm . Following the precooling stage, the mechanical oscillator freely evolves over a phase-space angle θ = π/2. Then, after the generation of the optical-mechanical entanglement we aim to verify, indicated by the blue figure eight, the mechanics rotates through a general angle θ. Including an optional optical squeezer U sq increases the amount of optomechanical entanglement generated. The jagged arrows indicate optical loss or mechanical decoherence, while the meters represent optical homodyne measurements. Finally, the dashed red lines in the verification stage indicate the different paths the optical pulses may take. (b) The corresponding experimental design. For clarity we don't include the optional optical squeezer. Mechanical precooling is achieved by a direct phase-quadrature measurement of the precooling pulse (P) after the optomechanical interaction, and so the optical switch is set to position 1 during this stage. Verification of the optical and mechanical blocks of the optical-mechanical covariance matrix is achieved by directly homodyning the entangling (E) and verification pulses (V), respectively, and so the optical switch is again set to position 1. To access the correlations block, the entangling and verification pulses interfere on a beamsplitter and the optical output modes are homodyned. To realize this interference, the entangling pulse is sent down a delay line-with the optical switch in position 2-and then the verification pulse is sent towards the beamsplitter by moving the optical switch to position 3. LO refers to the coherent local oscillator pulse entering the homodyne detector. 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 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 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 X v (θ) = (X v (θ), P v (θ)) 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 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 Fig. 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 bs ( π 4 , 0) and S bs ( π 4 , π 2 ). 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 bs ( π 4 , 0) and S bs ( π 4 , π 2 ) operations as (X v (θ), X l (θ)) T and (X v (θ), X l (θ)) T , respectively. Then the C ver block are determined by . (13) 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 [62][63][64][65] and for bipartite Gaussian states [66][67][68]. 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 bỹ where∆ = det A + det B − 2 det C. A necessary and sufficient condition for bipartite entanglement in Gaussian continuous variable systems is given by ν − < 1/2. The logarithmic negativity of the two-mode state with covariance matrix σ is then calculated as 3.2. Optical squeezing and conservative approaches to entanglement estimation 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 Refs [49][50][51].
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 E η , 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. The application of squeezers to protect quantum coherence from interactions with the environment, as discussed in Refs [69][70][71], inspires including this optional squeezing operation in our entanglement protocol. Eq. (8), which describes the generation of optical-mechanical entanglement, is therefore modified to Here, the squeezing unitary is U sq = exp r 2 [(a † ) 2 − a 2 ] , with real squeezing parameter r. Moreover, the channel E ηcav represents losses from the optical cavity and E ηdet is the channel describing detection inefficiencies. We assume that the optical squeezer can be turned off during the precooling Table 1: List of proposed experimental parameters for the optical-mechanical and mechanical entanglement protocols. The parameter values are based on proposed improvements to the sliced-silicon nanobeam system that appears in Refs [49][50][51].
Here, we assume that the mechanical oscillator and its environment are in equilibrium, hencen =N = 500-this corresponds to a temperature of approximately 0.1 K. However, in Ref. [49], a reduction in the thermal occupation fromn = 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. 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 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 superpositon as studied in Refs [69][70][71]. 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 [72,73]. 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 , 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 E η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 Fig. 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.
Importantly, as is evident from the lower plot of Fig. 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  Table 1. Here, the mechanical state has evolved, and decohered, through a phase-space angle of θ = π/2 after the time of entanglement generation. The maximum achievable logarithmic negativity is E N = 2.12, which requires χ = 3.12 and r = 0.60. Without the use of an optical squeezer, a value of E N = 2.01 is obtained with χ = 2.94. The squeezing parameter r sym symmetrizes the optical mode after the optomechanical interaction, while r = r opt maximizes E N (σ) as a function of χ. The contour plot shows r sym > r opt for all χ, which demonstrates that an asymmetric optical phase-space distribution is optimal for maximizing entanglement. Lower: Plot showing logarithmic negativity as a function of interaction strength χ in the presence of both optical loss and mechanical decoherence. The plot also shows that two conservative approaches, which ensure that the logarithmic negativity is not overestimated, can verify a siginificant amount of entanglement, i.e. E N (σ ver ) > 1, for a range of χ. Here, θ refers to the angle through which the mechanical oscillator evolves in phase-space after the optomechanical interaction. The conservative approach in time, which accounts for different times between the entangling and verification pulses, is denoted by (t), while the conservative approach in time and noise, which also accounts for lack of knowledge about the optical noise on the verification pulse, is denoted by (t, P in ). At large χ the two conservative approaches show the same limiting behaviour as the contribution of the optical noise becomes less relevant. The arrow indicates that better approximations of σ can be made from σ ver by applying the inverse of the mechanical decoherence map to each element as noted in the main text and explored further in Appendix B.
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 N (σ) can be reached than in Fig. 2. We refer the reader to Appendix D for a discussion of this case.
The lower plot in Fig. 2 also demonstrates how accurately the opticalmechanical 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 N (σ) ≥ E N (σ ver ), 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 Eq. (13). We refer to this approach as conservative in time, and the conservative nature of this approach is demonstrated in the lower plot of Fig. 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 occupationn 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 Fig. 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.

Mechanical entanglement
Our proposal for preparing an entangled state of two mechanical oscillators builds upon the optical-mechanical entanglement described in the previous section. In this section, we describe schemes for converting entanglement between two optical-mechanical systems to entanglement between two mechanical oscillators using optical homodyne measurements. In particular, we introduce two such complementary schemes. One scheme employs an optical interferometric design and therefore necessitates the introduction of two optical modes, while the other scheme comprises a non-interferometric setup. For each scheme, we compute the amount of entanglement that can be prepared and verified.

Interferometric scheme
The interferometric protocol to entangle two mechanical modes is shown in the circuit diagram of Fig. 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 X l1 = (X l1 , P l1 ) T and X l2 = (X l2 , P l2 ) 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 om ⊗ U om = U ⊗2 om . After these interactions, the decoherence channel E ηcav describing cavity losses acts on the two-mode optical subspace, and optional unitary squeezing operations U ⊗2 sq may be applied to each optical mode to increase the final amount of entanglement generated. Following this, the optical modes experience the loss channel E η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 Fig. 3(b) and discussed in further detail below. During this free evolution the mechanical oscillators decohere via the channel E γ . Including mechanical decoherence, the interferometric entanglement protocol is therefore given by From Eq. (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 (ϕ, ψ) = {(0, π/2), (π/2, 0)}. 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 Fig. 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 1 = θ/ω m and t 2 = φ/ω 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 Eq. (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  Proposed experimental setup to implement the protocol. By virtue of the optical switches and controllable delay lines, the setup allows for mechanical precooling, entangled state preparation, and verification. Each mechanical mode is initially precooled by an interaction with a precooling pulse, followed by an optical phase-quadrature measurement. This is achieved by setting the optical switches to position 1. Mechanical entanglement is generated by first establishing pairwise optical-mechanical entanglement, and then converting this to mechanical entanglement via optical interference and homodyne measurements. Setting the optical switches to position 2 allows for this. The blocks of the covariance matrix which correspond to each mechanical oscillator are reconstructed by interacting each mechanical mode with a verification pulse and directly homodyning the optical outputs, which requires the optical switches to be set to position 1. Instead, if the optical switches are set to position 2, then the mechanicalmechanical correlations can be measured. Here, the other system parameters correspond to those listed in Table 1. Since the time of entanglement generation, each mechanical oscillator has rotated, and decohered, over a quarter of a mechanical period. The maximum logarithmic negativity is E N = 1.07, which corresponds to χ = 3.15 and r = 0.57. Without the use of a squeezer the protocol achieves E N = 0.98 with χ = 2.97. The curve corresponding to the optimal squeezing parameter terminates at the point where, even with squeezing, no entanglement can be generated. We note that the formula for r sym is the same as in the optical-mechanical case. Lower: Logarithmic negativity in the presence of both optical loss and mechanical decoherence after state preparation, and conservative approaches for entanglement verification. The plot shows the case where the two mechanical modes have freely evolved through the same phase-space angles, θ = φ, and hence have decohered by the same amount. Here, r opt is the optimal squeezing parameter to maximize the logarithmic negativity. The arrow indicates that with an inverse transform, better approximations to the entanglement generated at θ = φ = 0 can be made from the statistics gathered using the conservative approaches.
switches are placed after the optomechanical interactions in Fig. 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 Eq. (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. After this operation, the resulting phase quadrature operators at the beamsplitter outputs are given by P v± (θ, φ) = (P v1 (θ) ± P v2 (φ))/ √ 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 .
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 beamsplittercontrollable delay lines may be inserted before and after the optomechanical interactions in the interferometer of Fig. 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. Fig. 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 [74,75]. But as this entanglement measure is also a monotonically decreasing function ofν − , it is an equivalent measure to E N . The lower plot of Fig. 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 Fig. 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.
The interferometric scheme we have introduced above offers an alternative strategy to the mechanical entanglement scheme presented in Ref. [32]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-statecharacterization procedure using these short optical pulses.

Non-interferometric scheme
The non-interferometric scheme for entanglement preparation and verification using pulsed optomechanics is shown in the circuit diagram of Fig. 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 E η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.
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 om = e iχXlXm to U om = e iηcavχXlXm .
After the optomechanical interaction with the second mechanical oscillator, the light passes through the optical loss channels E ηcav and E η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 E γ , is represented by Eq. (19).
Focusing 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 Fig. 5. The A ver and B ver blocks are therefore given by an equivalent expression to Eq. (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 Plots analysing the behaviour of the logarithmic negativity obtained using the non-interferometric scheme are shown in Fig. 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 Here, a single homodyne detector and controllable delay line suffice to realize the precooling, entangling, and verification stages of the protocol. For each mechanical mode, the precooling stage consists of an interaction with the precooling pulse followed by a phase-quadrature measurement on the output pulse. For the first mechanical mode, this is achieved using the switch settings 2 and 3, respectively, while for the second mechanical mode switch settings 1 and 4 are required. Mechanical entanglement is generated by choosing the switch settings 2 and 4 and homodyning the optical output, which corresponds to Bayesian inference of a joint mechanical quadrature. In the verification procedure, the same combination of paths as in the precooling stage will yield the blocks of the covariance matrix corresponding to each mechanical oscillator, while the switch setting 2 and 4 allows the mechanical-mechanical correlations to be accessed.  Table 1. Here, each mechanical mode has decohered over a quarter of a period since the time of entanglement generation. The maximum of the logarithmic negativity is E N = 1.53, which corresponds to χ = 2.97 and r = 0.15. Without the use of an optical squeezer the protocol achieves E N = 1.51 with χ = 3.00. As compared to the previous entanglement schemes, the asymmetry in the optomechanical interactions leads to a different form for r sym , which is given by the real-positive solution of a quartic equation in e 2r . Also in contrast to the previous schemes, the first point in the contour plot with E N (σ) > 0 corresponds to finite value of r opt , which we highlight with the yellow circle. Lower: Logarithmic negativity of the entangled state, showing the effect of mechanical decoherence and optical squeezing on E N (σ), and conservative approaches for mechanical state verification. As the mechanical modes evolve and undergo decoherence, the optical squeezing operations become less useful-demonstrated by the convergence of the orange and red curves. In addition to conservative approaches, the physical requirement that φ ≥ θ leads to further decay of the elements of σ ver .
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 Fig. 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.

Comparison between the entanglement schemes
In this section, we compare the three entanglement schemes outlined above. In particular, we discuss the relative importance of the squeezing operations, the fidelity of the verification procedure, and the sensitivity to the cavity coupling efficiency. Focus will be placed primarily on the interferometric and non-interferometric mechanical entanglement schemes to establish which one is favourable for preparing and verifying entangled mechanical states.
In Fig. 7, we compare the logarithmic negativity generated in each of the three entanglement schemes and also analyze the relative importance of optical squeezing in each scheme. More specifically, in this figure, we plot the percentage increase in the logarithmic negativity due to the use of optical squeezers, and find that squeezing is most useful in the interferometric scheme and least useful in the non-interferometric scheme. It is interesting to note that as the mechanical oscillators evolve and decohere after the instance of entanglement generation, the percentage increase in logarithmic negativity rises for the optical-mechanical and interferometric entanglement schemes, while the percentage increase drops towards 0% for the non-interferometric scheme. This is due to the active role the second optomechanical interaction has on the optical phase-space distribution in the non-interferometric scheme. For the parameter values listed in Table 1, we establish that a squeezing parameter of r ≈ 0.60 is required to reach maximum logarithmic negativity for the optical-mechanical and interferometric entanglement protocols, which is experimentally feasible. However, it is, of course, experimentally simpler to implement the entanglement protocol without the use of optical squeezers, in which case the non-interferometric scheme is advantageous.
In the conservative approaches for state verification, the C ver elements in the non-inteferometric protocol must be measured at later times than those in the interferometric protocol. Despite this requirement, the logarithmic negativity of the verified state E N (σ ver ) is much higher in the non- Figure 7: Comparison of entanglement schemes with and without optical squeezers. The maximum logarithmic negativity E N (σ) generated using the three different entanglement schemes: optical-mechanical, interferometric, and non-interferometric. The squeezing parameter is taken to either be r = 0 or the value which maximizes logarithmic negativity r = r opt . The phase-space angle θ over which the mechanical modes rotate and decohere, after the time of entanglement generation, is taken to be 0, π/2, or 2π. For the mechanical entanglement schemes, we take the second mechanical phase-space angle to be φ = θ. The percentage increase in E N (σ), as a result of optimal squeezing, is listed above each bar.
interferometric scheme than in the interferometric scheme-as is evident from the lower plots of Figs. 4 and 6. Specifically, in the interferometric scheme, the conservative in time approach reaches a maximum logarithmic negativity of E N (σ ver ) = 0.90, as compared to the non-interferometric scheme which reaches E N (σ ver ) = 1.34.
To aid experimental development, here we now investigate the behaviour of the logarithmic negativity as the cavity coupling efficiency is varied, keeping all other parameters in Table 1 the same. Furthermore, we also set all the squeezing parameters to r = 0 for this analysis and assume that the input-cavity-coupling and output-cavity-coupling efficiencies are both given by η cav . We are therefore interested in finding the minimum value of η cav required to satisfy the condition E N (σ) > 0 in each of the three schemes. Table 2 lists this minimum value of η cav for different values of the phase-space angle through which the mechanical modes evolve and decohere after entanglement generation. We find that the optical-mechanical system can still be entangled, even after the mechanical mode has decohered over a full period, for η cav > 0.0065. This entanglement scheme is therefore currently experimentally feasible. The mechanical entanglement schemes lead to stricter requirements on the cavity coupling efficiency, with the non-interferometric scheme being favourable in this regard. However, if one calculates the total optical efficiency in each of the two schemes, we find that only relatively low efficiencies are required to generate entanglement. More specifically, by using Table 2 and Eqs. (17) and (19), one can calculate that in order to generate entanglement total optical efficiencies of 50% and 11% are required for the interferometric and non-interferometric schemes, respectively. Here, the total optical efficiency in the interferometric design is given by the efficiency of a single optical path, i.e. the product of the cavity coupling efficiency and the detection efficiency, while for the linear non-interferometric design the total efficiency is given by product of all the intensity efficiencies along the optical path. Moreover, Table 2 also shows the minimum cavity coupling efficiency required to verify entanglement in each of the schemes, i.e. to satisfy E N (σ ver ) > 0 in the case of no squeezing. These values correspond to minimum total optical efficiencies of 0.63%, 60%, and 16% to verify entanglement in the optical-mechanical, interferometric, and non-interferometric schemes, respectively. Table 2: Minimum cavity coupling efficiency for entanglement generation and verification. The table shows the minimum cavity coupling efficiency η cav required to retain an entangled state after each mechanical oscillator rotates through a phase-space angle θ whilst undergoing some decoherence, and the minimum cavity coupling efficiency required to verify the entangled state using the conservative in time approach. Other systems parameters correspond to those given in Table 1, apart from the optical squeezing parameter r-which we set to zero here. The total optical efficiency is given by the product of the intensity efficiencies along the optical path. To summarise this comparison between schemes, in the parameter range we consider, the non-interferometric design allows for both the preparation and verification of a higher amount of mechanical entanglement-both with and without the use of optical squeezers. The non-interferometric scheme also provides a route for generating entangled mechanical states if the cavity coupling efficiency is the major source of technical difficulty. In fact, the only regime we have studied in which the interferometric design is preferable, is when the output coupling efficiency from the cavities is set to unity, but the input coupling efficiency remains less than unity. In this case, the nonmonotonic behaviour in the logarithmic negativity with increasing interac-tion strength χ is removed as the optical squeezers act before any decoherence channels, and so the entangled state may then be made resilient to optical loss and mechanical decoherence-see Appendix D for a detailed discussion. However, if both input and output coupling efficiencies are set to unity, then the non-interferometric scheme allows for higher values of logarithmic negativity to be reached, and also allows the non-monotonic behaviour to be avoided.

Conclusion
In this work, we have introduced schemes for the preparation and verification of optical-mechanical and mechanical-mechanical entanglement in the unresolved-sideband regime of optomechanics. These schemes utilize short pulsed interactions and optical homodyne measurements to generate these types of entanglement, and by using optical squeezers and an additional precooling stage, we see an increase in the entanglement generated in the presence of both mechanical decoherence and optical loss. Moreover, we have also carefully studied the problem of state verification and have designed procedures for accessing each element of the covariance matrices, while respecting the requirement that each element must be measured at a specific time.
Two mechanical entanglement schemes are introduced and compared using a feasible set of parameters based on current state-of-the-art experiments. We concluded that the non-interferometric design is favourable to the interferometric design in several ways. Whilst both schemes are capable of generating a large amount of entanglement-i.e. E N > 1, the non-interferometric scheme is able to prepare and verify more entanglement in the parameter regime we consider. Furthermore, mechanical entanglement can be prepared and verified at a minimum total optical efficiency of 16% using the noninterferometric approach even without squeezing, while the interferometric scheme requires a higher optical efficiency and squeezing parameter. However, we see that preparing and verifying optical-mechanical entanglement can even be achieved with optical efficiencies much less than 1% and so can be realized with present-day experiments.
This work provides a promising path for the development of new quantum technologies, entangled resources for quantum-information-based tasks, and for the exploration of the foundations of physics using the tools of quantum optics.

31
Appendix A. Symplectic transformations, Gaussian operations, and entangled state preparation in phase space Appendix A.1. Symplectic transformations 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 SX. Here, S belongs to the real symplectic group-S ∈ Sp(2n, R). Therefore, under the action of a Gaussian unitary operation, the first moments and the covariance matrix of a state transform as follows: A necessary and sufficient condition for the covariance matrix σ to represent a Gaussian state is given by the uncertainty relation σ + iΩ/2 ≥ 0 [76].
As the symplectic matrix S preserves the symplectic form SΩS T = Ω, the symplectic transform preserves the Gaussian character of the quantum state.
Appendix A.1.1. Single-mode transformation Consider a single mode described by the quadrature vector X = (X, P ) T , with X = (a + a † )/ √ 2 and P = −i(a − a † )/ √ 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 S rot (θ) = cos θ sin θ − sin θ cos θ .
Squeezer -A squeezing operation is described by the unitary operator . For r ∈ R, the corresponding symplectic matrix is given by Appendix A.1.2. Two-mode transformations A two-mode quadrature vector may be written as X = (X 1 , P 1 , X 2 , P 2 ) T . Here, Beamsplitter -The unitary operation U bs = exp{α cos β(a † 1 a 2 − a 1 a † 2 ) + iα sin β(a † 1 a 2 + a 1 a † 2 )} acts as a general beamsplitter operation and corresponds to the symplectic transformation (A.5) Linearized-pulsed optomechanical interaction-The optomechanical interaction described by e iχX 1 X 2 corresponds to the symplectic transformation

. Decoherence
A phase-insensitive Gaussian decoherence channel acts on the quantum state ρ according to E g (ρ). The channel E g (ρ) 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 [60] This map allows one to model both optical and mechanical decoherence in the entanglement protocols. The description of the optical loss channel E η may be achieved by setting G = η1 2 and σ env to the vacuum optical state. Here, η is the intensity efficiency of the loss channel. Describing the mechanical decoherence channel E γ is accomplished by letting G = e −γt 1 2 and σ env be a thermal state with mean occupationN . Here, the decoherence channel acts on the mechanical modes for a time t and γ is the intrinsic mechanical decay rate.

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 Eq. (5), a Gaussian measurement on the optical mode can be described through the mapping Here, X meas corresponds to the measurement outcome and σ meas describes the Gaussian measurement [61]. In the limit of a projection onto an infinitely squeezed Gaussian state, which constitutes an optical homodyne measurement, the term (A + σ meas ) −1 becomes (Π ϕ AΠ ϕ ) MP . Here, MP refers to the Moore-Penrose pseudo-inverse and Π ϕ is the projection operator onto the X l (ϕ) quadrature Π ϕ = cos 2 ϕ cos ϕ sin ϕ cos ϕ sin ϕ sin 2 ϕ . (A.11) The Krauss operator corresponding to Π ϕ is given by the Hilbert-space pro- Appendix A.2. Precooling stage Eq. (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 wheren is the initial mechanical occupation, then the precooled state of the mechanics is given by Here, 2V x = 1 + 2N + 2e −γt (n −N ) + e −γt χ 2 (A.14) 2V p = (1 + 2N )(1 − e −γt ) + e −γt 1 + 2n 1 + ηχ 2 (1 + 2n) .
The preparation of an entangled state of light and a precooled mechanical oscillator, followed by optical and mechanical decoherence, is described by Eq. (16). Taking 1 2 1 2 ⊕ σ cool as the initial covariance matrix of light and mechanics, then the entangled state has the covariance matrix 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 E γ and the unitary operation described by U rot commute, the form of σ lm in a non-rotating frame may be obtained by computing σ cool , the covariance matrix of each light mode interacting with the mechanical modes is 1 2 1 2 , and the homodyne angles are chosen to be (ϕ, ψ) = (0, π/2), the final bipartite mechanical covariance matrix is given by , , e 2r e −γt ηχ 2 1 + η det (e 2r − 1) .
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.
Appendix A.5. Non-interferometric mechanical entanglement scheme Eq. (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 1 2 1 2 , then the final bipartite mechanical covariance matrix is The functions J a , J b , and J c are given by .
For r = 0, these functions simplify to the form Appendix B. Verification procedure and inverse maps Appendix B.1. Outline of the verification procedure The covariance matrix σ(t) of the entangled state after a time t since entanglement generation is given by 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 don't 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 Eq. (12), we have that and by using P v (θ) = P in + χX m (θ) one finds that, in the absence of decoherence, Similarly, from Eq. (13) we have that P v (θ) and X l (θ) are defined implicitly in the main text and are given by P v (θ) = 1 √ 2 (P v (θ) + X l ) and X l (θ) = 1 √ 2 (X l − P v (θ)), respectively. Therefore, we have 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 conservative in time and noise. (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 In the absence of decoherence, σ(0) ≡ σ(t) and so B ver,12 = B 12 (0).

Appendix B.2. Including optical loss in σ ver
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 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: .

41
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 /η 2 .
Assuming that η is a well-known quantity, none of the results in the reconstruction stage are affected.

Appendix B.3. Inverse map
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 Eq. (A.8), for the elements of σ ver defined at a single time, we have 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 or inversely Elements defined at two different times, t ij,1 and t ij,2 , may be written as Here, c nm,k are real coefficients and σ ij (t ij,k ) is given by 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 . (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 σ(0) ij 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 ) . This allows one to work out B 11 ( π 4 ), B 11 ( 3π 4 ),B 22 ( π 4 ), and B 22 ( 3π 4 ) using Eq. (B.13). Eq. (B.14) may then be used to find an expression for B 12 (0). Namely, 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 .
Appendix C. Single mode Gaussian state of fixed purity in a phaseinsensitive decoherence channel Consider the covariance matrix which characterizes a general single mode Gaussian state of fixed purity µ(σ) = 1/ √ 4Detσ = 1/ √ 4ab up to a local rotation. After the state passes through a phase-insensitive decoherence channel E g , with G = η1 2 and σ env corresponding to a thermal state with occupationN , the covariance matrix is given by Here, δ = (1−η)(1+2N ) 2η > 0. Maximizing µ(σ ) with respect to r gives e 2r = b/a and hence the input state σ which maximizes the output purity from the channel is σ = √ ab1 2 . Likewise, the von-Neumann entropy of the output state [77] is given by where ν is the symplectic eigenvalue of σ -given by The output entropy is minimized with e 2r = b/a, which again leads to a symmetric input state.
interferometric entanglement schemes is avoided. This is because the entangled state is protected from loss of information to the environment before it enters the decoherence channels. In this case, as the optical loss channels are now consecutive, and occur after the optical squeezing operations, they may be combined into a single loss channel with efficiency η = η cav η det . Below, we present a brief study of how optical squeezers are used in the opticalmechanical entanglement scheme to avoid non-monotonic behaviour in this case.
The PPT entanglement criterion reads and can be rewritten as Here, and σ refers to the optomechanical state formed by Eq. (16) when the order of the U sq and E ηcav operations is swapped. We may expand Λ in the parameter χ to study the PPT criterion in the limit χ → ∞. For brevity, this expansion is not reported here, but we find that, with an ansatz for the optical squeezing parameter of the form r = ln(1 + cχ 2 ), it is possible to ensure that lim χ→∞ Λ < 0 if V p < f (c,N , η, e γt ). Here, f is a function of the system parameters and represents the minimal cooling requirement to observe entanglement in the large χ limit-hence avoiding non-monotonic behaviour. Maximising the upper bound over all c gives the cooling requirement for entanglement to persist as χ → ∞ The c which maxmizes f is given by For a genuine mechanical covariance matrix we also require V p > 0. When decoherence processes dominate, f becomes negative and it is no longer possible to avoid separability in the large χ limit. We also note that this cooling condition does not depend on the position variance of the mechanical precooled state V x . This is because the optomechanical unitary generates phase-space displacements in an orthogonal direction to the position quadrature.

Appendix E. Krauss operator approach
In this appendix, we derive the Krauss operators that describe the entanglement protocols of Sections 3 and 4. This Krauss operator representation allows one to further understand the dependence of logarithmic negativity on homodyne measurement angles, and the balance between unitary and positive operations in the generation of entanglement.
As this approach will be completely equivalent to the symplectic formalism outlined in the main text, and so as not to obscure the key elements of this discussion, we neglect the effects of optical loss and mechanical decoherence. However, if one wishes, these deleterious effects may be included: each optical or mechanical loss channel introduces another Hilbert space one would need to consider when calculating the Krauss operator, and hence introduces an index to the Krauss operator indicating the number of photons or phonons which leak to the environment. A sum over all possible photon and phonon numbers detected by the environment would then be performed to describe the lossy Krauss operation. In what follows we take the amplitude of the coherent pulse to be real without loss of generality, α = X α / √ 2 ∈ R.
Appendix E.1. Optical-mechanical and the precooling stage The Krauss operator corresponding to the optical-mechanical entanglement protocol is given simply by Υ(X l , X m ) = U sq (r)e iχXlXm . (E.1) The two mode unitary e iχXlXm 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 η.

(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 P ent = exp −e −2r χ 2 (X 2 m1 + X 2 m2 )/4 cosh 2r , U ent = exp −ie −2r χ 2 (1 + 2 tanh 2r)X m1 X m2 /2 . (E.22) Fig. E.8 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 Figure E.8: Dependence of the logarithmic negativity, with and without optical squeezing, as a function of the optical homodyne angles ϕ and ψ in the interferometric scheme. Here, E N is plotted after the mechanical modes have rotated and decohered through a quarter of a period since entanglement generation. The contour plot with no optical squeezing has χ = 2.97, while the plot with optimal optical squeezing has r = 0.57 and χ = 3.15. These values are the same as those obtained from Fig. 4, which maximize the logarithmic negativity. Without optical squeezing, the logarithmic negativity reaches a maximum of E N = 0.98, while with optimal optical squeezing, E N = 1.07 is reached. There are a continuum of maximum in the region ϕ = [0, π) and ψ = [0, π), and we plot an extra π/4 outside this range to demonstrate the oscillatory nature more clearly.
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.