Entanglement-enhanced time-continuous quantum control in optomechanics

The cavity-optomechanical radiation pressure interaction provides the means to create entanglement between a mechanical oscillator and an electromagnetic field interacting with it. Here we show how we can utilize this entanglement within the framework of time-continuous quantum control, in order to engineer the quantum state of the mechanical system. Specifically, we analyze how to prepare a low-entropy mechanical state by (measurement-based) feedback cooling operated in the blue detuned regime, the creation of bipartite mechanical entanglement via time-continuous entanglement swapping, and preparation of a squeezed mechanical state by time-continuous teleportation. The protocols presented here are feasible in optomechanical systems exhibiting a cooperativity larger than 1.


I. INTRODUCTION
Quantum control plays a crucial role in modern quantum experiments across different fields. In optomechanics alone its applications range from feedback cooling of the mechanical motion [1], mechanical squeezing [2,3] and two-mode squeezing [4] to back-action elimination [5,6] with possible applictions in gravitational wave detection. Importantly for quantum information processing and communication, it can also be used to robustly generate entanglement between remote quantum systems, as has been demonstrated recently for spin qubits [7]. At the same time entanglement itself can be an essential component to facilitate control of quantum systems, e. g., as a resource for teleportation [8] when employed as a means for remote state preparation. In optomechanics, pulsed entanglement between a mechanical oscillator and the electromagnetic field [9] has recently been demonstrated in an electromechanical setup [10]; state preparation (and verification) of an arbitrary mechanical quantum state (e. g., a Fock state) is yet to be accomplished (see, however, [11]). Typical quantum control protocols are operated in a time-continuous fashion and often rely on continuous measurements which are capable of tracking the quantum state of the controlled system. The resulting measurement record-and the so-called conditional quantum state inferred from it-is then used as a basis for the applied feedback [12]. Thus, the success of the control protocol critically depends on the precision of the employed measurement. The most essential prerequisite for quantum limited feedback control turns out to be the regime of strong (linearized, thermal) cooperativity. This regime has been witnessed in several experiments in the past few years [13][14][15][16][17][18]. Recently, monitoring of a mechanical oscillator with a measurement strength matching its thermal decoherence rate (equivalent to a cooperativity above 1) and measurement-based feedback * sebastian.hofer@univie.ac.at cooling to an occupation number of several phonons has been demonstrated in [19].
In this article we explore protocols which exploit optomechanical entanglement as a resource for measurement-based time-continuous control of cavityoptomechanical systems. The presented protocols rely on the fact that for a laser drive appropriately (blue) detuned from the cavity resonance, the radiation-pressure interaction generates entanglement between the mechanical oscillator and the cavity output light. We will show that due to this fact it is in principle possible to feedbackcool the mechanical motion to its ground state, even though the optomechanical interaction effects heating of the mechanical mode in this particular regime. This stands in contrast to common feedback-cooling schemes arXiv:1411.1337v2 [quant-ph] 14 Nov 2014 which operate on cavity resonance [1,19]. The scheme presented here is based on standard homodyne detection to monitor a single quadrature of the output light, and feedback by (phase and amplitude) modulation of the driving laser. The feedback signal is calculated by applying linear-quadratic-Gaussian (LQG) control from classical control theory, which attains the optimal cooling performance for the chosen configuration. Beyond cooling, entanglement-based protocols for quantum control can be used to achieve more sophisticated quantum state engineering. We study in detail two optomechanical implementations of time-continuous Bell measurements [20]: Time-continuous teleportation allows for preparation of a mechanical oscillator in a general Gaussian (squeezed) state, while time-continuous entanglementswapping can be used to prepare two remote mechanical systems in an (Einstein-Podolsky-Rosen) entangled state. Both schemes generate dissipative dynamics which drive the mechanical system(s) into the desired stationary state. They are shown to work if the effective measurement strength is on the same order as the mechanical decoherence rate (i. e., for a optomechanical cooperativity of around 1), which is the same condition that holds for ground-state cooling [14,15], observation of back-action noise [13,18], and ponderomotive squeezing [16,17], all of which have been achieved experimentally. Although we here consider optomechanical systems only, the presented methods are very versatile, applicable to different (continuous and discrete) physical systems, and can be extended to describe more complex topologies, such as multiple interferometric measurements and quantum networks [20].
The results on time-continuous teleportation and entanglement swapping have been published in parts in [20]. Here we provide an extended treatment focusing on an optomechanical implementation and present a detailed derivation of the resulting equations.
The manuscript is organized as follows: In Sec. II we summarize and illustrate the central results concerning cooling, mechanical squeezing, and bipartite mechanical entanglement generation. We start this section by discussing the phase diagram of the optomechanical steady state, emphasizing its unique features for blue detuned laser drive. Sec. III presents in detail all technical aspects in the derivation of our results. Some background information about quantum stochastic calculus and LQG control is presented in the appendices A and B.

A. The cavity-optomechanical system
In this article we consider a cavity optomechanical system with a single mechanical mode oscillating at a resonance frequency ω m . The cavity has a resonance frequency ω c and a [full width at half maximum (FWHM)] decay rate κ, and is driven by continuous-wave laser light at a frequency ω l . In a linearized description and in a frame rotating with ω l , the system is described by the effective Hamiltonian [21] where c m and c l are bosonic annihilation operators of the mechanical and the optical mode respectively. ∆ l = ω 0 − ω c is the detuning of the driving laser with respect to the cavity, and g is the optomechanical coupling strength. In writing this Hamiltonian we implicitly assumed that the cavity is driven strongly, such that the radiation-pressure interaction can be linearized around a large classical intracavity amplitude. The coupling strength is then given by g = g 0 [2κP/ ω 0 (κ 2 + ∆ 2 )] 1/2 with the single-photon coupling g 0 and the input laser power P . To work with the linearized description we assume the existence of a unique classical steady state, thus neglecting effects of bistability [21,22]. The linearized radiation-pressure Hamiltonian H om = g(c m + c † m )(c l + c † l ) can be decomposed into two terms: a beam-splitter like interaction g(c m c † l + c † m c l ) which is resonant for ∆ l = −ω m (red detuned laser drive) and effects cooling of the mechanical motion, and a two-mode squeezing term g(c m c l + c † m c † l ) which is dominant for ∆ l = ω m (blue driving) and is responsible for creation of optomechanical correlations and entanglement. For a resonant drive (∆ l = 0) we retain the full Hamiltonian ∝ gx m x l , which is commonly associated with quantum non-demolition (QND) measurements of harmonic oscillators [23,24] and is the interaction typically used for measurement-based feedback control of these systems [25][26][27]. The optical and mechanical quadrature operators we define as which leads to canonical commutation relations [x i , p j ] = iδ ij . In the following it will be convenient to collect them into the vector operator X = (x m , p m , x l , p l ) T .

B. The optomechanical phase diagram
Optomechanical sideband cooling and entanglement creation in steady state have been analyzed in detail in the literature [28][29][30][31]. Both phenomena are captured by the standard optomechanical master equation (MEQ) for the quantum state ρ, given by [28] where L is the so-called Liouville operator. Here γ denotes the (FWHM) width of the mechanical resonance, andn the mean phonon number of the mechanical bath. Optical and mechanical decoherence is described by the Lindblad operators D[c]ρ = cρc † − 1 2 ρc † c − 1 2 c † cρ. As our system is Gaussian, its state is fully characterized by the first and second moments of X, i. e., the mean values X (t) = tr[Xρ(t)] and the symmetric covariance matrix (3) (Throughout this paper we will often omit the explicit time argument for the sake of brevity if no ambiguity exists.) The linear equations of motion of X and Σ are given by (App. B) The 4×4 matrices F and N describe the system's dynamics and noise properties respectively, and are algebraically connected to the Liouvillian in the MEQ (2). Provided the system is stable, it will in the long-term assume a steady state, lim t→∞ ρ(t) = ρ ss , where ρ ss is determined by the condition Lρ ss = 0. The stability of a linear system can be assessed by applying the Routh-Hurwitz (RH) stability criterion [32] which is fulfilled iff all eigenvalues of F have negative real parts. If a stable steady state exists the means X ss :=tr[Xρ ss ] = 0 will vanish, while the steady-state covariance matrix Σ ss is given by the solution to the so-called Lyapunov equation, which is obtained from (4b) by setting its left-hand side to zero, i. e., FΣ ss + Σ ss F T + N = 0. This equation can readily be solved to obtain steady-state properties such as the mean mechanical occupation number n ss m = c † m c m ss or logarithmic negativity E ss N [33,34]. In this section we will mainly be concerned with these steadystate properties of optomechanical systems.
The characteristic features of an optomechanical system's steady state can nicely be illustrated by plotting a phase diagram with respect to the laser detuning ∆ l and the optomechanical coupling g, as depicted in Fig. 2 for an optomechanical system in the resolved sideband regime (κ < ω m ) for a high-Q (Q = ω m /γ) mechanical oscillator. The grey background depicts the regions of instability, given by the corresponding Routh-Hurwitz criterion, where no steady state exists. The first thing to note is that the system is unstable in nearly all the right half-plane, i. e., for blue detuned laser drive, while for red detuning the system becomes unstable only for appreciably high optomechanical coupling. Centered around the first mechanical sideband at ∆ l = −ω m where the beam-splitter part of the optomechanical interaction is resonant, lies the region where n ss m < 1 (dashed purple line) and thus ground-state cooling is possible. Right at the border of stability, for a similar detuning, we find regions of large steady-state entanglement between the intracavity field and the mechanical resonator (colored in turquoise/blue) [31]. On the opposite side of the phase diagram, around the blue mechanical sideband at ∆ l = ω m , we also expect to observe optomechanical entanglement due to the effect of the optomechanical two-mode squeezing dynamics. However, there the formation of a steady state is inhibited by the optomechanical instability which is due to parametric amplification of the amplitude of both the mechanical and the optical mode [21]. The connection of laser cooling, entanglement generation, and the instability region has been analyzed in detail in [31]. Although no steady The turquoise/blue area shows optomechanical entanglement (logarithmic negativity E ss N ), with largest values close to the instability region. The orange (light gray) line encloses the regions where the conditional mean mechanical occupation number c † m cm ss c < 1 for a measurement of the optical phase quadrature, i. e., φ = π/2. The right axis shows the corresponding optomechanical cooperativity, given by C = g 2 /(n + 1)γκ. Lower Plot: Cut through the phase diagram at g = ωm/10, depicting the conditional phonon number c † m cm ss c for LO phases φ = π/2 (orange line) and φ = 0 (green line). The dashed purple line again shows the mean occupation number for the unconditional state for sideband cooling.
state exists for a blue-detuned laser drive, various alternative approaches permit to work with the resonantly enhanced two-mode squeezing dynamics of the optomechanical interaction. Pulsed optomechanical entanglement creation, for example,-which does not require to be operated in a stable regime-has been analyzed in detail in [9], and has been experimentally demonstrated (for an electromechanical system) in [10]. Working with a continuous-wave blue-detuned laser drive on the other hand, is still possible if we employ stabilizing feedback which inhibits the exponential growth of the optomechanical system's quadratures. One possible type of feedback is measurement-based feedback using homodyne detection, which we will consider in the following.
By adding a homodyne detector to our setup [measuring a single light quadrature defined by the local oscillator (LO) angle φ], we can condition the state of the optomechanical system on the resulting photo-current I(t), which leads to a stochastic master equation (SME) in the Itō sense (see [12,35] and Sec. III A) ρ c is the so-called conditional quantum state, which describes our knowledge of the system given a specific measurement record I(t). The effect of conditioning is described by the operator H is thus nonlinear in ρ c , as is expected for a measurement term. I(t) can be expressed as where W is a Wiener increment with dW (t) 2 = dt, and 0 < η < 1 is the efficiency of the detection. Here and in the following we denote by A c (t) = tr[Aρ c (t)] the expectation value with respect to the conditional state.
In contrast to the conditional state ρ c which solves a SME, we will call the solution of a standard MEQ [such as (2)] the unconditional state, which we denote by ρ.
Gaussian conditional quantum states are fully described by the mean vectorX(t) = X c (t) and covariance matrix defined with respect to ρ c . Their equations of motion are given by a linear stochastic differential equation and a (deterministic) matrix Riccati equation respectively, where H describes the homodyne measurement and M is related to the noise properties of the system (see App. B). K(t) is a time-dependent gain factor which depends onΣ(t). For a one-dimensional system [with a twodimensional phase space (x, p)] these equations allow us to give a simple graphic interpretation of the SME (5) in terms of a phase-space description (see Fig. 3): The conditional trajectoryX (blue line) is determined by the measurements I(t) and therefore follows a random walk in phase space. The covariance matrixΣ (turquoise ellipse) on the other hand evolves deterministically, independent of the measurement results. Averaging over all possible phase-space trajectories recovers the broad Gaussian distribution described by the standard MEQ (2) [or equivalently, equations (4)]. For an unstable system (e. g., in the blue detuned regime), the blue line will spiral outwards, leading to a growing unconditional covariance. The conditional covariance matrixΣ, however, may still possess a (finite) steady state. This is due to the fact that the exponential growth is tracked by the conditional mean, with respect to which the covariance matrix is defined. The steady-state conditional covariance ma-trixΣ ss can be found in analogy to Σ ss by setting the left-hand side of equation (9) to zero and by solving the resulting algebraic Riccati equation [12]. Having obtainedΣ ss we can easily evaluate the mean mechanical occupation number of the conditional state, depicted in Fig. 2(a) for a measurement of the phase quadrature of the light field, i. e., φ = π/2. We find a large region where c † m c m ss c < 0.4 for all detunings −ω m ∆ l ω m (orange line). In the region around the red sideband ∆ l ≈ −ω m this effect can mainly be attributed to passive sideband cooling of the mirror, which we discussed above. However, we now also find a region of low occupation on the opposite (blue) sideband at ∆ l ≈ ω m . In this region the reduction of the conditional phonon number, which at the same time means an increase of the purity of the mechanical state, is due to correlations between the mechanical oscillator and the light field. These correlations allow us to extract information about the mechanics from the homodyne measurement. We will see in the next section that in the sideband-resolved regime κ < ω m this effect is strongest for ∆ l ≈ ω m where the two-mode squeezing, entangling term of the optomechanical Hamiltonian is resonant.
To illustrate how the choice of the LO phase influences the conditional mechanical occupation we plot a cut through Fig. 2(a) at a fixed optomechanical coupling g = ω m /20 in Fig. 2(b). If we choose to measure the optical amplitude quadrature we find that on resonance we do not have a reduction of the conditional phonon number. For a detuned laser drive (|∆ l | ω m ) however, we again find regions of c † m c m ss c < 1. This is easily explained by noting that on resonance (∆ l = 0) only the optical phase quadrature couples to the mechanical oscillator, while the amplitude quadrature contains noise only. Measuring the amplitude quadrature therefore does not allow us to make inferences about the mechanical motion. In general there will be an optimal LO angle, depending on all system parameters (especially g, ∆ l , κ) at which we obtain maximal information about the mechanical motion. Thus, homodyne detection at this particular angle yields the minimal conditional occupation. Typicallyespecially in the weak-coupling regime where g < κthe optimal angle corresponds to the optical quadrature which is anti -squeezed by the optomechanical interaction and thus features the best signal-to-noise ratio. We will see in the next section how these features of the optomechanical phase diagram connect to feedback cooling of the mechanical oscillator.

C. Optomechanical Feedback Cooling
Now that we discussed the conditional optomechanical state in detail, the question arises whether it can be realized via feedback, i. e., whether we can prepare the unconditional state of the system such that it resembles the conditional state.
Consider the setup depicted in Fig. 1(a), where the results from a homodyne measurement of the cavity output light are fed back to the optomechanical system in a suitable manner, such that the mechanical system is driven to a low entropy steady state. This situation has been analyzed in [25][26][27]30]. However, the regime discussed for feedback cooling is typically restricted to resonant drive and the bad-cavity regime κ > ω m . In this section we will discuss that feedback cooling can also be effectively operated in the sideband-resolved regime κ < ω m , and even on the blue sideband ∆ l = ω m , which is normally affiliated with heating. Here we will show that we can harness the entanglement created by the optomechanical two-mode squeezing interaction for a measurement-based feedback scheme, which enables us to cool the mechanical motion to its ground state.
Feedback onto the mechanical system can either be effected by direct driving through a piezoelectric device [36], or by modulation of the laser input, as we will assume in the following. This type of optical feedback can be described by adding an additional time-dependent term to the Hamiltonian, where ε(t) = u x (t) + iu p (t) ∈ C is the complex amplitude of the feedback signal, and |ε(t)| 2 accordingly is the incoming photon flux. To choose an appropriate feedback strategy we employ quantum linear quadratic Gaussian (LQG) control [37], which is designed to minimize a quadratic cost function as described in App. B. Applied to feedback cooling the basic working principle is the following: From the measurement results of the homodyne detection we calculate the system's conditional state ρ c (t), whose evolution is described by (5). Based on this state we can then determine the optimal feedback signal ε(t) which minimizes the steady-state mechanical occupation number . This of course means that the final occupation number depends on the conditional state (more specifically on the covariance matrix Σ), and thus on the chosen LO angle for the homodyne detection as we discussed above. A suitable cost function for this problem is given by with h m > 0. Note that h also includes a contribution by the feedback signal ε, which precludes feedback strategies with unrealistically high feedback strength. The parameter h m therefore effects a trade-off between feedback strength and final occupation number n ss m : high values of h m result in low occupation number possibly requiring large |ε| and vice versa. The mean photon flux in the feedback signal can be calculated as described in App. B. For the parameters used in this section we find that on average |ε| 2 is small compared to the overall driving strength in typical experiments. Only in the region of κ → 0-where almost no photons enter the cavity-the required |ε| 2 may increase dramatically. We note that in order for LQG control to work correctly, certain observability and controllability conditions need to be satisfied [12], which is indeed the case for our system.
The final mechanical occupation is found by first calculating the steady-state variances (∆x m ) 2 and (∆p m ) 2 for a closed feedback loop as outlined in App. B. n ss m is then given by n ss Fig. 4 we plot the steady-state occupation numbers of the feedback-cooled mechanical mode against the laser detuning ∆ l , for the bad-cavity regime (upper plot) and the sideband-resolved regime (lower plot), for two different coupling strengths g. For each detuning the homodyne phase φ is chosen such, that the occupation number is minimized. 1 Note that we keep g constant while varying ∆ l (or κ). This means that the driving laser power has to be adjusted accordingly. In the bad-cavity regime κ > ω m we find that driving on resonance is favorable for both values of g. In this case the optimal LO phase is φ = π/2, as discussed in the previous section. This is the usual regime for feedback cooling [1,19], which is inspired by the idea of quantum non-demolition measurements [38], as they are commonly used in gravitational wave detection.
For micro-mechanical systems, however, the sidebandresolved regime κ < ω m is typically more relevant. In this regime the picture changes completely. For weak coupling (g < κ) we find two pronounced dips at both mechanical sidebands (∆ l = ±ω m ), where n ss m is locally minimal and clearly lies below the value on resonance. It is obvious from the figure that cooling works best on the red sideband (∆ l = −ω m ), where we have a cumulative effect from passive sideband and feedback cooling (see also Fig. 5). However, even on the blue sideband (∆ l = ω m )which is commonly associated with heating-we find an appreciable reduction of the mechanical occupation by feedback cooling. As we discussed in the previous section, we can attribute this effect to large optomechanical correlations, which allow for a good read out of the mechanical motion and thus a good feedback performance. If we increase the coupling strength to g = 2.5κ we see a peak appearing around the blue sideband (which we attribute to ponderomotive squeezing of the output fields), pushing the occupation number above the value at ∆ l = 0. For both regimes we plot graphs for two different detection efficiencies η = 1 (lossless detection) and η = 8/10. Clearly, non-unit detection efficiency leads to a noticeable degradation of feedback-cooling performance. Only at the red sideband and in the sideband-resolved regime, where the effect of sideband cooling dominates, the final occupation number is virtually unaffected. Figure 5(a) shows the mechanical occupation for three detunings ∆ l = 0, ±ω m plotted against g. For ∆ l = −ω m we show, additionally to the closed-loop feedback solution (red solid line), the solution for sideband cooling (red dashed line). While for ∆ l = 0 and ∆ l = −ω m the occupation number steadily decreases-in the depicted range-for growing g, for ∆ l = ω m a clear minimum is visible in the weak coupling regime at g ≈ κ/5. This minimum lies well below the value for ∆ l = 0 (but still above the value for the red sideband). This means that there exists a considerably large parameter regime where a detuned operation significantly improves the performance of feedback cooling. Note that all curves rise drastically in the strong-coupling regime, where g/κ 1 (not shown in the plot). In Fig. 5(b) we plot n ss m against cavity linewidth κ for constant coupling g. Again we find that feedback on the sidebands works best in the sidebandresolved regime, while in the bad cavity regime working on resonance yields (slightly) better performance. Again, the occupation number is minimized with respect to the homodyne phase φ at each point in the plot.
In summary we illustrated that feedback cooling in the resolved-sideband regime is a viable option for cooling the mechanical oscillator into its ground state. It turns out that in this regime driving the system on the blue mechanical sideband yields a lower mechanical occupation number than operating on resonance. As an extension of this protocol we will show in the next section, that a similar setup operating at the same working point can be used to remotely prepare a squeezed mechanical state via time-continuous teleportation.

D. Time-continuous optomechanical teleportation
Time-continuous teleportation is facilitated by what we call a time-continuous Bell measurement [20], as depicted in Fig. 1(b): The output field A of an optomechanical system (denoted by S) is mixed with a second field B on a beam-splitter. The resulting fields are then sent to two homodyne detection setups which measure the Einstein-Podolsky-Rosen (EPR) quadratures x + = x a + x b and p − = p a − p b where x a , x b and p a , p b are the amplitude and phase quadratures of the respective fields. The field B is prepared in a pure state of Gaussian squeezed white noise, which we denote by  5. (color online). Steady-state mechanical occupation number n ss m minimized with respect to φ for different driving frequencies ∆ l = 0, ±ωm, corresponding to a laser drive on the mechanical sidebands and on resonance (represented by different colors/gray levels). Solid lines represent feedback cooling, while the dashed line (for ∆ l = −ωm only) corresponds to sideband cooling without feedback. (a) n ss m against coupling g for fixed cavity decay rate κ = ωm/4 (sideband resolved regime). (b) n ss m against κ for g = ωm/10 (weak coupling regime). Other parameters for both (a) and (b): Q = 10 7 ,n = 3.5 · 10 5 , hm = 100, η = 1 |M , where M ∈ C characterizes the squeezing (see App. C); |M | describes the absolute increase/reduction of the anti-/squeezed quadrature, while arg (M ) determines the squeezing angle. Provided the optomechanical systemfield interaction creates entanglement between the mechanical mode and the outgoing light field, the state of B can be teleported to the mechanical oscillator by applying (instantaneous) feedback proportional to the measurement results of the Bell measurement (I ± ). This effectively generates dissipative dynamics which drive the mechanical system into a steady state coinciding with the input light state. In Sec. III B we derive the constitutive equations of motion (the conditional master equation and feedback master equation) for a generic system. In this section we will focus on the optomechanical implementation. Technical details are discussed in Sec. III D 1.
In order to successfully implement continuous teleportation in optomechanical systems we need to appropriately design our measurement setup. To do this we first need a clear picture of the system's dynamics: In the regime g κ ω m and for a blue drive with ∆ l = ω m the optomechanical interaction is H om ≈ g(c m c l + c † m c † l ). Under the weak-coupling condition (g κ) the cavity follows the mechanical mode adiabatically. We will see that in this regime we effectively obtain the required entangling interaction between the mirror and the outgoing field. Moreover, the mechanical oscillator resonantly scatters photons into the lower sideband at ω c = ω 0 −ω m . Spectrally, the photons which are correlated with the mechanical motion are therefore located at this sideband frequency. We thus set up our Bell measurement in the following way: Firstly, we choose the center frequency of the squeezed input light located at the sideband frequency ω c . Secondly, we now use heterodyne detection to measure quadratures at the same frequency.
In Sec. III D 1 we show that after adiabatic elimination of the cavity mode and a rotating-wave approximation, the evolution of the conditional mechanical state ρ (m) c in a rotating frame with ω m (neglecting the mechanical frequency shift by the optical-spring effect, see Sec. III D 1) is described by the SME where we defined γ − = γ(n + 1) + 2g 2 Re(η − ), γ + = γn + 2g 2 Re(η + ) and η ± = [κ/2 + i(∆ l ± ω m )] −1 . The second row describes passive cooling and heating effects via the optomechanical interaction, as has been derived before in the quantum theory of sideband cooling [28,29]. The third row describes the time-continuous Bell measurement, where the squeezing parameter M is encoded in µ = 1−α and ν = 1+α, with α = (N +M )/(N +M * +1) (App. C). The parameter N > 0 is connected to M via |M | 2 = N (N + 1). η is the detection efficiency as before.
The measured photo-currents of the Bell measurement are where dW ± are correlated Wiener increments whose (co-)variances are given by as is shown in Sec. III B. For the choice ∆ l = ω m we have η + = 2/κ and η − = 1/( κ 2 + 2iω m ). Thus I ± approximately correspond to measurements of the mechanical quadratures p m and x m respectively. We model the feedback as instantaneous displacements of the mechanical oscillator in phase space, where the feedback strength is proportional to the heterodyne currents I ± (t). This is described by Hamiltonian terms I ± (t)F ± , where F ± = F † ± are generalized forces. The feedback operators we choose to be F + = − 2g 2 κ η + x m and F − = − 2g 2 κ η + p m , which generate a displacement in p m and x m respectively. The prefactors of F ± (i. e., the feedback gain) we chose such that they match the measurement strength of the Bell detection. The corresponding feedback master equation (in the same rotating frame) can be written aṡ where = [1 + (4ω m /κ) 2 ] −1 quantifies the suppression of the counter-rotating interaction terms (i. e., the optomechanical beam-splitter). This suppression is large ( 1) in the sideband-resolved regime where κ ω m and small ( ≈ 1) in the bad-cavity regime (κ ω m ). The effective Lindblad terms are determined by λ i and J i , which are obtained (see App. D) from the eigenvalue decomposition of the positive matrix . Moving away from the ideal case, the protocol's performance is degraded by mechanical decoherence effects (γn > 0), counter-rotating terms of the optomechanical interaction which are suppressed by < 1, and inefficient detection (η < 1) which leads to imperfect feedback. Figure 6 shows the steady-state squeezing ζ transmitted to the mechanical mode for different parameters plotted against the optomechanical cooperativity C = g 2 /(n + 1)γκ. In the upper plot we assume perfect detection efficiency η = 1 and find that in this case there exists a critical value C crit (N ) = 1/4[ N (N + 1)] determined by the input squeezing N above which the resulting mechanical state is squeezed for any thermal occupation numbern. The lower plot clearly shows that this is no longer true if we assume non-unit detection efficiency η. We find that below a certain critical value η crit (N,n) we can no longer transfer squeezing to the mechanical oscillator, but we rather heat it instead. (This is even true for a zero-temperature mechanical environment, as illustrated in the plot.) In this general case it can be beneficial to chose a modified feedback gain, i. e., use feedback operatorsF ± = σF ± with σ = 1. In the parameter regime we consider however the resulting improvement negligible.

E. Time-continuous optomechanical entanglement swapping
We now replace the squeezed field mode with a second optomechanical cavity, as is depicted in Fig. 1(c). The goal of this scheme is to generate stationary entanglement between the two mechanical subsystems. This is again facilitated by a time-continuous Bell measurementmeasuring the output light of both cavities-plus feedback [20]. The implementation is akin to the teleportation protocol presented above: Both cavities are driven on the blue sideband to resonantly enhance the two-mode squeezing interaction, and their output light is sent to the Bell detection setup which operates at the cavity resonance frequency ω c . Feeding back the Bell detection results I ± corresponding to the x + and p − quadratures of the optical fields to both mechanical systems dissipatively drives them towards an entangled state. There is a slight complication, however. A single Bell measurement only allows us to separately monitor two of the four variables (x m,1 , p m,1 , x m,2 , p m,2 ) needed to describe the quantum state of the mechanical systems. 2 Combined with the fact that we drive the system on the blue side of the cavity resonance (and thus in an unstable regime) this means that we cannot actively stabilize the system and-depending on the driving strength and sideband resolution-no steady state may exist. To compensate for this we extend the setup by two additional heterodyne detectors, measuring x − and p + with outcomes I ∓ . The effective measurement strength of this stabilizing measurements with respect to the Bell measurement is set by the transmissivity υ of the beam-splitter in front of the heterodyne setup (see Fig. 1). Appropriate feedback of all measurement currents I ± , I ± (for simplicity labeled I i , i = 1, . . . , 4, below) to both mechanical systems finally allows us to stabilize them in an entangled state. Note that this setup effectively realizes two simultaneous Bell measurements of the pairs (x + , p − ) and (x − , p + ) with detection efficiencies υ and 1 − υ respectively. In the rest of this section the two optomechanical systems are assumed to be identical and all detectors to have the same quantum detection efficiency η.
In Sec. III D 2 we show that in an adiabatic approximation the conditional state of the two mechanical oscillators ρ (m) in a rotating frame can be described by the SME (setting ∆ l = ω m ) where we set (J 1 , J 2 ) = √ υ c m,+ , ic m,− , (J 3 , J 4 ) = √ 1 − υ(ic m,+ , c m,− ) and c m,± = c m,1 ±c m,2 . The Wiener processes W i are uncorrelated with unit variance, i. e., dW i dW j = δ ij dt, and correspond to the photo-currents 2 In the language of control theory this means that the complete system is not observable (see for example [12]).
The final steady state of this protocol depends on the feedback operators F i = F † i we apply. In ananolgy to the previous section we choose ( , which can realize independent displacements of all mechanical quadratures. This time we introduced an additional gain parameter σ which we can vary in order to optimize the amount of entanglement in the resulting steady state. With these choices the FME for optomechanical entanglement swapping takes the forṁ where the dynamics generated by the feedback is de- . We can now analyze the stability properties of the linear feedback system by evaluating the corresponding Routh-Hurwitz criterion. In the case of no stabilizing feedback (υ = 1) we find that the admissible optomechanical coupling is limited from above by 4g 2 /κ < 1/(1 − ), which only gives an appreciably high limit for values of ≈ 1 and thus in the bad-cavity regime. The stabilization is caused by the counter-rotating beamsplitter terms c m,i c † l,i +H. c. of the optomechanical Hamiltonian, which cool the mechanical systems. This cooling effect, however, diminishes the amount of generated steady-state entanglement. If we switch on the stabilizing feedback and thus choose υ < 1, we can rewrite the RH criterion in the form [3 + (4g 2 /κ) −1 + ] > 4υ > [(1 − ) − (4g 2 /κ) −1 ]/σ (where we assumed < 1). These inequalities are tightest in the limit → 0, g 2 /κ → ∞ where we have 3 > 4υ > 1/σ. For the rest of this section we choose υ = 3/4 which ensures stability of the feedback system for any values of g 2 /κ and -and consequently the sideband-resolution κ/ω m -as long as the feedback gain fulfills σ > 1/3. In the stable regime and for = 0, η = 1 we find a simple analytic expression for the steady-state logarithmic negativity, where we again introduced the cooperativity C = g 2 /(n+ 1)γκ. Generally we can-for each set of parameters (C,n, , υ, η)-maximize the entanglement E N with respect to the feedback gain σ. In Fig. 7 we plot the resulting steady-state values in terms of logarithmic negativity E N and EPR-variance where x φ m,i = c m,i e −iφ + c † m,i e +iφ and p φ m,i = x φ+π/2 m,i are rotated mechanical quadratures. A Gaussian state is entangled iff ∆ EPR < 2 [39,40]. In the first plot we assume a perfect detection efficiency η = 1 and consider different bath occupation numbersn. We again see that there exists a critical cooperativity C crit above which we are able to generate entanglement regardless ofn. From (19) we can deduce the expression C crit (υ, σ) = 1/[3σ(1+ 4υ − 2σ) − (1 + 4υ)]. (As is evident from the plot, the C crit is independent of ). For the parameters used in the plot (taking into account the optimization with respect to σ) we find C crit = 1/2. Again, counter-rotating terms decrease entanglement but are strongly suppressed by the sideband resolution. In Fig. 7(b) we take into account losses and non-unit detection efficiency, η < 1, which drastically reduces the amount of achieved entanglement. As before we find a critical loss value η crit (n, υ) (for the parameters chosen in the plot slightly above 65%) below which entanglement creation is prohibited.

III. DERIVATION OF CONDITIONAL AND FEEDBACK MASTER EQUATIONS
We present here a brief (and informal) derivation of the stochastic master equations (SME) and Markovian feedback master equations (FME) we use throughout the paper. A rigorous and complete account of the quantum stochastic formalism and quantum filtering theory can be found in the literature, e. g., [12,[41][42][43][44][45], and a brief summary of the most important relations is given in App. A.

A. The homodyne master equation
We consider a situation similar to Fig. 1(a), where a system S couples to the one-dimensional electromagnetic field A, which is initially in vacuum and is subject to homodyne detection. We first assume unit detection efficiency, but will discuss the case of inefficient detection at the end of the section. The system-field coupling is mediated by the Hamiltonian 3 where s is a system operator (e. g., a cavity creation or destruction operator), and the light field is described (in an interaction picture at a central frequency ω 0 ) by a(t) = dω a 0 (ω) e −i(ω−ω0)t [46], where a 0 (ω) is the (Schrödinger) annihilation operator of the field mode at ω. In a Markov approximation the field operators are δ-correlated, and fulfill the commutation relations Under this approximation we can introduce the Itō increments dA, dA † , which (formally) obey dA(t) = a(t) dt, etc., and the vacuum multiplication table given in App.
A. The Schrödinger equation for the full system (S + A) initially in the state |φ 0 = |ψ 0 S |vac A can be written in Itō-form as where we used the fact that dA(t)|φ(t) = dA(t)|φ 0 = dA(t)|vac = 0 [35]. A homodyne measurement of an electromagnetic quadrature x = a + a † with a result I x effectively projects the state of the light field onto the eigenstate |I x of x, where x|I x = I x |I x [47]. Projecting (23) onto |I x leads to the linear stochastic Schrödinger equation (SSE) [48] d|ψ with the forward pointing Itō-increment d|ψ c (t) = |ψ c (t + dt) − |ψ c (t) . |ψ c is the unnormalized system state which is conditioned on the homodyne photocurrent I x . I x is δ-correlated, i. e., I x (t)I x (s) = δ(t−s), and its probability distribution is given (for a fixed time t) by Υ t (I x ) = | I x |φ(t + dt) | 2 [12,47]. Using this, one can show that I x can be written as [12,35] where W is a classical Wiener process with dW (t) 2 = dt and E[dW (t)] = 0. Here the conditional expectation value should be read as A c (t) = ψ c (t)|A|ψ c (t) . We can introduce the classical stochastic process X defined by d X(t) = I x (t) dt, which is statistically equivalent to dA(t) + dA(t) † . This is due to non-demolition properties of the measurement operator, see [44]. It obeys d X(t) 2 = dt.
The according equation of motion for the unnormalized conditional stateρ c = |ψ c ψ c | can be deduced from (24), where we used Itō calculus as presented in App. A. Here the Liouvillian L is given by Equation (26) is the quantum analog to the classical Zakai equation [49]. Note that although the Liouvillian is trace-preserving (i. e., tr[Lρ] = 0), the second term in (26) does not possess this property. The equation for the normalized state ρ c (t) =ρ c (t)/tr[ρ c (t)] is then found by noting that where now A c = tr[ρ c (t)A] and we used tr[ρ c (t)] = 1.
Thus we find which is obtained by expanding tr[ρ c (t + dt)] −1 to second order in d X (which leads to a first-order expansion in dt). Using Itō multiplication rules this leads to which is the desired result [12,35]. The (nonlinear) measurement term is given by It is clear from the derivation that under made assumptions the SSE (24) is equivalent to the SME (30). The stochastic master equation is more general, however, as it can accommodate for additional, unobserved decay channels (such as photon losses/inefficient detection or coupling of a mechanical oscillator to a heat bath), as well as for mixed initial states. We can generalize the homodyne master equation in several ways: Above we assumed a measurement of a specific light quadrature x = a + a † . To measure a rotated quadrature x φ = ae −iφ + a † e +iφ we have to make the replacement s → e iφ s, which simply follows from replacing a → e −iφ a. We can as easily obtain the SME corresponding to heterodyne detection at a LO frequency ω lo = ω 0 by substituting s → e i∆ lo t s, where ∆ lo = ω lo − ω 0 . Below we will discuss the situation where we split up the field with a beam-splitter (with a splitting ration η : 1 − η) and perform two simultaneous homodyne measurements on its outputs. The measured modes A , B after the beam-splitter [see Fig. 1(a)] are related to A before the beam-splitter via [50] where A and B are both initially in vacuum and are uncorrelated such that dA(t) dB † (t) = 0, etc. Plugging this relation into (23) and projecting onto the quadratures e −iφ1 a + e iφ1 a † , e −iφ2 b + e iφ2 b † one can repeat above steps and find with uncorrelated Wiener processes W i , i. e., dW i dW j = δ ij . To model photon losses or inefficient photo-detectors, we average over, say, the second measurement process, and thus discard all information obtained from it. Due to the fact that E[dW 2 ] = 0, the equation of motion for the resulting conditional state-which is now conditioned on the measurement results of the first channel only-is obtained by dropping the last term in (32). The beam-splitter transmissivity is then identified with the efficiency of the photo-detection. Formally we can obtain the same result from (30) by replacing s → √ ηs in the measurement term, while keeping the Liouvillian unchanged.
B. Time-continuous teleportation

Conditional master equation
We now extend the derivation of the homodyne SME presented in the previous section to the Bellmeasurement setup depicted in Fig. 1(b). Again, a onedimensional field mode A [described by a(t)] couples to a system S via (21). A is assumed to be in the vacuum state. A second field mode, B [with a field operator b(t)] is prepared in a pure squeezed state, parametrized by M ∈ C, which we simply denote by |M . M describes the degree and angle of squeezing, as described in App. C. The Itō multiplication table for B is thus given by with the condition |M | 2 = N (N + 1). A and B are combined on a balanced beam splitter, whose output is sent to two homodyne detection setups, which are set up such that they measure the EPR operators We call the continuous measurement of the two quadratures x + and p − a time-continuous Bell measurement [20]. To find the corresponding SME, describing the state of S conditioned on measurements of x + and p − , we apply the same reasoning as in the previous section. We start from the Schrödinger equation (23) where µ = 1 − α, ν = 1 + α and dX + (t) = x + (t) dt, dP − (t) = p − (t) dt. Going from the first to the second line we used the fact that a † +b = (x + +ip − )/ √ 2. We emphasize that x + and p − commute and can be measured simultaneously. We can thus directly project equation (33) onto the EPR state |I + I − AB defined by x + |I + I − AB = I + |I + I − AB and p − |I + I − AB = I − |I + I − AB . This yields the linear stochastic Schrödinger equation where d X + (t) = I + (t) dt and d P − (t) = I − (t) dt are again classical processes which possess the same statistical properties as their quantum counterparts. The photocurrents I ± (analogous to the previous section) can be written as with Wiener increments dW ± . Comparison to the output of a single homodyne setup (25) shows that I ± correspond to two simultaneous homodyne measurements with half efficiency. The (co-)variances of dW ± are given in equations (14) and directly follow from the definition of dX + , dP − and the multiplication tables for dA and dB. We now repeat the procedure from the previous section which is now more involved due to fact that we have to deal with two correlated random processes. It is convenient to introduce the complex process dY (t) = µ d X + (t) + iν d P − (t), which obeys dY 2 = 2ζ dt:=−2N/M * dt and | dY | 2 = 2 dt. We can then write the Zakai equation corresponding to (34) as with L defined in (27). To normalize this equation we first calculate tr[ρ c (t + dt)] = 1 + 1/2 s c (t) dY (t) + H. c. , (37) which we use to obtain (by expanding to second order in dY ) tr[ρ c (t + dt)] −1 = 1 − 1/2 s c (t) dY (t) + H. c.
Combining this with (36) we find after some algebra It can easily be checked that this is a equation of the form (B1) and thus a valid Belavkin equation [51].
To conclude this section let us briefly discuss, as a slight variation of above setup, the situation where instead of the squeezed state |M B we use a displaced squeezed state |M, β B = D(β)|M B (see App. C) as an initial state of mode B, and thus as an input state for teleportation. Transforming the Schrödinger equation (23) into a displaced frame with D(β) shows that the structure of the SSE (24) and the SME (39) remains unchanged, if the measurement processes are replaced by appropriately displaced versions, d X + → d X + − 1/2(β + β * ) dt and d P − → d P − + i 1/2(β − β * ) dt. Consequently, the same transformation has to be applied to the currents I ± in equations (35). A more rigorous derivation of the results in this section is presented in [52].

Feedback master equation
We follow [53] to apply Markovian (i. e., instantaneous) feedback proportional to the homodyne currents I ± to the system S. This procedure consists of three steps: Converting the Itō-equation for the conditional state into Stratonovich form [35], adding a feedback term, and converting back to Itō-form in order to average over all possible measurement trajectories and to obtain an unconditional master equation. The feedback we model as generalized forces F ± = F † ± in the form of additional Hamiltonian terms which we write as K ± ρI ± = −i[F ± I ± , ρ].
We start by rewriting the SME (39) in terms of the complex Wiener increment dW y = µ dW + + iν dW − (this time including the detector efficiency η as discussed in Sec. III A), where we defined Gρ = (s − tr[sρ])ρ and G † ρ = (Gρ) † . The Stratonovich form of this equation is given by Here we used the fact that G † G = GG † . Adding feedback terms [ρ c ] fb = 1/2η[K + I + +K − I − ] and converting back to Itō form yields where we chose an ordering KG to get a trace preserving master equation [53]. We can now average over all possible measurement trajectories (i. e., over all classical processes X + , P − , W y ) to obtain an unconditional (deterministic) master equation. A longish calculation leads toρ = Lρ + 1 2 where we used the fact that 1 2 (K ± ) 2 = D[F ± ] and . w i are the (co-)variances of dW + , dW − and are given by (14). Using the identity this can be written in the more familiar forṁ If we consider again the situation where we use a displace squeezed state |M, β B as input, we make the replacements d X + → d X + − 1/2(β + β * ) dt and d P − → d P − + i 1/2(β − β * ) dt in equation (42). This only changes the third line as all products d X 2 , d X dW y , etc. are unaffected. After taking the classical average this yields an additional Hamiltonian term H coh = √ 2[Re(β)F + − Im(β)F − ] which has to be incorporated into L in the FME (43).
Note that equation (44) is not necessarily a Lindblad equation, as the prefactors to the operators D may in general be negative. It can easily be brought into Lindblad form, however, see App. D.
C. Time-continuous entanglement swapping

Conditional master equation
Consider the setup depicted in 1(c): Two systems S 1 and S 2 couple to field modes A and B (described by field operators a and b, both in vacuum) via interaction Hamiltonians analog to (21). A and B are combined on a 50:50 beam-splitter to form the combinations a ± b in the outputs. These outputs are sent to a pair of beamsplitters (with an uneven splitting ratio υ : 1 − υ) and subsequently to a total of four homodyne setups. If we label the modes incident on the homodyne detectors as We now choose the LO phases of the four homodyne setups such that they measure . These four measurements allow us to simultaneously monitor both quadratures of both systems (although with imperfect precision). The measurement of x + and p − , which we choose to have a relative strength υ set by the beam-splitting ratio, realize a Bell measurement as before, while the measurement of the conjugate quadratures x − and p + , with a strength 1 − υ, we will need for stabilization of S 1 and S 2 .
To derive the SME we apply the same logic as before. We start from the Schrödinger equation for the full system (S 1 + S 2 + field modes), with an initial state |φ = |ψ(0) S1S2 |vac field and an effective Hamiltonian H eff = H sys + H We then rewrite this in terms of dX ± and dP ± and project onto eigenstates corresponding to measurement outcomes I ± , I ± . With the definition s ± = s 1 ± s 2 we find where |ψ is unnormalized. As all electromagnetic field modes are assumed to be in vacuum we find that the measurement processes have unit variance, d X ± (t) 2 = d P ± (t) 2 = dt, and that they are mutually uncorrelated, i. e., d X + d X − = d X + d P + = 0, etc. This can be shown by expressing c i in terms of a, b and using Itō rules, where we have to take into account vacuum noise entering through the open ports of the second pair of beamsplitters (not explicitly introduced here). The homodyne currents are given by where the Wiener increments dW ± and dV ± obey a multiplication table corresponding to the one of d X ± and d P ± . Following the derivation from Sec. III A with four uncorrelated homodyne measurements with non-unit efficiency we can derive the corresponding SME with These results can alternatively be derived in a similar spirit but in a more formal way within the framework of quantum networks, as for example presented in [50].

Feedback master equation
In this entanglement swapping scheme all four homodyne currents, I ± (Bell measurement) and I ± (stabilizing measurements), are fed back to both systems. (For convenience we will in the following label the photo-currents by I i , i = 1, . . . , 4, according to the light modes C i they correspond to.) We again describe this by the operations i now act on the combined Hilbert space of S 1 + S 2 . Using the procedure from before it is straightforward to show that the corresponding FME can be written aṡ with Here we assumed that all detectors have the same efficiency η.

Time-continuous teleporation
Here we derive the stochastic and feedback master equations for the optomechanical teleportation setup outlined in Sec. II D, following the lines of Sec. III B with modifications accommodating the optomechanical implementation. The one-dimensional electromagnetic field A couples to the cavity via the linear interaction . As before, A is assumed to be in the vaccuum state, while B is in a pure squeezed state. In this section we refer to several different rotating frames: the frame of the driving laser rotating at ω 0 (which is our standard frame of reference), the squeezing frame which defines the central frequency for the squeezed input light at ω s , and the local oscillator frame at ω lo in reference to which all measurements will be made. We therefore have the relations with the definitions ∆ lo = ω lo − ω 0 , ∆ s = ω s − ω 0 . The squeezed input state is then, in the squeezing frame, defined by the eigenvalue equation b s (t) − αb † s (t) |M B = 0 with α = (N + M )/(N + M * + 1). The Schrödinger equation of the full system in the LO frame can be written as (neglecting for the moment the coupling to the mechanical bath as this can easily added in the end) where |φ is the state describing the complete system with an initial condition |φ 0 = |ψ 0 S |vac A |M B and δ = ∆ lo − ∆ s . If we now choose ∆ s = ∆ lo , i. e., δ = 0, we can rewrite this as d|φ = −iH eff dt + κ/2 (µ dX + + iν dP − ) c l e i∆ lo t |φ , dt, and µ = 1 − α, ν = 1 + α as before. By comparing this to Schrödinger equation (33) we can deduce that the heterodyne Bell measurement at ω lo is described by SME (39) together with the expression for the measurement currents (35) if we set s = c l e i∆ lo t . Thus the master equation together with the output equations where we neglected contributions from higher sidebands, introducing corrections on the order 1/(δt ω eff m ). With the identification s = −i g 2 κ η + c † m the set of equations (58), (59) is equivalent to the generic case discussed before. However, equation (58) additionally contains decoherence terms due to the coupling to the mechanical environment (γnD[c † m ]+γ(n+1)D[c m ]) and due to optomechanical back-action (2g 2 Re(η − )D[c m ]). For the choice F + = − g 2 κη + (c m + c † m ) and F − = i g 2 κη + (c m − c † m ) (where the prefactors are chosen to match the operator s), and after adding the appropriate decoherence terms, the FME for optomechanical teleportation can be written with U (0) = 1, where s is a system operator, and the effective Hamiltonian is given by H eff = H sys − i 1 2 s † s, where H sys is the Hamiltonian describing the evolution of S. A(t) and A(t) † are the bosonic annihilation and creation processes acting on the Fock space of the electromagnetic field. The increments dA, dA † are forwardpointing, dA(t):=A(t+dt)−A(t), and are (formally) connected to the singular field operators a(t) [46] introduced in Sec. III A by dA(t) = a(t) dt. The definition of the increments leads to the property that they commute with U for equal times, i. e., [U (t), dA(t)] = [U (t), dA † (t)] = 0. If we assume that initially the electromagnetic field is in the vacuum state, the increments obey the multiplication rules × dA dA † dt dA 0 dt 0 dA † 0 0 0 dt 0 0 0 More generally two quantum stochastic processes X(t), Y (t) obey the Itō product rule

d[X(t)Y (t)] = [dX(t)]Y (t) + X(t) dY (t) + dX(t) dY (t),
(A2) again with dX(t):=X(t + dt) − X(t), etc. The standard chain rule is modified in a similar way. For a differentiable function f , we have df (X(t)) = f (X(t)) dX(t) + 1 2 f (X(t)) dX(t) 2 , (A3) which in particular leads to f (X(t + dt)) = f (X(t)) + f (X(t)) dX(t) + 1 2 f (X(t)) dX(t) 2 . To convert between the Itō and the Stratonovich formulation we can use the following approach [35]. Consider the Stratonovich stochastic differential equation with linear operators A and B, a Wiener process W (t) with dW (t) 2 = dt, and a initial condition X(0). This equation has the formal solution where T denotes the time ordered product. We can now calculate the Itō increment dX(t) = X(t+dt)−X(t) and find where we expanded the exponential to second order and used Itō rules. All stochastic differential equations in this manuscript are assumed to be in Itō form unless noted otherwise [and denoted by an (S)]. In this section we briefly review the most important equations of quantum LQG control, following closely the presentation in [54]. Consider a Gaussian n-dimensional open quantum system coupling to m vacuum field channels, m ≤ m of which are subject to homodyne detection. (In the remainder of this section, we will assume that all channels are measured and thus m = m. 4 ) The joint evolution of system plus field is then given by (A1). The system's state conditioned on the outcomes of the homodyne measurements is described by the stochastic master equation (or quantum filter) where dW i are Wiener processes with dW i dW j = δ ij dt and the Hamiltonian is at most quadratic in the system's quadratures, which we collect into a column vector X = (X 1 , . . . , X 2n ) T . The canonical commutation relations can then be written as [X i , X j ] = iJ ij , where J is an skew-symmetric real matrix. We can parametrize L = (L 1 , . . . , L m ) T and H as L = ΛX and where R ∈ R n×n is symmetric,R ∈ C n×m and u(t) is a m-dimensional input signal, which will later be used as a control input. We can describe the system in terms of a vector quantum Langevin equation and an output equation [54] dX(t) = FX(t) + Gu(t) dt + dV (t), dY (t) = HX(t) dt+ dA(t) + dA(t) † , with the definitions F = J[R + Im(Λ † Λ)], H = Λ + Λ † , G = J(R +R * ), and dV = iJ(Λ T dA † − Λ † dA), where dA = (dA 1 , . . . , dA m ) T . We assume the field is in the vacuum state ρ vac , such that dA i (t) dA j (t) = δ ij dt. The measurement currents from the homodyne measurements are (formally) given by I(t) = dY (t)/ dt. Using these definitions we can deduce the equations of motion for the conditional mean valuesX = tr[Xρ c ] and symmetric covariance matrixΣ = tr[XX T ρ c ] −XX T .
We find [54] dX(t) = FX(t) + Gu(t) dt Equations (B4) together with (B5) are known as the Kalman-Bucy filter in classical estimation theory [55]. Assuming a stable system [12], the steady-state solution of the conditional covariance matrixΣ can be found by setting the right-hand side of (B4b) to zero, and solving the resulting algebraic Riccati equation. If instead we are interested in the properties of the unconditional state, we can solve the Lyapunov equation obtained from (B4b) by dropping the last term. [The resulting equation can also be obtained from (B3a) by application of Itō calculus.] The goal of LQG control is to control a system in a way that minimizes a quadratic cost function. In this paper we only deal with the asymptotic control problem for t → ∞ as we are interested in the steady state of our systems. We therefore want to find a feedback strategy which minimizes the total cost [56] ∞ 0 E[h(X(t), u(t))] dt, where we introduced E[·] = tr[ρ(0)ρ vac (·)], the expectation value with respect to the initial state ρ(0) of the system and the vacuum state of the field. We choose h to be of the form h(X, u) = X T PX + u T Qu, where P ≥ 0 and Q > 0 are both real, symmetric matrices of appropriate dimensions. Under the assumption of certain stability conditions [12] the optimal feedback signal is given by [54] u(t) = −C(t)X(t), with Ω ss the solution of the algebraic Riccati equation F T Ω ss + Ω ss F + P − Ω ss GQ −1 G T Ω ss = 0.
In Sec. II C we need to calculate the steady-state covariance matrix of a linear system including optimal feedback. This can be achieved by first noting that (due to the separation principle [12]) we can write (B4a) as where d W is a Wiener process with d W (t) d W (t) T = 1 m dt (the so-called innovations process). We also need that [44,57] Re E[(X(t) −X(t))(X(t) −X(t)) T ] =Σ(t), (B12) E[(X(t) −X(t))X(t) T ] = 0, where the first relation follows from the definition ofX andΣ, and the second from the orthogonality principle [44]. We therefore find The steady-state solution of the symmetric covariance matrix of the controlled quantum system is thus given by Finally, we want to estimate the magnitude of the expected feedback signal. We quantify this by E[u T (t)u(t)].
In the steady state we find