Preparation and Verification of Two-Mode Mechanical Entanglement Through Pulsed Optomechanical Measurements

We propose a protocol how to generate and verify bipartite Gaussian entanglement between two mechanical modes coupled to a single optical cavity, by means of short optical pulses and measurement. Our protocol requires neither the resolved sideband regime, nor low thermal phonon occupancy, and allows the generation and verification of quantum entanglement in less than a mechanical period of motion. Entanglement is generated via effective two-mode mechanical squeezing through conditioning position measurements. We study the robustness of entanglement to experimental deviations in mechanical frequencies and optomechanical coupling rates.

So far, optomechanically-induced mechanical entanglement has been achieved in systems that require the mechanical frequency ω to be much larger than a pulse or cavity bandwidth κ [30,31,32]. This resolved sideband regime can be challenging to achieve as it requires a quality factor Q em of the electromagnetic mode to be higher than the ratio of electromagnetic frequency ω em over mechanical frequency ω. For systems that do not resolve motional sidebands (κ ω), an alternative can be sought through conditioning the mechanical state with snapshot measurements of the mechanical position, using optical pulses shorter than the mechanical period [33]. Provided that the mechanicalposition-induced phase shift can be well resolved above the optical shot noise, significant mechanical cooling or squeezing by measurement can be achieved [34]. Indeed, using this technique, mechanical squeezing approaching the quantum regime has been achieved, thus providing a promising route for mechanical entanglement [35,21,22].
In [21], a protocol is proposed for mechanical entanglement and verification of two spatially-separated optomechanical resonators, which can be essential for entanglement distribution applications. Routes for generation and verification of entanglement by combining pulsed measurements on the separated systems are then identified. Reference [22] proposes to create two-mode entanglement of a single device using measurements in the so-called stroboscopic limit, with a protocol that does not explicitly depend on the mechanical mode frequencies. The sensitivity of the expected degree of entanglement to experimental deviations such as optical matching of spatially-separated devices, or optomechanical coupling mismatch between the mechanical modes, remains unexplored.
In this work, we explore the ability to generate and verify entanglement of two mechanical modes of a single optomechanical system, in a minimal two-pulse measurement sequence, lasting for less than a mechanical period of both resonators. As such, our methods differs from those using measurements on separate cavities [21] or using stroboscopic measurements [22] of a single device. These two mechanical modes can be two mechanical normal modes of a single body, or separately engineered, for instance with two membranes in a Fabry-Pérot cavity [36,37,38], or a sliced photonic crystal nanobeam with two flexural modes of different frequencies [39] as sketched in Fig. 1a. We show in section 2 that the requirements to generate entanglement is as demanding as a pulsed-measurement induced single mode ground-state cooling would require, in terms of optomechanical features, which is within reach with state of the art optomechanical systems [35]. Section 3 investigates the robustness of the protocol to experimental deviations. We find that our protocol necessitates a relative frequency tuning of the resonators at percent-level precision, and optomechanical coupling equality at the ten percent level, which make this entanglement scheme realistic from a practical perspective in various opto-and electromechanical platforms.

Qualitative explanation
The essence of a conditional pulsed optomechanical measurement at a time t, is that one can obtain a snapshot of information about the mechanical position of a resonator x(t) by measuring the outgoing light pulse with which it interacted. Reading out the phase quadrature P L of the light pulse allows one to measure mechanical position with shot noise imprecision x zpf /χ, where x zpf is the zero-point fluctuation standard deviation and χ ∼ η N p g/κ is the measurement strength [33], scaling with the optical linewidth κ, the vacuum optomechanical coupling g and the number N p of photons in the pulse interacting with the cavity, i.e. the number of incident photons minus input coupling losses. The efficiency η incorporates detection and outcoupling efficiency. For large enough χ, one can achieve backaction-evading measurement of the mechanical position, as light shot noise backaction is imparted into the unmeasured momentum quadrature.
Let us consider now that we have two mechanical modes b j=1,2 with frequencies ω j , quadratures 2i, and positions x j (t) = X j cos(ω j t) + P j sin(ω j t). Those two modes are assumed for now to interact with the same strength with the optical mode. Then, any conditional pulsed measurement will give information about the total mechanical position x(t) = x 1 (t) + x 2 (t), leading to a joint measurement of the two resonators at once. This implies strong measurementinduced correlations between the two mechanical modes, and eventually entanglement if the resolution of the collective position is lower than the amplitude of the collective zero-point motion (ZPM). Namely, entanglement of the two-mechanical modes will be achieved if variances of the collective mechanical quadratures are such as where the collective mode quadrature basis (X + , P + , X − , P − ) links to the canonical mode quadrature basis (X 1 , P 1 , X 2 , P 2 ) by a beamsplitter-like mixing X ± = (X 1 ± X 2 ) / √ 2 and P ± = (P 1 ± P 2 ) / √ 2. The measured mechanical position x(t) = X 1 cos(ω 1 t) + P 1 sin(ω 1 t) + X 2 cos(ω 2 t) + P 2 sin(ω 2 t) scans over time multiple collective quadratures because of the frequency mismatch of the two resonators ω 1 = ω 2 . Then, a measurement of the mechanical position at an arbitrary  Figure 1. (a) Optical pulses shorter than the mechanical period are sent to an optomechanical system. Copies of these pulses are derived and play the role of a local oscillator in the homodyne detection of the outgoing pulses. Adjusting the relative phase Φ between the two optical paths selects the optical phase quadrature P L detection. The two mechanical modes can be the two flexural modes of a sliced nanobeam photonic crystal cavity, two octaves detuned. (b) Circuit diagram of the mechanical entanglement and verification protocol. Two pulsed optomechanical measurements of the mechanical position delayed by π/2ω 1 allows to get information and then reduces the thermal uncertainty of the X + and P − collective quadrature respectively, leading to entanglement of the two resonators in less than a mechanical period. The verification protocol involves an extra pulsed measurement at various delays to reconstruct the covariance matrix of the two-mechanical-modes system. (c) Wigner functions of the partially traced two-mode system over the canonical mode basis (b 1 , b 2 ) and the collective eigenmode basis (b, c). Each plot corresponds to different steps of the entanglement protocol. We assume arbitrary outcome values of X + = 15 and P − = 0 and n th = 500. Measurements affect the collective mechanical quadrature statistics, without significantly affecting the mechanical modes taken individually.
time t reduces the uncertainty of an arbitrary collective quadrature X(t) = x(t)/ √ 2 given by Eq. (2). We label b(t) the time-dependent collective mode whose amplitude quadrature is X(t), and c(t) to complete the two-mode basis.
Starting from a measurement of the position at the time origin x(0) = √ 2X + , it is not always possible to measure the quadrature P − , which is necessary to fulfill Eq. (1). Indeed, the closer the mechanical frequencies, the longer we need to wait between measurements to measure significantly different collective quadratures, leading to a potential for decoherence to negatively influence the conditional state purity. On the other hand, too much different frequencies can hardly be designed to equally couple to the same optical mode. Then, achieving entanglement shall necessitate multiple measurements to compensate for these drawbacks [22]. In this work, we focus on a minimal two-pulse sequence that allows to generate entanglement, independently of the initial thermal occupation of the resonators. For instance, the situation were ω 2 /ω 1 = 3 allows to entangle and verify in less than a mechanical period of both resonators, while ensuring frequency in the same order of magnitude. This choice of mechanical frequency ratio leads to the two-pulse entanglement protocol schemed in Fig. 1b. Starting from uncorrelated mechanical resonators (step 0 ), the protocol consists of a measurement at the time origin x(0) = √ 2X + (step 1 ), followed by waiting a delay of a quarter mechanical period of the lowest frequency resonator τ = π/2ω 1 (step 2 ) that leads to a mechanical position of x(τ ) = √ 2P − . A second conditioning measurement at this point (step 3 ) generates two-mode squeezing of the mechanical resonators, for large enough measurement strength. Figure 1c shows Wigner functions of the partially-traced mechanical states in the canonical (b 1 , b 2 ) and collective (b(t), c(t)) mode basis, conditioned on the measurements, at the different steps of the entangling sequence. Starting from an initial wide thermal distribution ( 0 ), uncertainty of the amplitude quadrature of mode b(0) is strongly reduced by the first conditioning measurement ( 1 ). The two canonical mechanical modes b 1,2 now have some correlation, which we strengthen with the second entangling pulse. Right before the second entangling pulse ( 2 ), the mechanical resonators are out of phase, at opposite velocity. At this point, the mechanical collective mode b(π/2ω 1 ) has the statistics of the collective mode that was unaffected by the first conditional measurement, and its amplitude quadrature is X(π/2ω 1 ) = P − . After this second conditional measurement ( 3 ), the mechanical modes are two-mode squeezed. These successive measurements always purify and classically entangle the mechanical resonators, and quantum two-mode squeezing can be achieved if the measurement strength is large enough.

Quantitative description
We introduce X can = (X 1 , P 1 , X 2 , P 2 , X L , P L ) T and X col = (X + , P + , X − , P − , X L , P L ) T the quadrature vectors of this bosonic three-mode system in the canonical and collective mode basis, respectively. We note that all the quadratures are defined in the frame rotating at the frequency of their mode, such that their free dynamics is timeindependent. The linearized Hamiltonian in the frame rotating at the cavity frequency reads [40] where we assume, for now, that the linearized optomechanical coupling rates of the two mechanical modes are equal to g. In the case of a short optical pulse containing N p photons over a duration τ p 1/ω 1,2 at t = 0, the evolution operator right after the pulsed interaction writes where the measurement strength χ quantifies the amount of information about mechanical position one can extract from an optical phase quadrature P L measurement Note that we omitted in the evolution operator the deterministic kick on conjugated mechanical momentum P 1 + P 2 induced by radiation pressure of the mean optical field [33], as it only features a classical displacement. We note that the other orthogonal collective mode (X − , P − ) is not affected by this measurement.
Measuring the optical phase quadrature at t = 0 reduces uncertainty of the X + quadrature defined previously, while adding backaction noise on P + . More generally, a pulsed measurement at any arbitrary time t always results in the measurement of a combination of quadratures of the two mechanical resonators. We define the timedependent eigenmode quadrature basis of the Hamiltonian R are unitary rotation matrices in the mechanical basis, at the average mechanical frequency ω = (ω 1 + ω 2 ) /2 and relative mechanical frequency Ω = (ω 1 − ω 2 ) /2. We denote by b(t) (c(t)) the mechanical collective mode with quadratures X(t) and P (t) (Y (t) and Q(t)), explicitly affected (unaffected) by the interaction. Hamiltonian (3) writes H(t) =hgX L (t)X(t) such that any pulsed measurement of the optical phase quadrature at a time t gives information about X(t) while adding backaction on P (t), and leave the other collective mechanical mode unchanged: The pulsed interaction generates strong correlations between the outgoing light and the mechanical modes. The covariance matrix σ = 1 2 {(X − X ), (X − X ) T } of the three-mode system reads where S U (t) is the symplectic transformation associated to the operator U (t) [41,42]. The first additive term is backaction noise and the second is measurement-induced optomechanical cross-correlations : When measuring the optical phase of the output light, we get information about the mechanical motion. This measurement M changes our knowledge about the mechanical modes, via the transformation [41] where σ HD describes the homodyne measurement and σ is decomposed over the mechanical (σ m ) and optical (σ L ) basis, and their cross-correlations (σ mL ). If we consider to have an ideal non-lossy homodyne detection locked on the phase quadrature of the outgoing light field, we have [41] (σ L + σ HD ) −1 → diag 0, 1/σ 22 Finally, the effect of a single measurement at time t can be summed up in the mechanical covariance matrix in the canonical basis as σ m = P (t)M (P (t) −1 σ m P (t)) P (t) −1 , where P (t) = R ω (t) R Ω (t). Then, the entangling sequence of Fig. 1b leads to the exact following covariance matrix in the canonical mode basis: where n i = 1 2 + n th,i is the total occupation of the mechanical mode i, and σ z the z−Pauli matrix. The asymmetry of this covariance matrix between the two mechanical modes (A 1 = A 2 ) originates from the initial unequal thermal occupation of the two modes (n 1 = n 2 ). The covariance matrix in the canonical basis can be explicitly written as a two-mode squeezed state by parameterizing its eigen-decomposition as diag (n e e −re , n e e re , n e e re , n e e −re ) with n e e ±re = α + + β ± α 2 − + β 2 /2, α ± = n 1 ± n 2 1 + 2 (n 1 + n 2 ) χ 2 and β = χ 2 1 + 4n 1 n 2 1 + 2 (n 1 + n 2 ) χ 2 .
In the relevant regime of initial high thermal occupations n i {1, χ} these eigenvalues simplify to n e e −re = 1 4χ 2 and n e e +re = 2n 1 n 2 n 1 + n 2 .
The first term corresponds to the measurement resolution associated with the X + and P − collective quadratures measurement. The second term corresponds to the standard deviation of the unmeasured collective quadratures, with thermal occupation averaged over the two mechanical modes.
The duration of the pulse is considered to be negligible compared to the decoherence decay rate. However, decoherence of the mechanical modes is taking place in between measurement pulses, i.e. during the free evolution of the mechanical modes. We take it into account by considering the mixing of the prepared state with a thermal state at rate Γ 1,2 : with G = diag e −Γ 1 t , e −Γ 1 t , e −Γ 2 t , e −Γ 2 t and Γ i n i the heating rate of mechanical mode i [43]. For sake of clarity, we focus on the relevant regime of high-quality thermal resonators, i.e. with mechanical quality factors Q i = ω i /Γ i 1 and n 1 = 3n 2 = n 1. In this regime, the full entangling sequence of Fig. 1b leads to the following covariance matrix eigen-decomposition : where we have developed at the leading order in 1/Q i . This covariance matrix can be parametrized with the relevant two-mode squeezing parameters diag (n + e −r + , n + e r + , n − e r − , n − e −r − ). The collective quadratures X + and P − (n ± e −r ± terms) are squeezed below the ZPM amplitude, provided that the measurement strength exceeds unity and that decoherence does not induce noise on the X + quadrature larger than the ZPM amplitude in between the two pulses. The two other collective quadratures (n ± e +r ± terms) are only affected by measurement quantum backaction, which is usually very weak with respect to the wide thermal variance. We note that P − is in principle not affected at all by decoherence because the covariance matrix above is written right after the second measurement.
We describe the amount of entanglement by means of the logarithmic negativity E N . This quantity estimates the maximal number of singlet Bell states one could get with an entanglement distillation protocol. Besides the advantage of being measurable [44], this quantity is widely considered in both discrete [45,46] and continuous variable [47] entanglement characterization. It characterises the deviation from positiveness of the two-mode system after a single-mode transposition has been applied. After transforming the covariance matrix (16) back to the canonical basis, an analytical derivation of the logarithmic negativity of a bipartite Gaussian system (16) leads to [44] E ent N = max 0, −log 2 2 n − n + e −r − −r + .
Then, entanglement can be efficiently generated as soon as the geometric mean of the two squeezed quadrature variances is lower than the ZPM amplitude. From now on, for sake of simplicity, we suppose the mechanical quality factors of both resonators to be equal Q 1 = Q 2 = Q, meaning that the thermal decoherence rates Γ i n i are equal too. Figure 2a shows the logarithmic negativity as a function of the mechanical quality factor, for different measurement strength. In the relevant approximated regime of Eq. (16), the criteria to achieve generation of entanglement are  Figure 2. (a) Estimated collective variance from the verification protocol against mechanical quality factor Q. Color contrast scales with the measurement strength χ, valued from 1 (shaded color) to 5 (dusky color). Noise variance of the X + and P − quadratures can be strongly reduced even below the ZPM amplitude, while the noise variance of the X − and P + quadratures are slightly affected by backaction, because of the relatively large thermal phononic occupation n th = 10 4 χ 2 . Thermal decoherence between pulses tends to thermalize the squeezed thermal state. (b) Logarithmic negativity of the state generated after entanglement sequence (continuous line) and the evaluated state by verification sequence (dashed line). The verification protocol always underestimates the amount of entanglement present in the two-mode mechanical state, and tends to optimally estimate it in the low-decoherence limit.
These requirements are not more demanding than those already required for conditional ground state cooling or squeezing of a single mechanical mode using conditional pulsed measurements [33]. Recent experiments with a single flexural mode of a sliced nanobeam photonic crystal structure with moderate mechanical properties (ω m /2π 3 MHz, κ/2π 20 GHz, g/2π 25 MHz, η 1 %, N p = 60, Q 3 × 10 5 ) have already reached a measurement strength of χ = 0.11 [35], highlighting that achieving the criteria above is within reach of near-term experiments. We note that a Brownian model may be used to model the mechanical decoherence, in which the mechanical momentum is displaced through interactions with the bath, instead of the symmetric beamsplitter model described by Eq. (15). If one does this, the criteria to generate entanglement will be less stringent than Eq. (18) as the optical phase quadrature measurement contains information about mechanical position only, which is unaffected by a Brownian model until free evolution of the mechanical mode evolves the momentum noise into position.
We finally note that two-mode squeezing can in principle be achieved with a single pulse measurement, in case of low thermal occupation. Indeed, if one measures the X + quadrature only, the conditioned covariance matrix eigen-decomposition writes diag (1/4χ 2 , n/2, n/3, n). One can show that this two-mode squeezed state is quantum entangled if the challenging criteria E ent N > 0 ⇔ n < 4χ 2 is fulfilled.

Verification of entanglement
A complete verification protocol of the entanglement could require one to reconstruct the full two-mode covariance matrix of the mechanical state. However, we can only access a restricted ensemble of collective quadratures X (t) by means of mechanical position measurement. Nevertheless, we use the fact that the covariance matrix (16) is diagonal in the X col basis to selectively detect each of its diagonal components [48]. The X + (P − ) quadrature variance σ ent X + (σ ent P − ) can be measured independently and selectively by waiting a delay t = π/2ω 1 (t = π/ω 1 ) from the last entangling pulse. Namely, multiple measurements would distribute over a Gaussian with width σ ver X + = D σ ent X + , π/2ω 1 (σ ver P − = D σ ent P − , π/ω 1 ). The two other collective quadratures cannot be measured selectively. One option is to measure the variance of the superposition σ X − +P + by waiting a delay t = π/4ω 1 from the last entangling pulse, and assume that the thermal distribution of the superposition is the same as the distribution of X − and P + taken individually : σ ver X − = σ ver P + = D σ ent X − +P + , π/4ω 1 . This assumption is exact in the ideal case of resonators detuned by two octaves. We will investigate in the next section to which extent this approximation is valid for entanglement verification for detuned resonators.
The reconstructed mechanical covariance matrix in the canonical mode basis from verification is then written as σ ver m = diag σ ver X + , σ ver X − , σ ver P + , σ ver P − . In the relevant parameter regime of Eq. (16), it writes Figure 2a shows the collective quadrature variances evaluated by verification protocol. Collective quadratures that were conditionally measured during the entangling sequence are drastically reduced below the thermal level. We note that both variances ∆X + and ∆P − are equally evaluated because of the same time delay -i.e. decoherence thermal mixing -between their respective conditional measurement pulse and verification pulse. The logarithmic negativity of this measurable covariance matrix will always underestimate the actual amount of entanglement right after the entangling sequence E ver The criteria to measure entanglement are then 2πn where the first condition can be interpreted as the requirement for the conditional quadrature variance to remain below the ZPM amplitude after half a mechanical period waiting delay, starting from an infinitely squeezed state. The logarithmic negativity of the evaluated mechanical state together with the entangled state is plotted in Fig. 2b. The requirement on the mechanical quality factor is more stringent to detect entanglement than to generate it because of mechanical decoherence and rethermalization between the entangling and verification pulses.

Robustness of the protocol to experimental imperfections
The protocol we propose requires having two mechanical resonators tuned at two octaves from each other, that moreover couple to an optical cavity with the same rate. These requirements can be experimentally challenging if one imagines to entangle two higher order modes of the same mechanical resonator for example. We discuss in this section how stringent these requirements are, by considering deviations of the mechanical frequencies (section 3.1) and the coupling strength (section 3.2) from the ideal protocol. Figure 3a shows numerically computed logarithmic negativity of the mechanical state generated, for a wide range of two-pulse entanglement sequence, and over a wide range of mechanical frequencies ratio. Plots are obtained in the case of a modest measurement strength χ = 2. Each blue area is centered on an efficient two-pulse entanglement sequence characterized by the doublet (k, l) ∈ Z + satisfying :
These two conditions mean that the second pulse should be applied after a delay τ after which the two mechanical oscillators are in opposite phase and at maximum velocity.
All these protocols are in principle equivalent. However, increasing values of k require longer waiting times, thus causing more impact of thermal decoherence. Moreover, we can see qualitatively that the width of the blue regions in the ω 2 /ω 1 direction decreases, with increasing values of ω 2 /ω 1 . We have chosen to focus in the former section on the (k = 0, l = 0) protocol, because of its lowered sensitivity to decoherence and the practical advantage of close mechanical frequencies. We moreover observe numerically that the entanglement generated by this protocol is the most resilient to experimental detuning of the resonators. Last but not least, the width of the blue regions indicates that one can compensate for any mechanical frequency detuning by shifting the delay between the two pulses of the entanglement sequence. Figure 3b compares the logarithmic negativity reachable by the entanglement sequence with the one obtained with the verification sequence, in the case of the (k = 0, l = 0) protocol. The entanglement sequence shows strong robustness to mechanical detunings, as it can be compensated by adjusting the delay between the entangling pulses, by following the black line of Fig. 3a. Similarly, entanglement verification can be adapted by adjusting the delays of the verification pulses, such that entanglement can be optimally verified. Figure 3b shows that the requirement on the mechanical frequency tuning to verify entanglement is more stringent than for its generation, on the order of ±5% for a measurement strength of χ = 2. Moreover, the verification protocol imprecision requirement is not centered on the ideal case where ω 2 = 3ω 1 . This asymmetry is because ω 2 ≥ 3ω 1 protocols are slightly more affected by decoherence because the optimally adjusted pulse delays are longer than π/2ω 1 . We interpret the fact that verification sequence is much more sensitive to frequency mismatch than the entanglement sequence because of the restricted number of component of the covariance matrix we can actually measure, while assuming all other non-diagonal terms vanish.

Optomechanical coupling imprecision
In this section we investigate the effects of a mismatch between the measurement strengths for the two mechanical modes, which we model using the more general Hamiltonian This can be taken into account in by applying an extra change of basis to write the mechanical eigenmodes of the Hamiltonian. Let's write g = (g 1 + g 2 )/2 and g = (g 1 −g 2 )/2. The Hamiltonian then writes H(t) =hX L (gX(t) + g Y (t)) / √ 2. Then, any measurement that was planned to focus on X + for instance will always be associated with measurement of X − as well, with backaction on both P + and P − quadratures. This tends to reduce the squeezing of the two modes. The eigenmodes of the Hamiltonian are defined using an extra transform breaking the symmetry between the two modes (23) Figure 4a shows the influence of this deviation on the logarithmic negativity of the mechanical state generated and verified by the (k = 0, l = 0) protocol, for ideally twooctaves-detuned resonators. The protocol generates maximal entanglement when the measurement strengths are equal. In the relevant regime of low decoherence and large thermal occupation, the amount of entanglement scales with the measurement strength mismatch χ = χ 2 /χ 1 − 1 as At first order, entanglement tends to be stronger because of the larger average measurement strength of the two resonators, displacing the maximal reachable entanglement to higher values of χ 2 /χ 1 . Nevertheless, the maximum of entanglement is reached within a χ range of ± log 2 (2χ 2 1 ) /2χ 4 1 , about ±10 − 30 % for moderate measurement strengths considered in Fig. 4a.
Noting that the two-pulse entanglement protocol (0,0) is robust to experimental deviations, Fig. 4b extends the idea of two-pulse protocol to mechanical resonators with any measurement strength and frequency detuning. For each two-mode configuration, we computed the maximally reachable entanglement by means of a two-pulse sequence, if resonator ω 1 has a measurement strength χ 1 = 2. We observe that the maximal entanglement is set by the lowest measurement strength of the two resonators, and that many configurations actually leads to efficient mechanical entanglement. By adjusting the verification sequence to address different frequency mismatches, these results show that such a two-pulse protocol can be implemented in many experimental configurations to generate mechanical entanglement.

Conclusion
We propose and analyze in this work the use of two-pulse conditional measurement for bipartite Gaussian entanglement of two mechanical resonators coupled to an optical cavity. We focused on a specific protocol where the two resonators are detuned by two octaves, and couple to the optical cavity with equal strength. This specific protocol allows generation and verification of entanglement in less than a mechanical period, and shows robustness to experimental deviations of the mechanical frequency mismatch or optomechanical coupling mismatch. We showed that the two-pulse protocol can be extended to wider range of systems, by carefully choosing the two-pulse entanglement sequence and the verification sequence. The criteria to generate entanglement in terms of measurement strength and mechanical quality factor is reachable with current state of the art optomechanical designs. Moreover, this type of protocol may be extended to higher order entanglement generation in a local network of resonators, exploiting for instance higher order mechanical modes.