Quantum manipulation of a two-level mechanical system

We consider a nonlinearly coupled electromechanical system, and develop a quantitative theory for two-phonon cooling. In the presence of two-phonon cooling, the mechanical Hilbert space is eﬀectively reduced to its ground and ﬁrst excited states, allowing for quantum operations at the level of individual phonons and preparing nonclassical mechanical states with negative Wigner functions. We propose a scheme for performing arbitrary Bloch sphere rotations, and derive the ﬁdelity in the speciﬁc case of a π -pulse. We characterise detrimental processes that reduce the coherence in the system, and demonstrate that our scheme can be implemented in state-of-the-art electrome-chanical devices.

It was recently shown [48] that these problems can be alleviated for carefully designed systems. Here, we consider a similar setup and demonstrate that one can realistically achieve two-phonon cooling [49], which is a process that decreases motional quanta by annihilating two phonons at the time. Contrary to the QND detection of the phonon number studied in Ref. [48], the two-phonon cooling is frequency selective and only talks to a single mechanical mode, making it more experimentally accessible. As a paradigmatic example of non-Gaussian operations, we describe how to create and manipulate a mechanical qubit, and use the two-phonon cooling to deterministically prepare states with negative Wigner functions. The performance of our protocol is characterized by a single parameter λ, which is uniquely determined by the experimental setup and, importantly, includes detrimental multi-mode contributions . We study the RLC circuit in Fig. 1(a), where one of the capacitor plates is a moving mem-brane. By changing the circuit's capacitance, its motion induces an electromechanical interaction. The considered asymmetric oscillatory mode ensures the suppression of the linear optomechanical coupling, making the higher-order (i.e. non-Gaussian) interactions dominant. The small magnitude of the quadratic coupling, however, obliges us to consider all possible residual linear effects that are detrimental to the twophonon cooling. Following Ref. [48], we identify two malicious contributions. The first originates from fabrication imperfections, and couples the "symmetric" electrical mode (indicated by subscript "s" in Fig. 1) with the membrane. The second arises from redistribution of charges on the capacitor plates, which we model by introducing a loop in Fig. 1(b) with parasitic resistances R and inductances L. The resulting "asymmetric" electrical mode (subscript "a") interacts linearly with the membrane. This work is organized as follows. In Sec. 2, we present the Master equation describing the system's dynamics and determine the conditions (summarized in the parameter λ) under which the two-phonon cooling is feasible. In Secs. 3 and 4, respectively, we discuss experimental implementations and introduce a protocol for coherently manipulating the membrane's quantum state within the reduced Hilbert space provided by the two-phonon cooling. Specifically, we consider a π rotation of an arbitrary initial mechanical state and characterize its fidelity in terms of λ. We summarize our results and present an outlook for future works in Sec. 5.

Two-phonon cooling
Here, we derive the two-phonon cooling rate Γ 2 and compare it with the thermalization rates from detrimental linear effects. The optomechanical interaction described by the couplings g 1 and g 2 is derived by expanding the inverse capacitances of the two halves in Fig. 1 The sign uniquely identifies the capacitor in the circuit, as indicated in Figs. 1(a) and (b). The operatorx is the mechanical displacement away from equilibrium, and ω s ≡ are the frequencies of the two electrical modes, with ω a ω s . Furthermore, δ takes into account the effects of nano-fabrication inaccuracies. Without these (δ = 0), the sign difference in the expansion suppresses the linear interaction between the membrane and the symmetric mode. The Hamiltonian of the system in Fig. 1 Here,Q µ andΦ µ correspond, respectively, to the charge sum (µ = s) or difference (µ = a) of the two capacitors, and the magnetic fluxes flowing in the circuit. The mechanical annihilation (creation) operator isb (b † ), and the membrane frequency is ω m . The remaining parameters are shown in Fig. 1(b) and described in its caption. The interaction Hamiltonian in Eq. (1) is derived for the symmetric mode being driven with a strong field of frequency ω s − 2ω m . This allows disregarding all membrane modes except for the chosen asymmetric one. With this driving, the only resonant process is ∝ g 2Q 2 sbb , where a single photon from the drive and two phonons are converted into an excitation of the symmetric electrical mode. Using the rotating wave approximation, we thus neglect a term ∝ g 2Q 2 sb †b , as well as the dynamical contribution from the quadratic coupling between the asymmetric electrical and mechanical modes. The latter follows from the absence of resonance with the asymmetric electrical mode, which is satisfied whenever L 0 L. Finally, we have retained two linear terms ∝ δg 1Q 2 s (b +b † ) and ∝ g 1QaQs (b +b † ) in Eq. (1) that are responsible for the detrimental heating processes described in the introduction.
From the Hamiltonian, we derive the master equation for the mechanical system (see App. A for details): whereρ m is the membrane's density matrix, and L[Ô;ρ m ] = 2Ôρ mÔ † − Ô †Ô ,ρ m for any operatorÔ.
In deriving Eq.
(2) is valid in the weak coupling regime, where excitations transferred from the mechanical to the electrical subsystems cannot be transferred back coherently. This occurs when Γ Furthermore, we have assumed the resolved sideband regime γ ω m to write the rates as above and neglect the sideband at frequency ω s + 2ω m , associated with twophonon heating.
The relevant processes in the system are schematically represented in Fig. 1(c). The twophonon cooling ∝ Γ 2 arises from the quadratic term in Eq. (1), while "standard" heating and cooling of a system coupled to a Markovian reservoir is described by the rates Γ  1 . The first originates from the residual linear coupling δg 1 between the symmetric and the mechanical modes. The latter originates from cross-coupling between symmetric and anti-symmetric modes described by the three-body term ∝ g 1 in Eq. (1). In these processes, a single phonon can either be annihilated or created by interacting with a photon. The membrane is then both cooled down at a rate ∝ Γ  1 . The difference arises because the sideband at frequency ω s − ω m is enhanced by being closer to resonance compared to the sideband at frequency ω s + ω m . The asymmetric mode, however, is assumed to be far detuned ω a ω s , and therefore the associated heating and cooling rates Γ (b) 1 are the same. We note that both the desired two-photon cooling Γ 2 and the undesired heating processes Γ 1 , are proportional to |α| 2 and can be controlled by the strength of the drive. To understand whether the two-phonon cooling is domi- When the parameter λ r (λ m ) is larger than one, the non-Gaussian dynamics from the quadratic coupling prevails over the single phonon processes induced by the residual linear coupling δg 1 (multi-mode term proportional to g 1 ). Therefore, when λ 1, the two-phonon cooling is the dominant process and the membrane is effectively confined to the two lowest states, thus representing a qubit. This mechanical qubit's coherence time λ/Γ 2 corresponds to the timescale at which un-wanted processes act. We remark that the intrinsic mechanical damping is also detrimental for the two-phonon cooling, but its rate is not enhanced by the drive |α| 2 . While it is in principle possible to speed up all other processes to make it negligible, the small magnitude of the couplings and limitations on the driving power set requirements on the mechanical Q factor necessary to achieve such a situation (see Sec. 3 for details).

Experimental feasibility
A realistic estimate of λ can be derived from recent experiments [14,16,22,26,27,[58][59][60][61]. We consider a rectangular monolayer graphene membrane with area 1 × 0.3 µm 2 oscillating at ω m = (2π)80 MHz and suspended d 0 = 10 nm above a conducting plate, forming the capacitor [see Fig. 1(a)]. Following Refs. [48,62], we derive the ratio g 2 /g 1 = π 2 x zpm /(8d 0 ), where x zpm ≡ /(2mω m ) is the zero point motion's amplitude of a membrane of mass m. For the electrical circuit, we assume ω s = (2π)7 GHz, γ = (2π)150 kHz and δ ∈ [10 −3 , 10 −2 ], which leads to varying λ. To include the influence of both detrimental contributions we set Γ The number of photons in the cavity can be regulated to ensure that the system is in the weak coupling regime Γ 2 γ. This condition is equivalent (once Γ 2 is plugged in) to N mw 4γ 2 /g 2 2 = 9 × 10 10 . Here, N mw is the number of microwave photons in the cavity and we assumed g 2 ∼ (2π)1 Hz, corresponding to a stray capacitance C s = 100C 0 (see below). The associated power P up is calculated to be P up = ω s N mw (4ω 2 m + γ 2 )/γ = 448 mW. The input power P in must therefore be upper bounded by P in P up in order to remain in the weak coupling limit. We remark that, being the input drive off-resonant, only a very small fraction of the power is dissipated into the device. The corresponding intra-cavity dissipated power associated with P up found above is P up γ 2 /(4ω 2 m ) = 397 nW. Below, we estimate that for N mw ≤ 10 10 , corresponding to an intra-cavity dissipated (input) power of 44 nW (50 mW), a mechanical quality factor Q = ω m /γ m 10 5 is sufficient to be dominated by the driving-induced processes (see also App. E). Therefore, we set γ m = 0 in the following and explain how to generalize our results to non-negligible values of γ m at the end of this section.
To understand how the two-phonon cooling works, we simulate Eq. (2) for the parameters above and an initial thermal state ofn average phonons (for the evolution of an initial coherent state, see App. D). Fig. 2(a) shows the time dependence of the mechanical Fock states' amplitudes forn = 4 and λ = 2000. From an initial thermal distribution, highly populated states decay into lower ones until the ground and first excited states remain. Horizontal dashed lines indicate the expected populations of perfect twophonon cooling; (n + 1)/(2n + 1) for |0 and n/(2n+1) for |1 . In Fig. 2 As it is possible to see from Fig. 2(b), F scales as 1/λ, reflecting the compromise between the two-phonon cooling and the detrimental processes described above. Since a longer time is required for decaying into the ground and first excited states, lower fidelities are reached for higher n. However, the probability ∼ 25/(216λ) of the steady state being outside the {|0 , |1 } subsector [dashed black line in Fig. 2(b)] is independent of the initial state, highlighting that the two-phonon cooling is always the dominant process. The shadowed area corresponds to the range λ ∈ [34, 3400] found above, demonstrating that non-Gaussian operations are within reach of current technology.
In the remainder of this section, we describe experimental limitations originating from a parasitic capacitance C s , from extra heating induced by intracavity photons, and from the mechanical rate γ m . Starting with the first, the effect of C s is to limit the coupling strengths g 1 and g 2 , which are proportional to 1/(C 0 + C s ). A stray capacitance that is hundreds or even thousands of times bigger than C 0 is expected in realistic scenarios [22,61,[65][66][67][68]. Assuming C s = 100C 0 [65], the linear coupling for the example parameters considered above becomes g 1 ∼ (2π)7 kHz, and the quadratic g 2 ∼ (2π)1 Hz. Since the optomechanical couplings enter in λ as a ratio, stray capacitances do not have a direct effect on the realization of the two-phonon cooling. In the model considered so far, both a limited value of the Q factor and the presence of a stray capacitance can thus be compensated with a stronger drive.
In realistic examples, however, one cannot simply increase the driving strengths to arbitrary values to neglect the effects of the environment. Whether the system can tolerate high powers is a complex question which depends on the details of the experimental setup. As explained above, the photon flux that needs to be sent into the circuit to overcome the intrinsic heating depends on the mechanical quality factor Q and the absolute value of the couplings. With C s = 100C 0 [65] and Q ∈ [10 5 , 10 7 ], the required number of microwave photons N mw in the cavity can vary from 10 10 to 10 8 [48]. The question is then how much heating is induced by these photons. As an example, in Ref. [69] heating effects are observed starting from N mw 10 8 . We believe, however, that compared to that work the heating may be reduced for our system. The experimental results indicate that the heating is due to electric fluctuations in the main electrical modeΦ s . This heating should be reduced by the vanishing linear coupling δg 1 (due to symmetry) . We may now determine the contribution from γ m to λ. Since λ is defined as the ratio between the two-phonon cooling rate Γ 2 and the combined rate at which the system's state looses its coherence, it is possible to include γ m into λ by redefining (4) All conclusions drawn for λ above and in the following Sec. 4, hold upon the substitution λ → λ . Therefore, whether or not the mechanical damping γ m can be ignored depends on the ra- By increasing the flux of photons |α| 2 (and consequently N mw ), the numerator is enhanced, which reduces the mechanical quality factor Q required for neglecting the intrinsic damping γ m .
While several aspects must be considered for realistic experimental scenarios, the flexibility of the system allows overcoming imperfections such as stray capacitances and intrinsic heating. Two detrimental contributions to the desired system dynamics that are challenging to be estimated are mechanical dephasing and additional heating effects coming from the driving. As commented above, there are qualitative arguments suggesting that the latter does not jeopardize the twophonon cooling. Dephasing, on the other side, is expected to be less or at most comparable to the intrinsic decay rate γ m [70,71]. However, by its nature the two-phonon cooling is resilient against it, making dephasing detrimental only for the coherent manipulation of the mechanical qubit -see the following Sec. 4. Importantly, even in scenarios in which the dephasing may be the dominant decoherence rate in the system, our proposal allows for implementing qubit rotations (see below) and thereby spin echo techniques [72] that can mitigate its detrimental effects.

Quantum manipulation
Above, we considered the two-phonon cooling and the condition to reduce the membrane's infinite dimensional Hilbert space onto a two level system. Here, we introduce a protocol to perform arbitrary rotations on the Bloch sphere of the re-sulting mechanical state. We assume the system to be initialized in a suitable state and later read out using e.g., Gaussian operations. Initialization can be done via ground state cooling [12][13][14]73] (a technique that can be implemented in the parameter regimes considered in Sec. 3 using the linear coupling introduced below), or by collapsing the wave function via detection. The latter requires λ 1 [48], similarly to the two-phonon cooling.
We assume that the two-phonon cooling is constantly running and is dominant compared to Γ 1 , i.e., λ 1. Furthermore, as sketched in Fig. 3(a), we modify the mechanical Hamiltonian to include an asymmetric linear drive: H m →Ĥ m + (Ωb + Ω * b † )/2. This driving term can be obtained by introducing two additional input fields that are off-resonant with the previously studied system dynamics and whose frequency difference equals ω m (for additional informations, see App. B). By changing their strengths and phases, it is possible to control Ω, which determines the rotation on the Bloch sphere. We choose Ω such that Γ (r) These conditions, which can always be satisfied when λ 1, ensure that the dynamics is constrained to the two lowest energy levels (despite the linear drive coupling the system to the second excited state) and that the coherent drive is dominant compared to the undesired thermalization processes.
Although Eq.
(2) can be used to compute the exact time evolution of the system, we derive an approximate solution to grasp the essential physics. We truncate the mechanical Hilbert space to the ground (|0 ) and first two excited 1 }, |2 can be adiabatically eliminated [see Fig. 3(a)] and we geṫ whereρ ij is the {i, j} component of the density matrixρ m ,ρ x =ρ 10 +ρ 01 ,ρ y = i(ρ 10 −ρ 01 ), andρ 00 = 1 −ρ 11 . To characterize the rotation, we consider a π-pulse applied to an arbitrary state η|0 + β|1 . The drive Ω is activated for the time π/|Ω| required to switch |0 into |1 and vice versa. The quality of the evolution depends on the initial state and λ. By expanding the fidelity F in terms of λ, we find where Ω = 2 (Γ (r) This value of Ω reflects the compromise between exceeding the thermalization rates Γ 1 , and having a sufficiently low error rate |Ω| 2 /(4Γ 2 ) from driving level |2 .
In Fig. 3(b), we plot the fidelity of the πpulse against λ for Γ  a lower bound for the real fidelity.
As shown in Fig. 3(b), the system performs near perfect non-Gaussian operations, provided λ is sufficiently large. However, even for limited λ, it is possible to deterministically produce interesting non-classical states. By applying the π-pulse to a system initially cooled to its ground state, we find that λ 0.25 is sufficient to show negativity of the Wigner function W (in contrast to the phonon number QND detection which requires λ 50 to see clear non-classical features [48]). Fig. 3(c) demonstrates that highly negative values can be obtained for λ = 20. Despite it representing a mixture between |0 and |1 , the negativity of the Wigner function confirms the non-classical properties of the state for much lower values of λ. Specifically, for λ = 0.25 the overlap F of the resulting state with the desired single phonon state is ∼ 56%, which is sufficient for negativity (Ref. [74] shows that the Wigner function is negative for any state with more than 50% chance of being in the first excited state). For higher λ, W resembles more and more the one of a single phonon, and the overlap monotonically increases. As a reference, for λ = 10 3 we get F ∼ 95%.

Conclusions and outlook
We propose a protocol to realize non-Gaussian interactions in electromechanics. This can considerably enlarge the class of states which can be prepared in mechanical systems, and thereby the phenomena that can be studied. Our setup relies on two-phonon cooling of the mechanical subsystem, enabled by the quadratic coupling between the membrane and the symmetric electrical mode. Interestingly, the conditions under which the two-phonon cooling is dominant are the same as for performing QND detection of the phonon number [48]. As opposed to the QND measurement, however, the two-phonon cooling is frequency selective, addresses only a single mode, and makes all others dynamically irrelevant.
Even for modest values of λ 1, the twophonon cooling allows for preparing states with negative Wigner functions. For large values of λ, the mechanical system can in principle be used for quantum information processing. The twophonon cooling provides a non-Gaussian operation that, in combination with Gaussian interactions, allows for multi-mode complex operations [39,75]. For instance, a bilinear interaction of the kindb 1b † 2 +b 2b † 1 (subscripts refer to individual mechanical modes), in combination with the twophonon cooling, provides a √ SWAP gate [76]. Along with the qubit rotations considered here, this constitutes a universal gate set, thus allowing for creating arbitrary multi-mode mechanical states, restricted within the two lowest energy levels of all oscillators.
While there are different ways to realize non-Gaussian interactions in mechanical systems, these generally require the strong coupling regime g 1 γ [40,50], or resort to the interaction with additional quantum systems such as superconducting devices [43,47]. From one side, our scheme does not rely on additional quantum systems, making it simpler to realize and protecting it from additional detrimental effects coming from, e.g., these systems' thermal reservoirs. On the other side, the suppression of the linear coupling alongside the high tunability of electromechanical setups, allows for greatly lowering the constraints that must be satisfied to exploit non-Gaussian interactions in electromechanical sys-tems (for details, see App. E). For these reasons, we expect this work to greatly facilitate the development of future devices that can overcome the limitations posed by Gaussian interactions.

Code Availability
The code and nummerical simulations are available upon request.
[4] Y Lai, S Pirotta, G Urbinati, D A Derivation of the Master equation Q s symmetric chargeΦ s symmetric flux Q a asymmetric chargeΦ a asymmetric flux (2)] to describe the electromechanical circuit in Fig. 1(b). Here, we provide the details of its derivation. As discussed in Sec. 2, the full HamiltonianĤ of the system corresponds to the sum of four termsĤ =Ĥ m +Ĥ s +Ĥ a +Ĥ int , wherê The HamiltonianĤ describes the dynamics of three interacting subsystems: the mechanical membrane of frequency ω m , the symmetric electrical mode of frequency ω s ≡ [C 0 (L+2L 0 )] − 1 2 , and the asymmetric electrical mode of frequency ω a ≡ (C 0 L) − 1 2 . Throughout our work, we assume ω m ω s ω a . For completeness, in Eqs. (7) we include all terms, such as the noise contributions, which were omitted in Eq. (1). Specifically,V R0 ,V R1 andV R2 are the Johnson-Nyquist [77] noises associated with the resistor R 0 and the two resistors R, respectively [see Fig. 1].F m , which is associated to the mechanical damping γ m , describes a random force acting on the membrane as a white noise. The additional terms that appear in Eq. (7d) but are absent from Eq. (1) arise from the expansions of the two halves of the capacitor C[±x] and are the ones proportional to δg 1Q 2 a , g 2Q 2 a , and g 2Q 2 sb †b in Eq. (7d). We explain in the following that they are negligible in the considered parameter regime.
Assuming a strongly driven coherent fieldV in at the frequency ω in = ω s − 2ω m , it is advantageous to linearize the operators usingV where . Here A in is the amplitude of the driveV in with γ t = Z out /(L 0 + L/2), and .

(10)
For convenience the phase of A in is fixed such that α ∈ R. Furthermore,δ V in is the Johnson-Nyquist noise associated with the incoming field. To resemble a realistic experiment, we assume that the circuit is driven by a semi-infinite transmission line of impedance Z out as in Ref. [48]. The results in Sec. 2 are derived assuming matching impedances between the circuit and the transmission line: Z out = R 0 +R/2. For generality, we keep R 0 , R, and Z out as distinct parameters in the following.
Plugging Eqs. (8) into Eqs. (7), we obtain the linearized HamiltonianĤ lin which can be used to determine the Heisenberg equations for the system operators. Before doing so, let us comment on the underlined terms in Eq. (11). In the third and fourth rows, we collected all contributions that can be neglected, as they are not enhanced by the strong drive and/or are offresonant. The two terms in the fifth row can be absorbed into the remainders by a shift in the rest position of the membrane [1], and a redefinition of the mechanical frequency ω m . The equations of motions resulting from the dynamically relevant terms in Eq. (11) are thuṡ where we have included decays into Markovian reservoirs associated to each subsystem. Specifically, γ m is the decay rate into the mechanical reservoir, γ l = R/L is the decay rate into the reservoir of the asymmetric electrical mode, and γ = γ r + γ t is the decay rate into the reservoirs (associated with R 0 + R/2 and Z out , respectively) of the symmetric electrical mode. Here, γ r = (R 0 + R/2)/(L 0 + L/2) and γ t is defined above.
The system in Eqs. (12) combines nonlinearly coupled differential equations, and does not have a simple analytical solution. However, under reasonable assumptions, it is possible to select the leading contributions to the dynamics -which will lead us to the desired master equation -and to remove negligible contributions. To do so, it is convenient to switch to the frequency domain. For an operator O(t), we define its Fourier series asÔ where Ω k = 2πk/τ (k ∈ Z) are the allowed frequencies, and τ is a sufficiently long time interval. The Fourier coefficientsÔ[Ω k ] are then defined bŷ The differential equations in Eqs. (12) can now be Fourier transformed into polynomial equations by using the property (dÔ/dt) We can thus find relations forQ a andδ Q s : while Eq. (12a) forb becomes Here, we have defined Resorting to the Fourier transform does not in itself allow us to find solutions to the equations of motion in Eqs. (12). Indeed, operator products in the time domain become convolutions in the frequency domain, which can be taken care of only by considering infinite many terms. However, in the limit where the mechanical subsystem is characterized by a narrow bandwidth γ m , the frequency domain is useful for selecting the dominant dynamical terms [48]. In particular, it is reasonable to assume that the only relevant contributions in the sums in Eqs. (15) and (16) Here, all noises need to be evaluated at the frequency indicated in the associated square brackets. This last equation describes the membrane's dynamics, including the effects of the reservoirs of the electrical subsystems, but excluding the symmetric and the asymmetric electrical modes. From Eq. (17), it is possible to derive its corresponding effective HamiltonianĤ eff aŝ which is also independent from the symmetric and asymmetric electrical modes.

B Mechanical beam splitter interaction
To coherently manipulate the mechanical state, in Sec. 4 we have employed a beam splitter Hamiltonian of the formĤ bs = Ωb + Ω * b † 2 . (22) In this section, we explain how to experimentally implement such interaction by using two input fieldŝ V 1 andV 2 at frequencies ω 1 and ω 2 such that ω 1 − ω 2 = ω m . We assume that bothV 1 andV 2 are off-resonant with all other system dynamics, so that it does not interfere with the two-phonon cooling.
The Hamiltonian considered here is of the formĤ =Ĥ m +Ĥ o +Ĥ int , whereĤ m = ω mb †b is the mechanical term considered throughout this work, and Here, subscript "o" is used to indicate an electrical mode with intrinsic inductance and capacitance L o and C o , respectively, that is coupled to the membrane via the linear coefficientg 1 . For practical purposes, this electrical mode can be the symmetric mode, such that L o → L 0 + L/2, C o → 2C 0 and exploiting the inevitable fabrication imperfections we setg 1 → δg 1 . In this case, via the superposition principle, we can analyse the effects from the drivingsV 1 andV 2 independently from the other dynamics studied in Sec. 2 and in App. A, and include them afterwards in the effective HamiltonianĤ bs in Eq. (22).
Following the same steps as in App. A, we assume that the input fieldsV 1 andV 2 are strongly driven with averages given by = 1, 2). As a consequence, the average charge becomes where α i (i = 1, 2) is defined by Here, γ o is the decay rate of the mode "o" into the transmission line used to drive the cavity with the input fieldsV 1 andV 2 . We remark that, in principle, the drives considered here also affect the other electrical modes, such as the symmetric and asymmetric studied in the rest of this manuscript. In practice, these contributions can be made negligible by tuning the frequencies ω 1 and ω 2 .
With the purpose of deriving the beam splitter Hamiltonian in Eq. (22) we use the mean field approximation for the chargeQ o . By plugging Eq. (24) into Eq. (23b) we find that the mechanical subsystem is subjected to an effective interaction of the form In the interaction picture, for ω 1 − ω 2 = ω m the time dependence e i(ω 1 −ω 2 )t is compensated by the free evolution of the mechanical annihilation operatorb. Therefore, upon the substitution Ω =g 1 α 1 α * 2 /2, Eq. (25) is the beam splitter Hamiltonian. This means that from the membrane's point of view, the dynamics introduced by the two fieldsV 1 andV 2 driving the electrical mode "o" can be viewed as the desiredĤ bs used in the main text for the coherent manipulation.
We note that exploiting the fabrication imperfection to achieve a linear couplingg 1 → δg 1 is perfectly compatible with the general requirement of a negligible δg 1 to prevent linear heating. The undesired heating is proportional to the vacuum fluctuations and suppressed by being far off resonance. In contrast the linear coupling is a resonant process enhanced by the driving field making is much stronger for the same asymmetry parameter. A similar argument can be made for the linear cooling assumed for initialization of the system before the coherent drive, although in this case the enhancement only comes from being on resonance. The black line represents the sum j>1 j|ρ m |j , and the purple dotted-dashed represent P 0+1 = P 0 + P 1 . Dashed lines are the expected populations of a perfect two-phonon cooling process (see Sec. 3), which is achieved transiently when 1 Γ 2 t 10 2 . For longer times, detrimental linear processes kick in and drive the system into its steady state, as highlighted by the dotted lines. All curves are obtained by simulating Eq. (2); the parameters used are the same as in Fig. 2.

C Dynamical contributions from the residual linear couplings
In this section, we analyze the effects of residual linear couplings on the system's dynamics. We consider the same settings as in Secs. 2 and 3. In the interesting regime λ 1, the two-phonon cooling is dominant and the system promptly relaxes into its two lowest energy states. This can be seen in Fig. 2, from which it is evident that in a timescale t ∼ Γ −1 2 we obtain the desired outcome of a perfect two-phonon cooling process. Here, we investigate the membrane's dynamics afterwards, i.e., on a timescale t λ/Γ 2 in which detrimental heating effects disturb the mechanical two-level system's coherence.
The resulting dynamics on longer time scales are shown in Fig. 4. In this plot, we use the same parameters as in Fig. 2, but simulate the system for longer times to demonstrate the detrimental effects of the linear couplings. Since the two-phonon cooling is always the dominant process, it is possible to neglect all Fock states except the ground and first excited mechanical states. The effective dynamics is then characterized by the system of Eqs. (5) with Ω = 0. By setting the time derivatives to zero, the steady states are found to be and P 0 = 1 − P 1 . In case Γ (r) 1 , we obtain P 1 and P 2 to be 5/24 and 19/24, respectively, which are the dotted lines in Fig. 4. Full lines are obtained by numerically simulating Eq. (2), which yields the exact membrane's dynamics, and we include the purple dotted-dashed line, representing the total population P 0+1 = P 0 + P 1 of the two lowest mechanical states.
As it is possible to see in Fig. 4, as long as the two-phonon cooling is active, there is no leakeage outside the subspace of the ground and first excited states. This can be understood by looking at the purple dashed-dotted line that reaches its asymptote 1 for t Γ −1 2 . However, while the {|0 , |1 } subspace is always protected by the two-phonon cooling, the coherence within is still limited by detrimental linear effects. In fact, on a timescale t ∝ (Γ

D Time evolution of pure states
In Secs. 2 and 3 we considered the mechanical state before the two-phonon cooling to be a thermal state. This choice is motivated by it resembling realistic experimental scenarios. At the same time the density matrices of thermal states are diagonal, making the numerical computation simpler. This allowed us to push our simulations to include more states, increase the precision and to run for longer periods of time.
In this section, we investigate the opposite scenario. Namely, we consider an initial highly coherent mechanical state whose density matrix has a large off-diagonal component. First, we demonstrate that the steady state of a perfect two-phonon cooling process is such thatρ ij = 0 for all i, j = 0, 1, i.e., all entries of the mechanical density matrixρ m outside the lowest two levels' subspace are zero. Then, we consider the initial mechanical state and show that a perfect two-phonon cooling process applied to it results in a mixed state that retains large coherencesρ 01 andρ 10 between |0 and |1 .
To prove that a perfect two-phonon cooling process asymptotically (for t → ∞) eliminates all entrieŝ ρ ij (i, j = 0, 1) ofρ m , we start from the master equation in Eq. (2). By setting γ m = Γ This demonstrates that not only the diagonal, but also all the off-diagonal elements of the mechanical density matrix outside the two lowest energy levels subspace, asymptotically reach zero when the two phonon cooling is active.
The question is then which state within the {|0 , |1 } subspace the membrane reaches after a perfect two-phonon cooling process. This depends on the initial mechanical state. As we have seen in Sec. 3, a thermal state will be projected onto [(n + 1)|0 0| +n|1 1|] /(2n + 1), that has no off-diagonal elements, withn being the initial thermal occupation. If we use |Ψ k in Eq. (27) instead, it is possible to prove that the resulting two-level density matrixρ m has the following elementŝ ρ 01 =ρ 10 = 1 2k Qualitatively, Eq. (28a) follows from symmetry and normalization. Since the initial state |Ψ k in Eq. (27) has equal contributions from even and odd Fock states, the mechanical density matrix at the end of the cooling must have equal populations in its two lowest energy levels. Determining the off diagonal elementsρ 01 andρ 10 is more complicated. With Eq. (28b) derived by induction, the reason is that only fractions of the off-diagonal elements of the initial density are transferred intoρ 01 andρ 10 and the rest are lost. Specifically, these fractions depend on the decay rates of the two Fock states |i and |j involved in the decay of |i j|. More energetic states (with higher i and j) are associated with a larger loss. This difference in decay rates reduces the coherence between the states, i.e. the rates at which states decay provide "which-way" information about the state, which disturbs the coherence. This effect is enhanced with more decays and as a resultρ 01 andρ 10 in Eq. (28b) are monotonically decreasing functions of the parameter k.
Albeit this last fact, as proven by the asymptotic limit for k → ∞ in Eq. (28b), even for very large values of k the coherence of the initial state survives to a large degree. For k = 3 we find that ρ 01 =ρ 10 0.451, which we confirm by numerical simultion. In Fig. 5 we present a simulation of the whole system with |Ψ 3 [see Eq. (27)] as initial state and the two-phonon cooling being the dominant process. As it is possible to see in panel (a), the probabilities of being in each Fock state (full lines) are reduced to zero from 1/6, except for P 0 and P 1 that both reach approximately the value in Eq. (28a). This is further demonstrated by the density matrices depicted in panels (b) and (c), describing the with the two-phonon cooling active. We assume λ = 2000 and initialize the system in the superposition state |Ψ 3 = (|0 + |1 + |2 + |3 + |4 + |5 )/ √ 6. The time Γ 2 t = 2.5 is where P 0 P 1 . In (b) and (c), we present the density matrices at t = 0 and at Γ 2 t = 10 respectively. These panels demonstrate how the two-phonon cooling concentrates the mechanical state into the two lowest energy levels. Furthermore, starting from an initial coherent state, it prepares a coherent superposition in the {|0 , |1 } subspace as evidenced by the off-diagonal elements of the density matrix. All curves are derived by simulating Eq. (2); the parameters used are the same as in Figs. 2 and 4. system at the initial and final times, respectively. Specifically, the latter is compatible with the values ofρ 01 orρ 10 in Eq. (28b), in agreement with the above discussion.

E Quadratic coupling vs strong coupling regime
An important question to be addressed is whether the system reaches any strong coupling regime simultaneously with our non-Gaussian operations from the quadratic coupling. This would imply that the Master equation in Eq. (2) would not be valid. To address this concern, we would like to make three distinctions about relevant kinds of strong couplings. First, there is the (non-Gaussian) single photon strong coupling regime, characterized by g 1 γ (γ being the decay rate of the main electrical mode, which is assumed to be more damped than the mechanical one) [50]. Achieving this regime is known to be extremely challenging. We refer to [48] for demonstrating that this limit is far more demanding than λ 1 (also this regime is driving independent). Second, there is the (Gaussian) linear strong coupling regime, which can be obtained -if the coupling is resonant -when g 1 α γ [40]. In our settings, the linear coupling is suppressed both by symmetry and because the main electrical mode is far off resonance. As long as λ 1, the effect of the linear coupling is thus smaller than that of the quadratic, hence we would actually reach the quadratic strong coupling regime before the linear one. Finally, the quadratic strong coupling regime is characterized by g 2 α γ. To be dominated by the two-photon cooling as opposed to the linear heating, we need g 2 2 |α| 2 /γ γ m (2n m + 1). The combination of these two requirements again puts a requirement on the quality factor Q of the mechanical oscillator. The parameters considered in this work are Q ∈ [10 5 , 10 7 ], γ/(2π) = 150 kHz and γ m (2n m + 1) ≤ 1 kHz (see Sec. 3). Hence, there is ample margin to adjust the driving strength such that γ g 2 2 |α| 2 γ γ m (2n m + 1), as required for our master equation in Eq.
(2) to work. Therefore, we believe that the quadratic interaction considered by us is a very promising route to show non-Gaussian effects in electro-mechanics, and the analysis we provided is accurate in the assumed parameter regimes. This is also supported by numerical simulations of the composite system before tracing out the electrical symmetric and asymmetric modes [Eqs. (12) in App. A].