Driven optomechanical systems for mechanical entanglement distribution

The interaction between mechanical devices and light has been the focus of very intensive research activity, in the last six years, directed towards the demonstration of quantum control achieved and operated over quasi-macroscopic systems undergoing quite non-trivial dynamics [1]. This has culminated in the very recent experimental demonstrations of ground-state cooling of massive mechanical modes achieved through radiation pressure-induced dynamical back-action [2], the realization of important technological steps directed towards the construction of quantum-limited position sensors and transducers [3], as well as the engineering of nonlinearities strong enough to mimic, in a fully mechanical context, processes and dynamics that were so far typical of quantum optical settings [4]. The perspective of exploiting massive mechanical systems in order to prepare quantum states and enforce non-classical dynamics has spurred a few theoretical proposals for the achievement of entanglement involving both optical and mechanical systems (from here on referred to as optomechanical entanglement for brevity) [5]. Such a scenario has been enlarged to consider schemes for the generation of fully mechanical entangled states [6, 7] and quantum correlations in hybrid systems involving light, mechanical modes and collective atomic ones [8]. Such achievements are endowed with intrinsic fundamental interest as they open up the way to the experimental study of mechanical versions of the long-craved Schrödinger cat and Einstein–Podolski–Rosen state, as well as the investigation of exotic decoherence mechanisms [9] including the quantum gravitational collapse model [10], the continuous spontaneous localization process [11] and gravitationally induced decoherence [12]. Moreover, the considerable flexibility and high potential for hybridization of optomechanical systems leave sufficient space for a significant re-design of the way a device for quantum communication is thought of. Taking inspiration from the (by now well grounded) paradigm of quantum networks for communication made out of local processing nodes, mutually interconnected by quantum channels [13, 14], one can think about a web of remote and non-interacting optomechanical systems that share multiphoton light fields ensuring long-haul connectivity. In a sense, the embryo of such an architecture was already given by Zhang et al [15], where the coupling between a two-mode squeezed vacuum (TMSV) state and two independent optomechanical cavities was explored for mechanical Hamiltonian engineering. This scheme has been pushed further in [7], where an extensive study of the efficiency of entanglement distribution from purely optical TMSV resources to all-mechanical receivers, including the proposal of a mechanism for the measurement of such entanglement, was presented. It should be noted that alternative configurations for the distribution of entanglement in a mechanical network have been proposed, based on optical postselection [16] or the use of atomic information carriers propagating in the close proximity of monolithic microresonators [17].

In this paper, we go significantly beyond the proof-of-principle study reported in [7] and analyse the performance of mechanical entanglement distribution by means of more general multimode optical resources to determine the extent of the validity of such a paradigm. We start by digging deep into the two-mode resource case, which entails the minimal configuration for the entanglement distribution and, as such, the basic building block for any optomechanical network. By working with generally mixed two-mode Gaussian states, we show that the mechanism very much relies on the symmetry of the optical resource's state: small asymmetries in the statistics of the two modes composing the optical entanglement carrier correspond to a spoiled mechanical entanglement distribution efficiency. The distribution mechanism is optimized by the use of TMSV states whose degree of squeezing (and hence entanglement) is determined by the working point of the mechanical receivers. A posteriori, this demonstrates that the study in [7] indeed entails a configuration of optimal entanglement distribution. We then move to the simplest multipartite setting, consisting of a (generally mixed) three-mode driving field guiding the dynamics of three remote and independent optomechanical cavities. We show that, under conditions of fully symmetric multimode driving fields, tripartite inseparability is indeed in order as signalled by the violation of the positivity of the partial transposition criterion [18]. Our work embodies the first quantitative investigation on the influence that the parameters characterizing a multimode driving field have on the performance of mechanical entanglement distribution. Our results allow for the identification of the optimal resources to use at a given optomechanical working point as well as, conversely, the way an optomechanical network should be adjusted in order to fully exploit the entanglement brought about by the driving field.

The remainder of this paper is organized as follows. In section 1, we introduce the model used for the entanglement distribution among mechanical receivers that are part of simple optomechanical devices, and sketch the analytical procedure to determine the stationary-reduced mechanical state resulting from the interaction with the driving field. Section 2 discusses the two-mode case and the corresponding efficiency of the entanglement distribution process. This analysis is then extended to a three-site case in section 3, where symmetric and genuinely tripartite entangled fields are used in order to get a three-mode inseparable mechanical state. Finally, in section 4 we draw our conclusions and discuss a few open points that we aim at addressing in future.

As a final remark of this introduction, it is important to stress that the working points used for our simulations adhere perfectly to the current experimental state of the art, while the sort of optical states used in order to drive our small optomechanical networks are all well within the reach of linear optics experiments aiming at the engineering of continuous-variable states.

We now provide a general description of the mechanical entanglement distribution scheme without limiting the number of sites of the network that we consider. Figure 1 shows a pictorial representation of the experimental configuration at hand, where the optical circuitry needed in order to generate an optical N-mode entangled resource is condensed in an unspecified 'resource generator'. On the other hand, the mechanical network consists of N remote and mutually independent optomechanical resonators. In order to provide a simple explicit example, we consider the case of Fabry–Perot cavities each endowed with a light mechanical end-mirror, which is one of the configurations currently utilized experimentally [1] and has been exploited to provide the first experimental evidence of radiation pressure-based passive cooling [19], and demonstrate strong optomechanical coupling [20]. Throughout the paper, for the sake of simplicity, we assume identical optomechanical nodes. Therefore, each cavity has length L, natural frequency ω_c, photon decay rate κ and is terminated by a mechanically movable mirror of effective mass m, free frequency (under harmonic approximations) ω_m and mechanical quality factor Q. Cavity j = 1, 2, ..., N is driven by both a large-intensity classical field (pump power η and frequency ω₀) and mode M_j of the (weak) optical entangled resource. The drivings populate the cavity field, which in turn is coupled, via radiation pressure at a bare rate χ = ω_c/L [21], to the mechanical mirror. The free energy of the whole optomechanical network thus reads

$$\begin{eqnarray} &&\fl \hat{H}=\sum ^N_{j=1}\hat{H}_{j}=\sum ^N_{j=1} \Bigg[\hbar (\delta -\chi \hat{q}_j)\hat{a}^\dag _{j}\hat{a}_j\nonumber\\ &&+\,\frac{1}{2m}\hat{p}^2_j+\frac{1}{2}m\omega _m^2\hat{q}^2_j+\hbar {E}\hat{y}_j\bigg], \end{eqnarray} \tag{ 1 }$$

where δ = ω_c − ω₀ is the cavity–pump detuning, ${E}=\sqrt{\frac{2\kappa \mu }{\hbar \omega _0}}$ is the amplitude of the displacement undergone by the jth cavity field due to the strong pump, $\hat{a}$ ($\hat{a}^\dag$) is the corresponding annihilation (creation) operator (with $\hat{y}_j=\mathrm{i}(\hat{a}^\dag _j-\hat{a}_j)/\sqrt{2}$ the out-of-phase quadrature of the cavity field) and $\lbrace \hat{q}_j,\hat{p}_j\rbrace$ the position and momentum operator of the jth mechanical oscillator satisfying the canonical commutation relations $[\hat{q}_j,\hat{p}_k]=\mathrm{i}\hbar \delta _{jk} (j,k=1,\ldots ,N)$. The dynamics induced by this Hamiltonian is captured by the quantum Langevin equations of the whole optomechanical system (we introduce the vector of operators for the subsystem j $\hat{O}_j=(\hat{q}_j,\hat{p}_j,\hat{x}_j,\hat{y}_j)$ with $\hat{x}_j=(\hat{a}^\dag _j+\hat{a}_j)/\sqrt{2}$)):

$$\begin{equation} \partial _t\hat{O}_j=\frac{\mathrm{i}}{\hbar }[\hat{H},\hat{O}_j]+{\rm noise\ terms} (j=1,\ldots ,N), \end{equation} \tag{ 2 }$$

which incorporate the mechanical damping at rate γ_m = ω_m/Q, the cavity photon decay and two noise contributions that deserve a more detailed introduction. The first is the incoherent Brownian dynamics induced by the interaction between the jth mechanical mode and the ensemble of phononic modes excited in the substrate onto which the mechanical oscillator is fabricated [22]. In principle, this mechanism gives rise to non-Markovian evolutions, as encompassed by the two-time correlation function of the (zero-mean) Langevin force operator $\hat{b}_j$ describing Brownian motion [23]:

$$\begin{eqnarray} &&\langle \hat{b}_j(t)\hat{b}_j(t^{\prime })\rangle =\frac{\hbar \gamma _m}{2\pi }m\int \omega {\rm e}^{-{\rm i}\omega (t-t^{\prime })}\left[\coth \left(\frac{\hbar \omega }{2k_{\rm B}T}+1\right)\right]{\rm d}\omega.\nonumber \\ \end{eqnarray} \tag{ 3 }$$

Here T is the temperature of the phononic background modes (with ohmic spectral density) with which each mechanical oscillator is assumed to be in thermal equilibrium and k_B is the Boltzmann constant. Evidently, the noise affecting the mechanical mode is, in general, not delta-correlated. However, for large mechanical quality factors, $\langle \hat{b}_j(t)\hat{b}_j(t^{\prime })\rangle$ becomes approximately delta-correlated, thus reinstating a Markovian description of the associated noise. The validity of such an approximation is guaranteed across this paper.

Zoom In Zoom Out Reset image size

The second source of noise is due to the open nature of the cavities, which determines the effect of the weak entangled driving fields on the dynamics of the mechanical oscillators. In fact, in our model, the pumped optomechanical network is affected by an N-mode quantum-correlated environment that we use in order to distribute the mechanical entanglement. The details of the correlations shared by the environmental modes depend on the form of the state that we decide to use as an entangled resource for the distribution mechanism, and will be specified later on in this paper by assigning the appropriate vector of two-time correlation functions ${\bm C}_{\rm in}(t,t^{\prime })$ (with t and t' being two arbitrary instants of time in the system's evolution).

Finally, under the assumption of intense pump fields, we can focus on the dynamics of the quantum fluctuations $\delta \hat{O}_j$ around the mean value of each operator $\hat{O}_j$, which are thus expanded as

$$\begin{equation} \hat{O}_j\simeq \overline{O}_j\,+\,\delta \hat{O}_j (j=1,\ldots ,N). \end{equation} \tag{ 4 }$$

Each $\overline{O}_j$ is determined by the steady-state Langevin equations [22]. It should be stressed that the optomechanical subsystems are mutually independent, so that each light–mechanical mode interaction occurs regardless of similar events taking place at other sites of the network. In this sense, by neglecting the sources of noise, the time-evolution operator of the whole system would read $\hat{U}(t)=\otimes ^N_{j=1}{\rm e}^{-{\rm i}\hat{H}_j t}$, with $\hat{H}_j$ ruling the coupling between the mode M_j and oscillator j, as specified in (1). The quantum correlations shared by the input noise link the various optomechanical subsystems to each other.

While details on specific and interesting instances of such correlated environments will be given in the next sections, here we briefly sketch the way one can track down the evolution of the system. First, the Langevin equations for the quantum fluctuations introduced above can be compactly cast into the vector form $\partial _t[\delta \hat{O}_j(t)]={K}_j\delta \hat{O}_j(t)+\hat{\mathcal N}_j(t)$ with the dynamical matrix

$$\begin{eqnarray} K_j= \left(\begin{array}{@{}cccc@{}}0 & 1/m & 0 & 0\\ -m\omega ^2_m & -\gamma _m & \sqrt{2}\hbar \chi \Re ({c_{\rm ss}}) & \sqrt{2}\hbar \chi \Im ({c_{\rm ss})}\\\ms -\sqrt{2}\chi \Im (c_{\rm ss}) & 0 & -\kappa & \Delta \\\ms \sqrt{2}\chi \Re (c_{\rm ss}) & 0 & -\Delta & -\kappa \end{array}\right)\nonumber \\&&\ms \end{eqnarray} \tag{ 5 }$$

and the vector of noise fluctuation operators $\hat{\mathcal N}_j=(0,\hat{b}_{j}(t),\sqrt{2\kappa }\delta \hat{x}^{\rm in}_j(t),\sqrt{2\kappa }\delta \hat{y}^{\rm in}_j(t))$, which depends on the Brownian-motion Langevin force and the quadrature fluctuation operators of the input mode j = 1, ..., N. We have introduced the steady-state amplitude of the field in each cavity c_ss = E/(κ + iΔ) and the detuning Δ = δ − ℏχ²|c_ss|²/(mω²_m). It should be noted that, notwithstanding the nonlinear character of the light–mechanical mode interaction, the description in terms of quantum fluctuation operators gives back a fully linear model that guarantees the solvability of the problem encompassed by (5) [5, 22, 7].

Here, we will be interested in the steady-state regime achieved at long interaction times after all transients induced by the light–mechanical mode interaction have faded away. We thus drop the time derivative in the Langevin equations and solve the associated time-dependent set of linear algebraic equations in the Fourier domain. In the following, we will concentrate on the properties of the mechanical system only. We thus report the solution of the Langevin equations for the position and momentum operators of the jth mechanical oscillator, which reads

$$\begin{eqnarray} &&\fl \delta \hat{\mathcal P}_j=\nonumber\\ &&\fl $\frac{m \omega [\sqrt{2\kappa } \hbar \chi ({c}_{\rm ss} {\delta \hat{\mathcal C}^{\rm in\dag }_j} (-\Delta +\mathrm{i}\kappa +\omega )+{\delta \hat{\mathcal C}}^{\rm in}_j c^*_{\rm ss} (\Delta +\mathrm{i}\kappa +\omega ))+\mathrm{i}\delta {\hat{\mathcal B}_j} (\Delta ^2+(\kappa -\mathrm{i}\omega )^2)]}{2{|c_{\rm ss}|^2} \hbar \Delta \chi ^2+m(\mathrm{i}{\gamma _{\rm m}} \omega +\omega ^2-\omega _{\rm m}^2) \left(\Delta ^2+(\kappa -{\rm i} \omega )^2\right)},\nonumber\\ \end{eqnarray} \tag{ 6 }$$

$$\begin{eqnarray} &&\fl \delta \hat{\mathcal Q}_j=\nonumber\\ &&\fl $\frac{\mathrm{i}\sqrt{2\kappa } \hbar \chi [{c}_{\rm ss} {\delta \hat{\mathcal C}^{\rm in\dag }_j} (-\Delta +\mathrm{i}\kappa +\omega )+{\delta \hat{\mathcal C}^{\rm in}_j} {c}^*_{\rm ss} (\Delta +\mathrm{i}\kappa +\omega )]-\delta {\hat{\mathcal B}_j}[\Delta ^2+(\kappa -\mathrm{i}\omega )^2] }{2 |{c}_{\rm ss}|^2 \hbar \Delta \chi ^2 +m(\mathrm{i}{\gamma _{\rm m}} \omega +\omega ^2-\omega _{\rm m}^2) (\Delta ^2+(\kappa -\mathrm{i}\omega )^2)}.\nonumber\\ \end{eqnarray} \tag{ 7 }$$

In the equations above, $\delta \hat{\mathcal O}_j$ stands for the Fourier-transformed fluctuation operator $\delta \hat{O}_j$ (the dependence on frequency is understood) and $\delta \hat{x}^{\rm in}_j=(\delta \hat{c}^{\rm in\dag }_j+\delta \hat{c}^{\rm in}_j)/\sqrt{2}$, $\delta \hat{y}^{\rm in}_j=\mathrm{i}(\delta \hat{c}^{\rm in\dag }_j-\delta \hat{c}^{\rm in}_j)/\sqrt{2}$. Quite clearly, the dynamics of the jth mechanical oscillator depends crucially on the input noise fluctuation operators: from (6) and (7) it is evident that any correlation function of the mechanical variables will be determined by the input noise correlators and, in turn, the entanglement that the engineered environment brings about. In the next sections, we specify the set of input correlators and determine the entanglement that is distributed among the elements of the small mechanical network that we are working on. Finally, throughout our simulations we will use the set of parameters for the optomechanical system given in table 1, which match the current experimental state of the art, thus making our proposal very close to feasibility.

Here, we restrict our study to the case of N = 2 sites of the network, which is thus driven by a two-mode quantum-correlated state of light. The linearity of the dynamics of the quantum fluctuations ensures that the corresponding map preserves the Gaussian nature of an input optomechanical state. This observation allows us to limit our attention to just this class of continuous-variable states. In fact, Gaussian states are most easily generated in quantum optics laboratories. Moreover, Gaussian resources appear to maximize the efficiency of schemes for the transfer of entanglement from continuous- to discrete-variable systems, when compared to experimentally accessible non-Gaussian states such as photon-subtracted or photon-added fields [14]. Finally, a key assumption of our model is that each mechanical oscillator is initially in thermal equilibrium with its own background of phononic modes; we take the elements of the mechanical network as initialized in thermal states with mean occupation number $\overline{n}_j=[{\mathrm{e}}^{\frac{\hbar \omega _m}{k_{\rm B}T}}-1]^{-1}$. The Gaussian-preserving nature of the problem allows us to use the formalism of covariance matrices: given an N-mode system described by the 2N long set of quadrature operators $\hat{Q}=\lbrace \hat{q}_1,\hat{p}_1,\ldots , \hat{q}_N,\hat{p}_N\rbrace$, we define the covariance matrix elements $({\bm \sigma })_{kl}=\langle \hat{Q}_k\hat{Q}_l+\hat{Q}_l\hat{Q}_k\rangle /2 (k,l=1,2,\ldots ,2N)$. The covariance matrix of a physical state should satisfy the Heisenberg–Robertson uncertainty principle ${\bm \sigma }+\mathrm{i}\Omega _N/2\ge 0$, with Ω_N = ⊕^N_{j = 1}iσ_y the so-called symplectic matrix and σ_y the y-Pauli matrix. The Gaussian state of the system is determined by the assignment of ${\bm \sigma }$.

We now need a general description of the two-mode resource that we are going to employ. The covariance matrix of any (pure or mixed) two-mode Gaussian state can be put in the so-called standard form parameterized as [24]

$$\begin{equation} {\bm \sigma }=\frac{1}{2} \left(\begin{array}{@{}cccc@{}}a&0&c_+&0\\ 0&a&0&c_-\\ c_+&0&b&0\\ 0&c_-&0&b \end{array}\right) , \end{equation} \tag{ 8 }$$

with

$$\begin{eqnarray} &&\fl a=s+d,\quad b=s-d, \nonumber\\ &&\fl c_\pm =\frac{\sqrt{(f_d-h_d)^2-4g^2}\pm \sqrt{(f_s-h_d)^2-4g^2}}{4\sqrt{s^2-d^2}} \end{eqnarray} \tag{ 9 }$$

and f_x = 4x² + (g² + 1)(λ − 1)/2(x = d, s). The state corresponding to ${\bm \sigma }$ is physical and entangled for s ⩾ 1, |d| ⩽ s − 1, 2|d| + 1 ⩽ g ⩽ 2s − 1 and λ ∈ [ − 1, 1]. Here, s and d set the local entropy of each mode, g determines the global purity of the state, while the degree of entanglement shared by the two modes grows with λ. The assignment of the input optical covariance matrix ${\bm \sigma }_{\rm opt}$ specifies the set of two-time correlators ${\bm C}_{\rm in}(t,t^{\prime })$ needed to solve the mechanical dynamics. Indeed, a covariance matrix such as (8) is obtained by taking an optical state having only non-zero correlators as

$$\begin{eqnarray} &&\fl {\bm C}_{\rm in}(t,t^{\prime })=(\langle \delta \hat{c}^{\rm in\dag }_1(t)\delta \hat{c}^{\rm in}_1(t^{\prime })\rangle , \langle \delta \hat{c}^{\rm in}_1(t)\delta \hat{c}^{\rm in\dag }_1(t^{\prime })\rangle ,\nonumber\\ &&\langle \delta \hat{c}^{\rm in\dag }_2(t)\delta \hat{c}^{\rm in}_2(t^{\prime })\rangle , \end{eqnarray} \tag{ 10 }$$

$$\begin{eqnarray} \fl && \langle \delta \hat{c}^{\rm in}_2(t)\delta \hat{c}^{\rm in\dag }_2(t^{\prime })\rangle ,\langle \delta \hat{c}^{\rm in}_1(t)\delta \hat{c}^{\rm in}_2(t^{\prime })\rangle , \langle \delta \hat{c}^{\rm in\dag }_1(t)\delta \hat{c}^{\rm in\dag }_2(t^{\prime })\rangle ) \end{eqnarray} \tag{ 11 }$$

$$\begin{eqnarray} &&\fl =\left(a-\frac{1}{2},a+\frac{1}{2},b-\frac{1}{2},b+\frac{1}{2},\mathrm{e}^{-\mathrm{i}\omega _m(t-t^{\prime })}\frac{(c_+{-}c_-)}{2},\right.\nonumber\\ &&\left.\mathrm{e}^{\mathrm{i}\omega _m(t-t^{\prime })}\frac{(c_+{+}c_-)}{2}\right)\delta (t-t^{\prime }). \end{eqnarray} \tag{ 12 }$$

Physically, this corresponds to taking a weak driving field whose modes have frequency ω₀ + ω_m [7]. The covariance matrix ${\bm \sigma }_{\rm mech}$ resulting from the assignment of a ${\bm \sigma }_{\rm opt}$ equal to (8) maintains the standard form and reads

$$\begin{eqnarray} {\bm \sigma }_{\rm mech}=\left( \begin{array}{@{}cccc@{}}\langle \delta \hat{Q}_1^2\rangle &0& \langle \delta \hat{Q}_1\delta \hat{Q}_2\rangle &0\\ 0& \langle \delta \hat{P}_1^2\rangle &0& \langle \delta \hat{P}_1\delta \hat{P}_2 \rangle \\ \langle \delta \hat{Q}_1\delta \hat{Q}_2\rangle &0&\langle \delta \hat{Q}_2^2\rangle &0\\ 0&\langle \delta \hat{P}_1\delta \hat{P}_2 \rangle &0&\langle \delta \hat{P}_2^2\rangle \end{array}\right),\nonumber\\ \end{eqnarray} \tag{ 13 }$$

where each element is calculated using (6) and (7) and then transforming back to the time domain [7]. We have generated a random sample of input two-mode states by taking parameters uniformly distributed across their ranges of variation, which are in turn determined by setting an upper bound to the possible values of s. Correspondingly, we have calculated the output entanglement distributed between the mechanical oscillators. This has been quantified using logarithmic negativity [25], which for Gaussian states is defined as ${\mathcal L}_{\mathcal N}=\max [0,-\ln 2\tilde{\nu }_-]$, where $\tilde{\nu }_-=\min {\rm eig}|\mathrm{i}\Omega _2\tilde{\bm \sigma }|$ is the minimum symplectic eigenvalue of the covariance matrix $\tilde{\bm \sigma }=P_{1|2}{\bm \sigma }P_{1|2}$ associated with the partially transposed state of the mechanical network [24] (P_1|2 = diag[1, 1, 1, −1] is the matrix which, by inverting the sign of momentum of the mechanical oscillator 2, realizes partial transposition at the level of covariance matrices [26]). The results are displayed in figure 2(a), where we compare the input optical entanglement (determined calculating the logarithmic negativity of the two-mode input state) and the corresponding output mechanical entanglement associated with (13) and calculated at the detuning Δ = ω_m that maximizes it. Although the range of entanglement of the input states extends uniformly from 0 to more than 3, the mechanical entanglement can only be distributed when the input one is larger than a non-zero threshold determined by the working parameters at the optomechanical level and not the input field's properties. We have also found that the system favours symmetric input states: the number of entangled mechanical states obtained when considering input states characterized by d = 0 is much larger than the analogous number for asymmetric states. This results from the uniformity of the parameters within the optomechanical network and is an effect that can be cured by matching the optical asymmetry with an equal mechanical one. Most importantly, the entanglement distributed in the mechanical network is constrained to stay below a curve that maximizes the degree of output entanglement at assigned values of the input one. We have found that the optical resource giving rise to such boundary should have d = 0 and g = λ = 1, i.e. it should be a pure, symmetric state (thus equivalent to a TMSV) with covariance matrix

$$\begin{eqnarray} {\bm \sigma }_{\rm opt}=\frac{1}{2} \left(\begin{array}{@{}cccc@{}}s&0&\sqrt{s^2-1}&0\\ 0&s&0&-\sqrt{s^2-1}\\ \sqrt{s^2-1}&0&s&0\\ 0&-\sqrt{s^2-1}&0&s \end{array}\right).\nonumber\\ \end{eqnarray} \tag{ 14 }$$

The behaviour of the boundary suggested by the numerical exploration is confirmed by an exact analytical approach which we sketch here. First, we write the input optical entanglement in terms of the parameter s characterizing the optimal resource as

$$\begin{equation} \epsilon =-\case{1}{2}\,\ln |{{1-2{s}^2+2s \sqrt{{s}^2-1}}}|\rightarrow s=\cosh \epsilon . \end{equation} \tag{ 15 }$$

This expression for s can then be replaced in the analytical one for $\tilde{\nu }_-$, the minimum symplectic eigenvalue of the partially transposed covariance matrix of the mechanical network. While it is possible to get an exact closed expression for $\tilde{\nu }_-$ as a function of all the parameters of the optomechanical setting, this is too lengthy to be reported here. In figure 2(b), we show such an analytic bound for the parameters used in panel (a). We defer an analysis of the behaviour of such an expression against some of the key parameters of our system to later.

Zoom In Zoom Out Reset image size

Therefore, a TMSV embodies the ideal resource for the distribution of entanglement in the mechanical network, thus proving that the analysis in [7] is optimal. Moreover, we can determine the value of s that, at the optimal detuning Δ = ω_m, maximizes the transferred entanglement, thus identifying the optimal input resource for the distribution scheme. To do this, we look for the solution of the equation $\partial _s\tilde{\nu }_-(s)=0$ at g = λ = 1 and for a symmetric resource. Figure 2(c) shows the trend followed by both $\tilde{\nu }_-$ and its derivative with respect to s. At s = 2.501, the symplectic eigenvalue of the partially transposed mechanical state achieves its minimum: again, the analytic expression for such a solution is too uninformative to be shown here.

While, in line with expectations, we have checked that increasing mass and/or mechanical damping rate reduce the amount of entanglement between the mechanical systems; the effects of the deviations from ideality of the resource are yet to be assessed. This is relevant given that it is likely that the optical resource produced in a quantum laboratory would be mixed (due, for instance, to residual thermal nature at the two-mode state generation stage). We thus now address the effects that a non-unit purity of the two-mode optical resource has on the performance of our scheme. In order to do this, we have considered the range of values of (s, g) for a symmetric optical input state (with λ = 1) that leaves $\tilde{\nu }_-<1/2$. As shown in figure 3(a), the scheme tolerates a considerable degree of mixedness of the input resource: only at g = 5.01 (corresponding to a state purity of 0.2) and for the optomechanical parameters stated in table 1, does the distributed entanglement disappear. The value of s that optimizes the entanglement distribution increases with g, making the function $\tilde{\nu }_-(g,s)$ asymmetric with respect to s. This is seen in figure 3(b), where we highlight the region of the (g, s) space within which the mechanical entanglement is found.

Zoom In Zoom Out Reset image size

Another important parameter to consider is the initial temperature T of the mechanical systems. How warm should the network be in order to jeopardize the entanglement distribution mechanism? The process appears to be robust against thermal effects. Figure 4(a) shows that the decrease of the maximum value of the output mechanical entanglement is quasi-exponential and very flat at small temperatures. (The drop in the maximum entanglement from T = 10⁻⁶ K to 10⁻⁴ K is only quite marginal. Note that the horizontal axis is in logarithmic scale.) The entanglement persists up to T = 0.02 K, showing that even macroscopically occupied thermal states (at T = 0.01 K the mean thermal occupation number of a mode at frequency ω_m is larger than 10³) can accommodate the set entanglement (cf figures 4(a) (inset) and (b)).

Zoom In Zoom Out Reset image size

We now assume that a three-mode optical driving field is prepared to drive the dynamics of a three-site optomechanical network. While the principles of the dynamics at the single optomechanical site are the same as those for the two-site case, the description of the optical resource needs some care. The covariance matrix of a pure symmetric three-mode state is parameterized as [24]

$$\begin{equation} {\bm \sigma }=\frac{1}{2}\left( \begin{array}{@{}cccccc@{}}a_1&0&c^+_{12}&0&c^+_{13}&0\\ 0&a_1&0&c^-_{12}&0&c^-_{13}\\ c^+_{12}&0&a_2&0&c^+_{23}&0\\ 0&c^-_{12}&0&a_2&0&c^-_{23}\\ c^+_{13}&0&c^+_{23}&0&a_3&0\\ 0&c^-_{13}&0&c^-_{23}&0&a_3 \end{array}\right) \end{equation} \tag{ 16 }$$

with

$$\begin{eqnarray} &&c^\pm _{ij}=\frac{\sqrt{[({a_i}-{a_j})^2-({a_k}-1)^2] [({a_i}-{a_j})^2-({a_k}+1)^2]}}{4 \sqrt{{a_i}{a_j}}}\nonumber\\ \end{eqnarray} \tag{ 17 }$$

$$\begin{eqnarray} &&\pm\,\frac{\sqrt{[({a_i}+{a_j})^2-({a_k}-1)^2] [({a_i}+{a_j})^2-({a_k}+1)^2]}}{4 \sqrt{{a_i}{a_j}}}.\nonumber\\ \end{eqnarray} \tag{ 18 }$$

In order to ensure physicality of the corresponding state, the parameters a_j should satisfy the triangular inequality

$$\begin{eqnarray} &&\fl |{a_i}-{a_j}|+1\le {a_k}\le {a_i}+{a_j}-1 (i\ne j\ne k=1,2,3 \nonumber\\ &&{\rm with}\quad a_j\ge 1). \end{eqnarray} \tag{ 19 }$$

The analysis presented in section 2 has revealed that pure and symmetric optical states enable better performances of our scheme, at least for the two-site case. We now conjecture that such condition extends to the three-site situation and restrict our attention to the case of a_j = a (∀j); therefore, considering the covariance matrix of the pure three-mode optical state

$$\begin{equation} {\bm \sigma }_{\rm opt}=\frac{1}{2} \left( \begin{array}{@{}cccccc@{}}a & 0 &{x^+} & 0 &{x^+} & 0 \\ 0 & a & 0 &{x^-} & 0 &{x^-} \\ {x^+} & 0 & a & 0 &{x^+} & 0 \\ 0 &{x^-} & 0 & a & 0 &{x^-} \\ {x^+} & 0 &{x^+} & 0 & a & 0 \\ 0 &{x^-} & 0 &{x^-} & 0 & a \end{array} \right) \end{equation} \tag{ 20 }$$

with $x^\pm =( a^2-1 \pm \sqrt{1 - 10 a^2 + 9 a^4})/(8 a)$. Such a conjecture has been verified a posteriori by taking covariance matrices with lower degrees of symmetry (such as cases where the variances of only two out of three modes are identical) and checking that the corresponding performances of the protocol were inferior to those of the fully symmetric case. The state with covariance matrix given by equation (20) is genuinely tripartite entangled according to the measure given by the so-called minimum residual Gaussian contangle [27] ${\mathcal E}_{\rm res}$. The Gaussian contangle ${\mathcal E}_{a|b}$ between the pair of subsystems a and b is defined as the convex roof (performed over all the pure Gaussian states with covariance matrix ${\bm \sigma }_{\rm Gauss}$ strictly smaller than the covariance matrix of the state at hand) of the squared logarithmic negativity [24] and is a proper entanglement monotone. For a set of modes μ = {i, j, k, l, ..., z}, the Gaussian contangle satisfies the monogamy constraint ${\mathcal E}_{i|jkl..z}-\sum _{\mu =\lbrace j,k,l,\ldots ,z\rbrace }{\mathcal E}_{i|\mu }\ge 0$, together with the analogous relations obtained by permuting the mode indices at will. With this, for three modes (i, j, k), ${\mathcal E}_{\rm res}$ is defined as

$$\begin{equation} {\mathcal E}_{\rm res}=\min _{\hat{\Pi }}\hat{\Pi }[{\mathcal E}_{i|jk}-{\mathcal E}_{i|j}-{\mathcal E}_{i|k}] , \end{equation} \tag{ 21 }$$

with $\hat{\Pi }$ the operator that permutes the indices {i, j, k}. The residual entanglement of the symmetric state with covariance matrix $2{\bm \sigma }_{\rm opt}$ (the factor 2 is used to make the notation uniform with the formalism of [27]) can be computed analytically as ${\mathcal E}_{\rm res}={\rm arc\,sinh}^2[\sqrt{a^2-1}]-(1/2)\ln ^2[2a^2+4ax^-]$. Moreover, the state has non-null entanglement between any two modes, growing (but not indefinitely) with a. Driving fields of this type are within the grasp of current experimental possibilities [24].

The corresponding mechanical covariance matrix is the extension of (13) to three mechanical modes. Throughout this section, we maintain the very same values of the optomechanical parameters given in table 1. Needless to say, the output covariance matrix ${\bm \sigma }_{\rm mech}$ is physical as far as ${\bm \sigma }_{\rm mech}+\mathrm{i}\Omega _3/2\ge 0$ (which is equivalent to the condition on the symplectic spectrum $\min {\rm eig}|\mathrm{i}\Omega _3{\bm \sigma }_{\rm mech}|>1/2$). As for the entanglement, given the tripartite nature of the mechanical network at hand, the situation is more involved than in the two-site case. Indeed, we can distinguish the case of bipartite one-versus-one-mode entanglement, achieved by tracing out any one of the three mechanical modes, from the one-versus-two-mode entanglement, which we will use to infer the tripartite inseparable nature of the mechanical state.

Let us start with the first instance, which implies the trace over the degrees of freedom of one of the modes of our mechanical network. Needless to say, the enforced symmetry of the optical input state guarantees full symmetry of the mechanical one, thus making all the two-mode reductions identical to each other. Formally, this is equivalent to considering only the blocks in the mechanical covariance matrix ${\bm \sigma }_{\rm mech}$ containing the degrees of freedom of the modes that are not traced out. The entanglement is then calculated by using the very same formal apparatus described in section 2, and we can characterize the output mechanical entanglement of the two-mode-reduced states versus the analogous quantity for the input optical one (cf figure 5(a)). In analogy with the case of two-site networks, we observe both non-monotonic behaviour of the output mechanical entanglement and the existence of a threshold in the input optical entanglement for the process to be successful. However, at variance with what has been observed in section 2, the curve of the output entanglement is asymmetric with respect to the input one. Needless to say, as the reduced two-site system is mixed, we expect the distribution performances in this case to be contained within the optimal boundary embodied by a TMSV. This is shown in figure 5(b) and its inset.

Symmetry is fully reinstated when one looks at the one- versus-two-mode splitting of the network. The entanglement in such a configuration (which is again invariant under permutations of the mode indices due to the overall symmetry of the three-mode mechanical covariance matrix) is detected by the violation of the Heisenberg–Robertson condition by the partially transposed covariance matrix $P_{i|jk}{\bm \sigma }_{\rm mech}P_{i|jk}$, where P_i|jk is the matrix that inverts the sign of momentum of mode i, thus implementing the partial transposition with respect to such an element of the mechanical network. The entanglement in such a one-versus-two-mode split is then quantified as ${\mathcal L}_{\mathcal N}=\max [0,-\ln 2\tilde{\nu }_-]$ with $\tilde{\nu }_-=\min {\rm eig}|\mathrm{i}\Omega _3P_{i|jk}{\bm \sigma }P_{i|jk}|$ (as the partially transposed covariance matrix has only one symplectic eigenvalue that can be smaller than 1/2) and is plotted, against the input optical entanglement, in figure 5(c). Not surprisingly, the similarity with figure 2(b) is only apparent and the one-versus-two-mode entanglement is strictly contained within the one-versus-one-mode entanglement.

Zoom In Zoom Out Reset image size

With this at hand, we can proceed to the qualitative characterization of the tripartite entanglement within the system at hand. Following the classification by Giedke et al [18], a tripartite Gaussian state whose one-versus-two-mode bipartitions are inseparable is tripartite inseparable. Looking at the results shown in figure 5(c), this occurs, in our system and for the set of optomechanical parameters at hand, for quite a large set of input optical states, thus showing that the multipartite entanglement used as a resource in our protocol is indeed able to generate tripartite entangled states of the mechanical network, in a sense signalling the success of the protocol itself.

Yet, a point that remains to be addressed is the quantification of the genuine tripartite entanglement within the full mechanical state. As the state of the mechanical network is overall mixed, the calculation should be performed with care. Without entering unnecessary details related to the rather technical approach to the calculation of the residual genuine tripartite entanglement in the state of a system of bosonic modes, we mention that the results of such calculations show that the genuine tripartite entanglement is indeed transferred to the mechanical network, albeit in an inefficient manner (the average tripartite mechanical entanglement being two orders of magnitude smaller than the corresponding input one, in general). Yet, this is sufficient for the scope of this investigation, which aimed at providing a formal analysis of the distribution performance in the simplest instance of multipartite optomechanical network.

In this paper, we have addressed extensively the problem of entanglement distribution in an optomechanical network that is driven by a multipartite entangled optical resource. By using the most general parameterization of (generally mixed) Gaussian optical states, we have found the conditions under which the two-site mechanical entanglement is maximally distributed, identifying the role that symmetries and local/global state purity of the optical resource play in this process. We have demonstrated the efficiency of the protocol, which can be implemented using state-of-the-art parameters and the possibility of distributing the genuine tripartite mechanical entanglement.

Many questions remain to be addressed, for instance related to the ability of the protocol to pass general quantum correlations that go beyond the quantum entanglement. This is an interesting opportunity in light of the preliminary analysis performed in [7], where such quantum correlations have been shown to be more robust than entanglement to the increase of temperature in the mechanical setup. Moreover, it would be interesting to find whether optical postselection is able to improve the distribution mechanism, giving rise to more substantial genuinely multipartite mechanical entanglement. Finally, we intend to perform a classification of the distributed mechanical entanglement at assigned values of the global and marginal purity, along the lines of some of the studies in [14], so as to check if a hierarchy of mechanical states arises from the distribution process.

MP is supported by the UK EPSRC through a Career Acceleration Fellowship and the 'New Directions for EPSRC Research Leaders' initiative (EP/G004759/1). LM acknowledges the EU under a Marie Curie IEF Fellowship (grant agreement 300368). JL thanks the School of Mathematics and Physics, Queen's University Belfast, for a Purser Studentship.