Entanglement Rate for Gaussian Continuous Variable Beams

We derive a general expression that quantifies the total entanglement production rate in continuous variable systems, where a source emits two entangled Gaussian beams with arbitrary correlators.This expression is especially useful for situations where the source emits an arbitrary frequency spectrum,e.g. when cavities are involved. To exemplify its meaning and potential, we apply it to a four-mode optomechanical setup that enables the simultaneous up- and down-conversion of photons from a drive laser into entangled photon pairs. This setup is efficient in that both the drive and the optomechanical up- and down-conversion can be fully resonant.

To enhance the generation of CV-entangled radiation beams, resonant modes can be exploited (e.g. [9,10] in the microwave domain, or [11,12] in the optical domain). An essential feature of such setups is that the output spectrum is peaked, and strong entanglement may persist only for narrow frequency intervals. This raises the general question how to characterize the overall rate of CV entanglement generation.
In this work, we introduce an "entanglement rate" to quantify CV entanglement in such settings. Afterwards, we illustrate this measure by applying it to a four-mode optomechanical setup that allows the fully resonant, and thereby efficient, generation of entanglement.
The entanglement rate. -For any source of CVentangled beams, we can imagine to apply a frequency filter to each of the beams. Subsequently, the entanglement between filtered modes of the two beams can be calculated, e.g. using the logarithmic negativity [4,30] as an entanglement measure. This filtering analysis has been applied to several settings [19,24,28], for a detailed explanation see e.g. [19,20].
However, if the source emits only in a narrow band- width ( Fig. 1), it can be argued that while the entanglement between individual wave packets may be large, the rate for emitting these is small, producing only a limited total amount of entanglement per unit time. We thus face the question of how to quantify the entanglement rate for a source with an arbitrary spectrum. This will allow us to optimize that rate, which, depending on the intended application, may actually be the prime figure of merit.
A simple approach consists in defining the temporal extent τ of the filter function, quantifying the entanglement E, and then calculating E/τ . This was employed recently to analyze the first experiment on spatially separated CV entanglement in superconducting circuits [9].
However, examples show that this does not always cover the full entanglement (especially when several spectral peaks are involved) and there is some arbitrariness in the choice of filter function and the definition of τ . We now turn to introducing a general definition that remedies those problems.
The natural way to make these notions precise is by using the Wannier-basis of wave packets f m,n (t) that live on a regular grid both in time, t n = nτ , and in frequency-space, ω m = mδω, where δω = 2π/τ . The Fourier-transform F m,n (ω) ≡´e iωt f m,n (t)dt/ √ 2π of these wave packets is nonzero only in the interval ω ∈ [(m− 1 2 )δω, (m+ 1 2 )δω], where F m,n (ω) = exp(iωt n )/ √ δω. The F m,n (ω), and therefore also the f m,n (t), form an orthonormal basis. That makes it possible to decompose uniquely the output field into modes defined by these wave packets:â out (t) = m,nâ m,n f m,n (t), wherê a m,n =´f * m,n (t)â out (t) dt annihilates a photon in mode (m, n), and â m,n ,â † m,ñ = δ m,m δ n,ñ . This construction can be applied to each of the entangled beams [31].
We now calculate the logarithmic negativity E between two wave packets of the two output fields, at the same time slot t n and in suitable frequency slots. In a typical application with a stationary-state source producing beams via parametric down-conversion or four-wave mixing, energy conservation dictates that pairs of frequencies ω 1 and ω 2 are entangled only when they obey To prepare for our definition of the entanglement rate, we have to discuss the dependence of E on the filter time τ , especially for the limit τ τ c , where τ is much longer than the physical correlation time τ c of the source (τ −1 c is the width of the spectral peaks). First, this ensures that there will be no correlations between wave packets located at different time slots, such that we capture the full amount of entanglement. Second, we will now explain in this wave packet picture why E tends to a well-defined limit for τ → ∞, which is consistent with direct calculations [19].
Consider enlarging the filter time τ by a factor N , shrinking the filter frequency interval by δω = δω/N . In effect, the new wave packets of size τ = N τ encompass N of the old wave packets. That this coarsegraining keeps the correlators, and thus the entanglement E, unchanged can be understood already from a very simple consideration. Take a suitably normalized sum of N operators,X = N −1/2 N n=1X n and like-wiseŶ = N −1/2 N n=1Ŷ n . Then the correlator between these "averaged" operators will equal the original correlator: XŶ = X nŶn . This holds provided there have been no cross-correlations and X nŶn is independent of n. The same logic applies to our case, where the new modes are properly normalized averages over the old modes:â m ,n = n K m (n − n N )â m,n , with K encoding the overlap between the two basis sets. The detailed structure of K is not important for our argument, but essentially K is nonzero only in a range of size N , for |n − n N | N , and it is normalized: n |K m (n − n N )| 2 = 1. As shown in [31], this ensures a well-defined limit E(τ → ∞).
The entanglement rate for a given frequency slot m may now be defined naturally as E m /τ , the ratio between the entanglement contained in a pair of wave packets at that frequency and the time τ between those wave packets. (Here m would denote the index for beam 1 to which corresponds uniquely an index m 2 for beam 2, as explained above.) As the logarithmic negativity is additive, it makes sense to add up E m for all the small frequency intervals into which the emission of the source has been decomposed.
We thus arrive at the total entanglement rate, which we define as: Here E m has been the entanglement calculated for the frequency slot m, we have exploited that the limit E(ω) ≡ lim τ →∞ E m is well-defined, where ω = mδω is kept fixed, and we have used δω = 2π/τ to convert the sum into an integral. We may call E(ω) the "spectral density of entanglement" [32]. In the limit τ → ∞, the details of the Wannier basis have become unimportant, and E(ω) can now be calculated using any filter. This definition has a direct operational meaning, corresponding to splitting each of the two entangled beams into many frequency-filtered output beams, where the frequency resolution has been chosen fine enough, such that δω 1/τ c . The rate Γ E quantifies the total entanglement per unit time contained in the sum of those streams. Our definition also shows that E(ω) itself may be interpreted as the entanglement rate per frequency interval.
A specific optomechanical example, and its implementation. -We illustrate the features of the entanglement rate in a model describing the effective interaction between two optical modes (â + andâ − ) and a mechanical mode (b):Ĥ Similar Hamiltonians have been studied previously in the context of entanglement generation between light modes [19][20][21][22][23][24][25][26][27][28]. For instance, Wang and Clerk [26] studied intracavity entanglement, while Tian [27] investigated also the stationary output entanglement. Before proceeding, we note how to implement this model using three equidistant optical modes, enhancing the efficiency beyond previous suggestions. The optical mode spacing J is nearly resonant with the vibration frequency Ω, with a frequency mismatch δ = Ω−J (Fig. 2a). A laser drives the center optical mode at ω 0 , with a detuning ∆ = ω L − ω 0 . An optomechanical interaction of the kind g 0â † +â0 (b +b † ), and similarly forâ − , scatters photons up and down, into modesâ + andâ − , while simultaneously destroying (creating) phonons. When a phonon is virtually emitted and re-absorbed, an effective four-wave mixing process is induced, generating a pair of a + andâ − photons out of a pair ofâ 0 photons. Thus, two-mode vacuum squeezing (EPR entanglement) is produced. We assume that the drive is strong and we can replaceâ 0 by the coherent amplitude α = â 0 . This yields the Hamiltonian (2) with g ≡ 2g 0 α, provided we choose frames rotating at ω L ±J for the modesâ ± , and at J forb. Moreover, only nearly resonant terms are kept, which is allowed if J, Ω g, κ, where κ is the optical intensity decay rate.
We use standard input-output theory [38] for our analysis:ȧ Here Γ is the mechanical damping rate and j = ±.
the thermal occupation, and likewise forâ j,in (but without thermal noise). Solving Eqs. (3) and employinĝ a j,out (t) =â j,in (t) + √ κâ j (t), we find the linear relation between output and input fields in terms of a scattering matrix.
Results. -We first address some general features. The system can become unstable (Fig. 3a), both towards mechanical and optical oscillations. The optical stability boundary is approximately given by (∆ + g 2 /2δ)∆ + (κ/2) 2 = 0. Eliminating the mechanical mode (for δ κ, Γ, g, |∆|), we obtain an effective optical model, which yields the photon pair creation rate This diverges at the optical stability boundary.
We now focus on the output entanglement E(ω) and especially the entanglement rate. In contrast to the intracavity entanglement (discussed in [26,27]), we find that E(ω) is not bounded. This is similar to the difference between intra-cavity and output squeezing [39][40][41]. Numerical plots of E(ω) for the Hamiltonian (2) have been shown so far [27] only for the special case ∆ = δ = 0 and asymmetric optomechanical couplings. Entanglement of temporal modes was also discussed for a pulsed scheme [28] in the case δ = 0.
The output entanglement grows significantly at the stability boundary (Fig. 3b), even diverging in the effective optical model. This is typical near an instability but it comes at the price of linewidth narrowing, reducing the entanglement rate. In order to appreciate this, we now discuss the output intensity spectrum, Fig. 4a,b,c. We expect that any mechanical noise contributing to the optical output is deleterious for entanglement, which is confirmed by explicit calculation. Thus, we plot the spectrum in an instructive way, distinguishing optical and mechanical contributions, as obtained from the linear relationâ +,out (ω) = . . .â +,in (ω)+. . . (â † −,in )(ω)+. . .b in (ω). There are typically two peaks, separated by δ = Ω − J and containing primarily optical or mechanical noise, respectively. Near the optical instability (Fig. 4b), the optical peak gets strong and narrow.
When the optical mode spacing matches the mechanical frequency (δ = 0) and the laser is on resonance (∆ = 0), the two peaks merge and have a narrow linewidth set by the mechanical damping rate Γ ( Fig. 4c). It will turn out that the entanglement rate is maximized near this point. This is entirely counterintuitive: One might expect that at δ = 0 mechanical noise is injected into the output beams, destroying entanglement. However, we find that the optical noise can completely overwhelm the mechanical noise for strong driving, when the cooperativity is sufficiently large, C ≡ g 2 /κΓ n th . The output entanglement E(ω) is shown in Fig. 4d,e,f. Typically, E(ω) is maximal at the optical peak near ω = 0. For the important case ∆ = δ = 0, we find an analytical expression, depends on both driving (via C) and temperature. We checked that choosing different optomechanical couplings for the two modesâ ± will not increase E(ω = 0), in contrast to the intra-cavity case [26].
In Fig. 5a,b, we compare the maximal output entanglement E max ≡ max ω E(ω) and the entanglement rate Γ E =´E(ω)dω/2π. While E max becomes large near the optical boundary of stability, the entanglement rate there remains small, due to the narrow bandwidth. This will be a general feature in many similar systems. Rather, Γ E is optimal near δ = ∆ = 0. We found that, as the mechanical damping rate increases, the optimum shifts away from ∆ = δ = 0.
The entanglement rate depends on the full shape of E(ω), in particular the peak width(s). Crucially, those widths are distinct from those in the output spectrum, due to the nonlinear (logarithmic) dependence of E on parameters. For example, in the case ∆ = δ = 0 of greatest interest, the peak width (see Fig. 4f, Fig. 5c) is not set by the small mechanical linewidth Γ, unlike for the output spectrum discussed above. Thus, the values of the entanglement rate Γ E for the parameters explored here are much larger, of the order κ · max ω E(ω). For increasing temperatures, Γ E decreases, but only slowly, indicating robust entanglement in this setup: Γ E ∝ n −1 th , with the prefactor set by the cooperativity C (Fig. 5d).
Conclusions. -We have introduced an entanglement rate as a quantitative measure for the CV entanglement production per unit time in setups involving resonant modes. Moreover, we have studied an optomechanical setup with one mechanical and three optical modes that allows fully resonant production of optical entanglement.
The output entanglement and the entanglement rate are maximized for different parameter choices. The concept introduced here should be useful for analyzing setups in other domains, like cavities with a Kerr medium [11,12] or microwave resonators with nonlinearities [9,10].
Acknowledgement. -We have benefited from discussions with A. Clerk, D. Vitali  Note added: We recently learned of a related work by Wang, Chesi, and Clerk (in preparation), studying other aspects of output entanglement in this system.

Supplementary material
In this supplementary material we provide some more details regarding the calculation in our paper.

Time-frequency Wannier basis and entanglement rate
First, we present some more technical details on our wave packet based discussion of the entanglement rate. As defined in the main text, we start from a complete orthonormal basis of packets that are localized both on a time-grid (t n = nτ ) and a frequency-grid (ω m = mδω), where τ = 2π/δω. In frequency space, these basis functions are defined as if (m − 1/2)δω ≤ ω < (m + 1/2)δω, and F m,n (ω) = 0 otherwise. In time-space, this corresponds to the well-known Wannier-basis of sinc-shaped wave packets, Upon temporal coarse-graining, we combine N old wave packets into one new one, and at the same time the frequency-resolution becomes refined: τ = N τ and δω = δω/N . Thus, the new frequency index m can be thought of as a combination of the old index m and another index l = 0 . . . N − 1 that splits the old frequency interval into N pieces. The interval belonging to m is thus: ω ∈ [(m − 1/2)δω + lδω , (m − 1/2)δω + (l + 1)δω ]. The definition of F m ,n (ω) on this interval reads like the old one, except for the obvious replacements: F m ,n (ω) = e iωt n / √ δω , with t n ≡ n τ = n N τ . Both the old and the new basis are complete.
We now want to obtain the overlap integrals that relate the old basis to the new one. We find: m , n |m, n ≡ˆF * m ,n (ω)F m,n (ω)dω = K m (n − n N ) , with K m (k) = k = 0 : Note that the overlap K obviously depends on the frequency index m only via the refinement index l in m = (m, l), and that K = 0 if we were to calculate the overlap for any m that is not part of the original frequency interval defined by m. The shape of the overlap as a function of the temporal distance n − n N is that of a sinc function with a decay scale set by N . One can confirm that the overlap matrix elements calculated here are normalized, n |K m (n − n N )| 2 = 1.
When we think of the situation with two entangled beams, we imagine each of them is defined by its own fluctuating output field,â (σ) out (t), where σ = 1, 2 denotes the beam. Each of those can be decomposed into the kind of Wannier basis defined here, and we choose the same filter time τ for each of them. Regarding the frequencies, we want to exploit the fact that in a typical steady-state situation (like in parametric down-conversion or four-wave mixing), it is pairs of frequencies that are entangled. Thus, to each frequency ω 1 of beam 1 belongs an entangled frequency ω 2 of beam 2 (with ω total = ω 1 + ω 2 ). To simplify the subsequent notation, we want to shift and revert the frequency scale of beam 2 such that the two mutually entangled frequencies are always both denoted by the index m. In other words, while ω 1 = mδω, we have ω 2 = ω total − mδω. We note that, after going into a rotating frame that makes the Hamiltonian time-independent (as in the main text), we obtain ω total = 0. In addition, it turns out that due to this matching between opposite frequencies, the basis functions for beam 2 have to be changed accordingly, and the basis transformation for beam 2 is effected by K * instead of K.
We now want to confirm (as indicated in the main text), that such a coarse-graining does not change the entanglement E between wave packets, provided we are already in the regime of sufficiently large filter time, τ τ c . To this end, we just have to show that the correlators between modes do not change upon coarse-graining, i.e. we want to show in particular where the modes on the right-hand side are the temporally coarse-grained ones, and ω m lies within the interval defined by m. The precise location of ω m within that interval will not matter, since δω 1/τ c , so the spectrum of the source is already flat on that scale. In addition, we note that, for the steady-state situation we assume here, neither the left-hand-side nor the right-hand-side of (9) actually depend on the time-point n or n , respectively.
Employing the overlap calculated above, we find: Here we have used thatâ m,n2 . In going to the second line, we exploited the fact that different time slots are already uncorrelated (τ τ c ). We also used the normalization of K in the last step.
Regarding other correlators, such as â † m,n , we can say the following: For the model discussed in the main text, they can be confirmed to be zero by explicit calculation. For more general models, where however the Hamiltonian can still be brought to a time-independent form by a proper choice of rotating frame, we find that stationarity dictates that â (1) (ω) â (2) † (ω ) ∝ δ(ω + ω ), and likewise for all other correlators. In evaluating a correlator of amplitudes in the Wannier-basis, like â m,n are zero, and the proof is in analogy to what was shown above. In summary, the entanglement E will go to a well-defined limit as τ → ∞ (τ τ c ). (As a technical aside, in taking this limit, we assume the frequency index m is re-adjusted such that the center frequency of the corresponding slot is held fixed at some given ω, thus arriving at E(ω) in the limit τ → ∞.) The actual calculation of E(ω) is discussed below.

Hamiltonian and implementation with three coupled optomechanical cells
In this section, we give a derivation of the Hamiltonian (2) and discuss its implementation. We consider three equidistant optical modes with a splitting J coupled to the same mechanical modeb with frequency Ω via radiation pressure. One of the optical modes (here called a 0 ) is driven by an external laser at frequency ω L . Such a setup in general can be described by the following Hamiltonian ( = 1), where for brevity we do not display the coupling to the dissipative environment: Here ω ± = ω 0 ± J, and g (0) q ,q represents the generic (Hermitian) optomechanical coupling matrix. We now assume that the frequency mismatch δ ≡ Ω − J is much smaller than the mechanical frequency, i.e., |δ| Ω and that the mechanical sidebands are resolved, i.e., κ Ω, where κ is the optical damping rate. After transforming to the rotating frame with respect toĤ 0 = (ω L + J)â † +â+ + (ω L − J)â † −â− + ω Lâ † 0â 0 + Jb †b , and neglecting rapidly oscillating terms rotating at ±J by a rotating wave approximation (RWA), we find where ∆ = ω L − ω 0 is the laser drive detuning. For a sufficiently strong laser drive, we can linearize the dynamics by replacingâ 0 by a complex number α = −iΛe −iφ −i +κ/2 . Thus we find the Hamiltonian (2) as given in the main text provided that the couplings to the "+" and "-" mode turn out to be equal, i.e. g/2 ≡ g (0) 0,− α. Without loss of generality, we assume that g is real-valued. The RWA made above is valid when |g| Ω. It may be unavoidable that there is a small asymmetry in the two optomechanical coupling strengths, but since the entanglement generation involves both transitions simultaneously this does not make a crucial difference (we have confirmed this for the output entanglement which we discuss). A more detailed account on some further aspects of asymmetric optomechanical coupling strengths in the context of this kind of Hamiltonian may be found in [26,27], where it is pointed out that for the case of intra-cavity entanglement (which is not our concern here) having such an asymmetry can actually be beneficial.
We now turn to deriving the Hamiltonian (11) for a concrete setup consisting of three coupled optomechanical cells. In each of the cells, we assume a standard (local) coupling between photons and phonons. The microscopic Hamiltonian reads ( = 1) where the index l = 1, 2, 3 runs over the three sites. The second term describes photon tunneling between different sites. It can be taken into account by introducing optical normal modes, as defined byâ ± = ( K1â1+K2â3 √ with eigenfrequencies ω ± = ω 0 ± K 2 1 + K 2 2 , ω 0 respectively. In terms of the normal modes, the Hamiltonian can be written aŝ This Hamiltonian takes essentially the same form as the Hamiltonian (11) given above. We assume that δ = Ω − K 2 1 + K 2 2 Ω and κ Ω, transform to a rotating frame withĤ 0 = q=±,0 (ω q − ω 0 )â † qâq + 3 l=1 lb l , and apply a RWA to find After adding an external laser drive for theâ 0 mode, moving into a frame rotating with the drive frequency and applying standard linearization, this Hamiltonian is identical in form to the Hamiltonian (2) given in the main text. The relevant mechanical normal mode is given byb Note that the design is quite flexible in that it also applies to setups with unbalanced hopping rates K 1 = K 2 and that, in principle, a mechanical mode coupled only to either the first or the third site would suffice. Theâ 0 mode may be driven through an additional channel, provided that the decay rate into this is sufficiently small so as to ensure that the entangled photon pairs dominantly decay into another, outcoupling waveguide.
As pointed out in the main text, a similar configuration of levels may be realized in optomechanical membrane-inthe-middle setups. In such setups, the frequencies of the transverse normal optical modes of the cavity depend on the membrane position and tilt angle [35]. These parameters may be tuned so as to create triplet equidistant optical modes. The frequency separation can be comparable to the frequency of a vibrational mode of the membrane, which could be matched to the optical splitting by applying the optical spring effect [13]. In such a configuration, the (linear) optomechanical coupling strengths are set by the slopes of the optical bands [34]. As pointed out above, even if these turn out to be different, that will not impact our scheme in any important way.

Effective optical model
Next, we derive the effective optical model, i.e. the model obtained after integrating out the mechanics. We will discuss how it captures some essential features of the full model. Assuming a large frequency mismatch δ κ, Γ, g, |∆|, and low temperatures, i.e., √ n th + 1 < g/κ, we can adiabatically eliminate the mechanical mode and the mechanical bath. This can be accomplished by a polaron transformationĤ opt = eŜĤ lin e −Ŝ withŜ = . Thus, we find where, in the last step, we explicitly used that |∆| δ. The optical modes are now decoupled from the mechanical mode and we obtain a closed set of quantum Langevin equations for the optical modeṡ where â q,in = 0, â † q,in (t)â q ,in (t ) = 0, â q,in (t)â † q ,in (t ) = δ qq δ(t − t ) with q, q = ±. The output fields are related byâ ±,out (t) = √ κâ ± (t) +â ±,in (t). For notational convenience, we introduce the operator vectors A = (â +â † +â−â † − ) T and A in = (â +,inâ † +,inâ −,inâ † −,in ) T . The quantum Langevin equations can be cast in the following compact formȦ where The system is stable only if all eigenvalues of M have non-positive real parts. The boundary of stability is located where one of the real parts becomes zero, i.e., − κ 2 + −( g 2 2δ + ) = 0, or equivalently when (∆ + g 2 /(2δ))∆ + κ 2 /4 = 0. This corresponds to the dashed line in Fig. 3a. Note that the lower quadrants of the stability diagram are essentially identical due to inversion symmetry with respect to the point δ = ∆ = 0.
We can solve the quantum Langevin equations in Fourier space. We define the Fourier-transformed operators byâ In the frequency domain, the input noise correlation reads â q,in (ω)â † q ,in (−ω ) = 2πδ qq δ(ω + ω ). For convenience, we define the vectors A in/out (ω) = (â +,in/out (ω)â † +,in/out (−ω)â −,in/out (ω)â † −,in/out (−ω) ) T . The solution of the Langevin equation can be written as A out (ω) = S(ω)A in (ω) with the scattering matrix where I is the 4 × 4 identity matrix. The input noise correlations In the effective optical model, photons are only created in pairs. The pair creation rate Γ pairs is equal to the intensity (photons/sec) in any of the two output streams. Due to the choice of normalization in the input-output formalism, this is given by: This result shows that, depending on the sign of ( + g 2 2δ ) , we get an enhanced or decreased photon pair creation rate compared with = 0. In addition, since the denominator vanishes at the boundary of stability, the photon pair creation rate diverges there.

Calculation of E(ω)
Here, we review the definition of the logarithmic negativity and apply it to quantify the entanglement of the filtered optical output fields, first for the effective optical model introduced above, and then for the full model. By applying narrow frequency filters, we select only single-frequency components of each of the optical output fields. By energy conservation, correlations only occur between the ω-component of one of the fields (sayâ +,out ) and the −ω-component of the other (â −,out ). The field operators for these two single-frequency output fields are obtained asâ +,ω (t) =´t −∞ f (t − s)â +,out (s)ds,â −,−ω (t) =´t −∞ g(t − s)â −,out (s)ds . For convenience, we now will choose Lorentzian filter functions f (t) = 2 τ θ(t)e −( 1 τ −iω)t , g(t) = 2 τ θ(t)e −( 1 τ +iω)t with θ(t) the Heaviside step function, as opposed to the Wannier basis used in the main text. In the limit of small bandwidth 1/τ → 0, which is the only one we discuss, however, both reduce to Dirac δ-functions.
We can use the logarithmic negativity [30] to characterize the entanglement for the output light beams [19]. In order to evaluate it, we define the vector u = (x +, . The entanglement is determined by the covariance matrix V with matrix elements V ij = 1 2 u i u j + u j u i , where the operators involved in this product are all taken at equal times. Inserting the stationary solution, we find where n + =ˆ+ ∞ −∞ e −iωt â † +,out (t)â +,out (0) dt + 1/2 = |S 14 (ω)| 2 + 1/2 , The logarithmic negativity is defined as 2 ) 2 + |x| 2 being the smaller symplectic eigenvalue of the partial transpose of matrix V . Choosing ω = 0, we plot E(ω = 0) as a function of /κ in Fig. 3b.
To obtain the spectral density of entanglement in the full model, we solve the following system of quantum Langevin equations, which derive directly from the model Hamiltonian (main text), and where dissipation and fluctuations have been taken into account using the usual input-output formalism: where, in addition, we have b in = 0, b † in (t)b in (t ) = n th δ(t − t ), b in (t)b † in (t ) = (n th + 1)δ(t − t ), with n th = (exp( Ω/(k B T )) − 1) −1 . For the logarithmic negativity, we need to evaluate the same optical correlations as above. The results are analytical, but too complicated to be reported here. Simpler analytic results can be found for δ = ∆ = 0. In this case we find η − (ω, δ, = 0) = 4C(C + n th + 1 2 ) + 1 2 − 2C 1 + 4(C + n th + 1 2 ) 2 with the cooperativity C = g 2 Γκ . The entanglement E(ω, δ, = 0) is only determined by C and n th as given in the main text.
Comparison to two-and three-mode schemes Finally, we compare the efficiency of our four-mode setup (3 optical, 1 vibrational) to two-and three-mode schemes that have been discussed previously. We show that the four-mode setup is more efficient than either of these schemes.