Synthesizing multi-phonon quantum superposition states using flux-mediated three-body interactions with superconducting qubits

TL;DR

This paper proposes a scheme to control a quantum mechanical resonator using superconducting qubits, enabling the preparation of non-classical states and entanglement, advancing quantum technology integration.

Contribution

It introduces a tunable three-body interaction scheme for quantum control of mechanical resonators with superconducting qubits, including ground-state cooling and state preparation.

Findings

Demonstrates ground-state cooling of the mechanical resonator

Achieves high-fidelity non-classical state preparation

Enables qubit-phonon entanglement with negative Wigner functions

Abstract

Massive mechanical resonators operating at the quantum scale can enable a large variety of applications in quantum technologies, as well as fundamental tests of quantum theory. Of crucial importance in that direction, is both their integrability into state-of-the-art quantum platforms as well as the ability to prepare them in generic quantum states using well-controlled high-fidelity operations. Here, we propose a scheme for controlling a radio-frequency mechanical resonator at the quantum scale using two superconducting transmon qubits that can be integrated on the same chip. Specifically, we consider two qubits coupled via a capacitor in parallel to a superconducting quantum interference device (SQUID), which has a suspended mechanical beam embedded in one of its arms. Following a theoretical analysis of the quantum system, we find that this configuration, in combination with an in-plane magnetic field, can give rise to a tuneable three-body interaction in the single-photon strong-coupling regime, while enabling suppression of the stray qubit-qubit coupling. Using state-of-the-art parameters and qubit operations at single-excitation levels, we numerically demonstrate the possibility of ground-state cooling as well as high-fidelity preparation of mechanical quantum states and qubit-phonon entanglement, i.e. states having negative Wigner functions and obeying non-classical correlations. Our work significantly extends the quantum control toolbox of radio-frequency mechanical resonators and may serve as a promising architecture for integrating such mechanical elements with transmon-based quantum processors.

Tables

Table 1. Table 1: Parameter set used in the numerical simulations.

Parameter	Value
${\omega }_{\text{m}}/\left(2\pi \right)$	10 MHz
${\omega }_{1,2}/\left(2\pi \right)$	$\sim$7 GHz
$g/\left(2\pi \right)$	$0.3$ MHz
$B$	10 mT
${\Phi }_{\text{b}}$	$0.495{\Phi }_0$
$l$	14.7 $\mu$m
$X_{\text{ZPF}}$	33 fm
${\beta }_0$	1
$C_{\text{c}}$	9.7 fF
$E_{\text{J},\Sigma }^{\text{c}}/E_{{\text{J}}_i}$	$\sim 10$
$E_{{\text{C}}_i}/h$	320 MHz
$n_{\text{th}}$ (T)	$\sim 20$ (10 mK)
$T_1,T_2$	30 $\mu$s
$Q_{\text{m}}$	${10}^6$

Equations

$$H=$$

$$+\frac{C_c}{C_1{C}_2}{Q}_1{Q}_2-E_{\text{J}}^{\text{c}}(X)\cos \left(\frac{{\varphi }_1-{\varphi }_2}{{\varphi }_0}\right),$$

$$E_{\text{J}}^{\text{c}}(X)=E_{\text{J},\Sigma }^{\text{c}}\left[c_{\text{J}}\cos (\pi {\Phi }_{\text{b}}/{\Phi }_0)-s_{\text{J}}\sin (\pi {\Phi }_{\text{b}}/{\Phi }_0)\alpha X\right],$$

$$\hat{H}={\hat{H}}_0+{\hat{H}}_{\text{int}},$$

$${\hat{H}}_0=\hslash {\omega }_{\text{m}}{\hat{b}}^{†}\hat{b}+\sum _{i=1}^2\hslash {\omega }_i{\hat{c}}_i^{†}{\hat{c}}_i-\frac{E_{{\text{C}}_i}}{2}{\hat{c}}_i^{†}{\hat{c}}_i^{†}{\hat{c}}_i{\hat{c}}_i,$$

$${\hat{H}}_{\text{int}}=\hslash g({\hat{c}}_1^{†}{\hat{c}}_2+{\hat{c}}_1{\hat{c}}_2^{†})(\hat{b}+{\hat{b}}^{†})$$

$$-\sum _{i=1}^2\hslash g_i{\hat{c}}_i^{†}{\hat{c}}_i(\hat{b}+{\hat{b}}^{†})+{\hat{H}}_{\text{12}}^{{}^{\prime }},$$

$$g=\frac{\alpha \sqrt{Z_1{Z}_2}}{2{\varphi }_0^2}{s}_{\text{J}}{E}_{\text{J},\Sigma }^{\text{c}}\sin (\pi {\Phi }_{\text{b}}/{\Phi }_0)X_{\text{ZPF}},$$

$${\hat{H}}_{\text{12}}^{{}^{\prime }}=\hslash (J_{\text{C}}-J_{\text{L}})({\hat{c}}_1^{†}{\hat{c}}_2+{\hat{c}}_1{\hat{c}}_2^{†})+\hslash V{\hat{c}}_1^{†}{\hat{c}}_1{\hat{c}}_2^{†}{\hat{c}}_2,$$

$$U(t)|0_1{0}_{\text{m}}{1}_2⟩=\cos (gt)|0_1{0}_{\text{m}}{1}_2⟩-i\sin (gt)|1_1{1}_{\text{m}}{0}_2⟩.$$

$$|\psi {⟩}_{\text{f}}=\frac{1}{\sqrt{2}}\left(|0_1{0}_{\text{m}}{1}_2⟩-i|1_1{1}_{\text{m}}{0}_2⟩\right),$$

$$|\psi {⟩}_{\text{f}}=\left(\frac{1}{\sqrt{2}}(|0_1{0}_{\text{m}}⟩+|1_1{1}_{\text{m}}⟩)\right)|0_2⟩,$$

$$\dot{\rho }=\frac{i}{\hslash }[\rho ,\hat{H}]$$

$$+\sum _{i=1}^2\frac{1}{T_1}\mathcal{L}[{\hat{c}}_i]\rho +\frac{1}{T_2}\mathcal{L}[{\hat{c}}_i^{†}{\hat{c}}_i]\rho ,$$

$$\mathcal{L}=$$

$$H=\frac{P^2}{2m}+\frac{m{\omega }_{\text{m}}^2{X}^2}{2}+\sum _{i=1}^2\left[\frac{Q_i^2}{2{\tilde{C}}_i}-E_{{\text{J}}_i}\cos \left(\frac{{\varphi }_i}{{\varphi }_0}\right)\right]+\frac{C_c}{C_1{C}_2}{Q}_1{Q}_2-E_{\text{J}}^{\text{c}}\cos \left(\frac{{\varphi }_1-{\varphi }_2}{{\varphi }_0}\right),$$

$$E_{\text{J}}^{\text{c}}=E_{\text{J},\Sigma }^{\text{c}}|\cos (\pi {\Phi }_{\text{b}}/{\Phi }_0+\alpha X)|,$$

$$E_{\text{J}}^{\text{c}}\simeq E_{\text{J},\Sigma }^{\text{c}}\left[\cos (\pi {\Phi }_{\text{b}}/{\Phi }_0)(1-{\alpha }^2{X}^2)-\sin (\pi {\Phi }_{\text{b}}/{\Phi }_0)\alpha X\right],$$

$$E_{\text{J}}^{\text{c}}({\Phi }_{\text{b}},X)=E_{\text{J},\Sigma }^{\text{c}}[{\cos }^2(\pi {\Phi }_{\text{b}}/{\Phi }_0+\alpha X)+a_{\text{J}}^2{\sin }^2(\pi {\Phi }_{\text{b}}/{\Phi }_0+\alpha X){]}^{1/2},$$

$${\cos }^2(\pi {\Phi }_{\text{b}}/{\Phi }_0+\alpha X)\simeq {\cos }^2(\pi {\Phi }_{\text{b}}/{\Phi }_0)-2\alpha X\cos (\pi {\Phi }_{\text{b}}/{\Phi }_0)\sin (\pi {\Phi }_{\text{b}}/{\Phi }_0),$$

$${\sin }^2(\pi {\Phi }_{\text{b}}/{\Phi }_0+\alpha X)\simeq {\sin }^2(\pi {\Phi }_{\text{b}}/{\Phi }_0)+2\alpha X\cos (\pi {\Phi }_{\text{b}}/{\Phi }_0)\sin (\pi {\Phi }_{\text{b}}/{\Phi }_0).$$

$$E_{\text{J}}^{\text{c}}=E_{\text{J},\Sigma }^{\text{c}}[1+a_{\text{J}}^2{\tan }^2(\pi {\Phi }_{\text{b}}/{\Phi }_0){]}^{1/2}\cos (\pi {\Phi }_{\text{b}}/{\Phi }_0){\left[1-\alpha X\frac{(1-a_{\text{J}}^2)\tan (\pi {\Phi }_{\text{b}}/{\Phi }_0)}{1+a_{\text{J}}^2{\tan }^2(\pi {\Phi }_{\text{b}}/{\Phi }_0)}\right]}^{1/2},$$

$$E_{\text{J}}^{\text{c}}\simeq E_{\text{J},\Sigma }^{\text{c}}\left[c_{\text{J}}\cos (\pi {\Phi }_{\text{b}}/{\Phi }_0)-s_{\text{J}}\sin (\pi {\Phi }_{\text{b}}/{\Phi }_0)\alpha X\right]+\mathcal{O}[X^2],$$

$$\mathcal{O}[X^2]=E_{\text{J},\Sigma }^{\text{c}}\frac{s_{\text{J}}{\sin }^2(\pi {\Phi }_{\text{b}}/{\Phi }_0)}{2c_{\text{J}}\cos (\pi {\Phi }_{\text{b}}/{\Phi }_0)}{\alpha }^2{X}^2,$$

$$H_{\text{int}}^{\text{(flux)}}=$$

$$H_{\text{3-body}}=\alpha E_{\text{J},\Sigma }^{\text{c}}{s}_{\text{J}}\sin (\pi {\Phi }_{\text{b}}/{\Phi }_0)\frac{{\varphi }_1{\varphi }_2}{{\varphi }_0^2}X,$$

$$H_{\text{RP}}=-\alpha E_{\text{J},\Sigma }^{\text{c}}{s}_{\text{J}}\sin (\pi {\Phi }_{\text{b}}/{\Phi }_0)\frac{({\varphi }_1^2+{\varphi }_2^2)}{2{\varphi }_0^2}X,$$

$$H^{\{{\varphi }^2\}}=-E_{\text{J},\Sigma }^{\text{c}}{c}_{\text{J}}\cos (\pi {\Phi }_{\text{b}}/{\Phi }_0)\frac{{\varphi }_1{\varphi }_2}{{\varphi }_0^2},$$

$$H^{\{{\varphi }^4\}}=E_{\text{J},\Sigma }^{\text{c}}{c}_{\text{J}}\cos (\pi {\Phi }_{\text{b}}/{\Phi }_0)\left(-\frac{{\varphi }_1^2{\varphi }_2^2}{4{\varphi }_0^4}+\frac{{\varphi }_1^3{\varphi }_2+{\varphi }_1{\varphi }_2^3}{6{\varphi }_0^4}\right),$$

$$H^{\{{\varphi }^{4}X\}}=\alpha E_{\text{J},\Sigma }^{\text{c}}{s}_{\text{J}}\sin (\pi {\Phi }_{\text{b}}/{\Phi }_0)\frac{({\varphi }_1-{\varphi }_2{)}^4}{24{\varphi }_0^4}X,$$

$$H^{\{{\varphi }^2{X}^2\}}=-{\alpha }^2{E}_{\text{J},\Sigma }^{\text{c}}\frac{s_{\text{J}}{\sin }^2(\pi {\Phi }_{\text{b}}/{\Phi }_0)}{2c_{\text{J}}\cos (\pi {\Phi }_{\text{b}}/{\Phi }_0)}\frac{{\left({\varphi }_1-{\varphi }_2\right)}^2}{2{\varphi }_0^2}{X}^2,$$

$$H^{\{{\varphi }^4{X}^2\}}={\alpha }^2{E}_{\text{J},\Sigma }^{\text{c}}\frac{s_{\text{J}}{\sin }^2(\pi {\Phi }_{\text{b}}/{\Phi }_0)}{2c_{\text{J}}\cos (\pi {\Phi }_{\text{b}}/{\Phi }_0)}\frac{({\varphi }_1-{\varphi }_2{)}^4}{24{\varphi }_0^4}{X}^2,$$

$$E_{\text{J},\Sigma }^{\text{c}}{c}_{\text{J}}\cos (\pi {\Phi }_{\text{b}}/{\Phi }_0)\left(\frac{{\varphi }_1^2+{\varphi }_2^2}{2{\varphi }_0^2}-\frac{{\varphi }_1^4+{\varphi }_2^4}{24{\varphi }_0^4}\right),$$

Full text

I Introduction

The ability to control massive mechanical objects at the quantum level constitutes a very interesting task for many technological applications, ranging from microwave-to-optical conversion to quantum memories, as well as fundamental studies regarding the quantum-classical divide penrose1996gravity ; marshall2003towards ; kleckner2008creating ; blencowe2011light ; blencowe2007quantum ; stannigel2012optomechanical ; bochmann2013nanomechanical ; andrews2014bidirectional ; metcalfe2014applications ; aspelmeyer2014cavity . The rapid development of cavity optomechanics over the last decade has enabled the exploration of mechanical resonators in regimes where quantum effects become prominent. One approach relies on resonant coupling of acoustic phonons to microwave excitations via piezoelectric materials oconnell2010quantum ; chu2018creation , however, the amplitude of the lattice vibrations in these systems is very small, limiting their applications. On the other hand, in typical opto- and electromechanical setups, low-frequency mechanical resonators are controlled via parametric coupling to an optical or microwave cavity chan2011laser ; teufel2011sideband ; wollman2015quantum ; riedinger2018remote ; ockeloen2018stabilized . The mechanical elements in these systems are usually realised with metal drumheads or beams, which are characterised by large quality factors and large mass, making them also prime candidates for experimental tests of gravity-induced wavefunction collapse theories penrose1996gravity ; marshall2003towards ; kleckner2008creating ; blencowe2011light .

Following the pioneering work on acoustic resonators oconnell2010quantum , recently, quantum superpositions of the ground and first excited state were for the first time generated in a parametrically coupled mechanical resonator reed2017faithful . In this approach, a superconducting qubit is used to create the excitation, which decays into a microwave resonator on a different chip and subsequently transferred to the mechanical element via an effective linearised interaction. It is recognised, however, that this method has severe limitations as it relies on strong driving, which is challenging with qubits, and suffers from unavoidable losses during the state transfer between different chips, limiting the fidelity of the prepared state viennot2018phonon . A different scheme, implemented in the optical domain, uses entanglement and post-selective measurements to generate single-photon states Hong2017handbury , although the non-deterministic nature of the protocol in combination with low count rates limits the types of states that can be prepared.

A promising route towards high-fidelity mechanical quantum control is the ability to operate in the single-photon strong-coupling regime, where interaction times are faster than dissipation processes, which however still remains an experimental challenge for far-detuned parametrically coupled mechanical resonators pirkkalainen2015cavity ; viennot2018phonon . Operating in this regime is predicted to give rise to non-classical photon correlations rabl2011photon and non-Gaussian states nunnenkamp2011single ; nation2013nonclassical , as well as macroscopic mechanical superpositions liao2016macroscopic . Moreover, using the radiation-pressure coupling to qubits can lead to the creation of mechanical Schrödinger cat states armour2002entanglement ; mehdi2015quantum . Generating Fock states in this regime could also be enabled using an additional microwave resonator to create an effective tripartite coupling, as predicted in Ref. mehdi2015quantum , although this approach is limited to low state preparation fidelities mainly due to the limitations of capacitive coupling and stray qubit-resonator interaction. In spite of the experimental and theoretical advances in the field, high-fidelity quantum state preparation of mechanical systems appears to be limited to a small class of engineerable states which are, to a large degree, architecture-dependent.

Here, we analyse a new scheme for synthesizing generic mechanical states by employing tuneable three-body interactions between two superconducting qubits and a mechanical resonator in the single-photon strong-coupling regime. The coupling relies on embedding a suspended micrometer-long beam in one of the arms of a superconducting quantum interference device (SQUID) in combination with an externally applied magnetic field nation2008quantum ; nation2016ultrastrong ; shevchuk2017strong ; rodrigues2019coupling . We find that, by connecting two superconducting transmon qubits koch2007charge directly via this mechanical SQUID, a tuneable three-body interaction arises as the qubit-qubit flux-mediated coupling is modulated by the mechanical displacement. Importantly, the detrimental exchange-type interaction between the two qubits can be suppressed by adding a capacitor in parallel to the SQUID, as realised in Ref. kounalakis2018tuneable . Using state-of-the-art parameters, reported in recent experiments rodrigues2019coupling , we numerically demonstrate the possibility of high-fidelity coherent quantum control of the beam, from ground-state cooling to fast and high-fidelity preparation of mechanical Fock states, as well as maximally entangled states of qubits and phonons. Finally, we devise a protocol consisting of qubit flux-pulsing and post-selective measurements for synthesizing multi-phonon superposition states, extending the quantum control toolbox and the plurality of engineerable quantum states in radio-frequency mechanical resonators.

II Results

II.1 Motion-dependent qubit-qubit interaction

The proposed circuit, shown in Fig. 1a, comprises two transmon qubits coupled directly via a superconducting quantum interference device (SQUID) shunted by a capacitor, $C_{\text{c}}$. This tuneable coupling scheme has recently been realised in circuit QED setups using transmons kounalakis2018tuneable and LC resonators collodo2018observation . The coupling is controlled by tuning the Josephson energy of the SQUID, $E_{\text{J}}^{\text{c}}(\Phi_{\text{b}})=E_{\text{J},\Sigma}^{\text{c}}|\cos(\pi\Phi_{\text{b}}/\Phi_{0})|$, with an out-of-plane flux bias, $\Phi_{\text{b}}$, where $\Phi_{0}=h/2e$ is the magnetic flux quantum and $E_{\text{J},\Sigma}^{\text{c}}$ is the sum of the two junction Josephson energies in the SQUID. The mechanical part of the circuit consists of a beam of length $l$ that is embedded in one of the arms of the SQUID loop and can oscillate out of plane. By applying an in-plane external magnetic field $B$, the loop can pick up an additional flux $\beta_{0}BlX$ due to the beam displacement $X$, resulting in a flux-tuneable and motion-dependent Josephson energy $E_{\text{J}}^{\text{c}}=E_{\text{J},\Sigma}^{\text{c}}|\cos(\pi\Phi_{\text{b}}/\Phi_{0}+\alpha X)|$, where $\alpha=\pi\beta_{0}Bl/\Phi_{0}$ and $\beta_{0}$ is a geometric factor depending on the mode shape (for the first mechanical mode, considered here, $\beta_{0}\simeq 1)$ etaki2008motion ; shevchuk2017strong ; rodrigues2019coupling .

The Hamiltonian describing the circuit in Fig. 1a is

[TABLE]

where $\{X,~{}P\}$ and $\{\phi_{i},~{}Q_{i}\}$ are conjugate variable pairs describing the mechanical and the electrical degrees of freedom at circuit node $i$, and $\phi_{0}=\Phi_{0}/(2\pi)$ is the reduced flux quantum. The first four terms describe the uncoupled system of the mechanical resonator and two transmon qubits, where $m,~{}\omega_{\text{m}}$ denote the mass and frequency of the beam and $E_{\text{J}_{i}}$, $\widetilde{C}_{i}$ represent the Josephson energy and loaded capacitance of each transmon, respectively (see Supplementary Sec. S1). The last two terms describe the qubit-qubit interaction via their charge and flux degrees of freedom, $Q_{i}$ and $\phi_{i}$.

The core of this proposal relies on the fact that the dynamical displacement of the beam results in a modulation of the superconducting current through the coupling SQUID and its Josephson energy, which mediates the qubit-qubit interaction. Taking into account the finite asymmetry of the SQUID loop, which is present in any realistic scenario, this results in a motion-dependent Josephson energy,

[TABLE]

where $c_{\text{J}}=\sqrt{1+a_{\text{J}}^{2}|\tan{(\pi\Phi_{\text{b}}/\Phi_{0})}|}$ and $s_{\text{J}}=(1-a_{\text{J}}^{2})/c_{\text{J}}$ are correction factors due to the SQUID asymmetry $a_{\text{J}}$. Note that the above expression is valid for $\pi\Phi_{\text{b}}/\Phi_{0}\gg\alpha X$ and relies on the assumption that $\alpha X\ll 1$ (for a full derivation see Supplementary Sec. S1). For the parameters considered here, which are compatible with values reported in recent experiments using micrometer-long Al beams and sub-Tesla magnetic fields rodrigues2019coupling , we have $\alpha X\sim 10^{-5}-10^{-6}$, therefore this is a valid assumption.

II.2 Electromechanical system dynamics

The Hamiltonian of the system can be expressed in second quantisation form, as

[TABLE]

where $\hat{b}^{(\dagger)}$ and $\hat{c}_{i}^{(\dagger)}$ are bosonic operators describing the annihilation (creation) of phonons and qubit excitations, respectively (see Supplementary Sec. S2). The effective electromechanical frequencies are $\omega_{\text{m}}$ and $\omega_{i}=\frac{1}{\hbar}\left(\sqrt{8\widetilde{E}_{\text{J}_{i}}E_{\text{C}_{i}}}-E_{\text{C}_{i}}\right)$, where $\widetilde{E}_{\text{J}_{i}}$ is the modified transmon Josephson energy due to the coupler. The full quantum mechanical treatment of the circuit, including higher-order nonlinear interaction terms, is presented in the Supplementary Material.

The first term in Eq. (5) describes a three-body interaction involving hopping of qubit excitations together with mechanical displacements of the beam. The coupling strength is given by

[TABLE]

where $Z_{i}=\frac{\hbar}{e^{2}}\sqrt{E_{\text{C}_{i}}/2\widetilde{E}_{\text{J}_{i}}}$ denote the transmon impedances, and $X_{\text{ZPF}}~{}=~{}\sqrt{{\hbar}/{(2m\omega_{\text{m}})}}$ is the zero-point motion of the mechanical resonator. The next interaction term describes the radiation-pressure coupling of each qubit with the beam at a rate $g_{1(2)}~{}=~{}g\sqrt{Z_{1(2)}/Z_{2(1)}}$. The electromechanical coupling strengths are plotted in Fig. 1b as a function of the out-of-plane flux bias $\Phi_{b}$. Close to a half-integer flux quantum, $E_{\text{J}}^{\text{c}}(X)$ is maximally susceptible to mechanical motion which maximises the electromechanical coupling. Note that all couplings become zero exactly at a half-integer flux quantum due to the finite SQUID asymmetry $a_{\text{J}}$. In our calculations an asymmetry of $0.01$ is included, reflecting a $2\%$ spread in junction fabrication targeting.

The last term in Eq. (5) describes the qubit-qubit interaction

[TABLE]

where $J_{\text{C}}$, $J_{\text{L}}$ are exchange-type coupling strengths arising from the coupling capacitor and SQUID, respectively, and $V$ is the cross-Kerr coupling strength, which is minimised at $\Phi_{b}\simeq\Phi_{0}/2$  (see Supplementary Eq. S(33)). A significant advantage of this architecture for realising tripartite interactions, compared to relying exclusively on capacitive coupling mehdi2015quantum , is that any stray linear coupling between the qubits can be suppressed with the right choice of coupling capacitance $C_{\text{c}}$ kounalakis2018tuneable ; collodo2018observation . This makes the three-body interaction dominant and ensures the ability to manipulate the state of each qubit individually by local driving, which is crucial for the state engineering protocols discussed below.

II.3 Ground-state cooling

Mechanical resonators realised using vibrating beams and drumheads lie in the radio-frequency regime ($\sim 10$ MHz), where thermal fluctuations are dominant even at millikelvin temperatures, achieved with conventional cryogenic techniques. An essential element of control is, therefore, the ability to cool these systems to their quantum ground state before manipulating them further. Typically, electromechanical experiments employ a “cold” microwave cavity ($\sim 10$ GHz), which is coupled to the mechanical element via an effective linearised interaction in the many-photon regime, and cooling is enabled by red-sideband driving teufel2011sideband . However, in a system comprised of qubits, cooling the mechanical resonator via sideband driving can be a challenging task, requiring multiple tones and eventually limited by the critical number of photons in the Josephson junction gely2019observation ; lescanne2019observing .

Here, we show that it is possible to overcome the challenges of sideband cooling with qubits by employing a time-domain protocol to cool the mechanical resonator to its quantum ground state, using the three-body interaction. The scheme, depicted in Fig. 2a, consists of a sequence of single-qubit operations which, combined with the tripartite interaction, enable the transfer of the thermal phonons to the environment in a stroboscopic fashion, as described below. At first, we bring one qubit ($q_{1}$) to its excited state and then tune its frequency such that $\omega_{1}=\omega_{2}-\omega_{\text{m}}$. Since $g\ll\omega_{\text{m}}$ the interaction $(\hat{c}_{1}\hat{b}\hat{c}_{2}^{\dagger}+\text{H.c.})$ is resonant at this condition, such that a phonon combined with the excitation in $q_{1}$ can be transferred to the other qubit ($q_{2}$) after variable time $\Delta t_{\text{cool}}$. The cycle is then completed by reinitialising both qubits using active reset protocols riste2012feedback ; campagne2013persistent ; magnard2018fast . A similar scheme has also been realised recently for detecting microwave photons in a superconducting resonator, using tripartite interactions with a transmon qubit and a dissipative mode lescanne2019detecting .

In Fig. 2b we plot the average number of phonons and qubit excitations as a function of time after $\sim 100$ cooling cycles, for a mechanical beam oscillating at $\omega_{\text{m}}/(2\pi)=10$ MHz and a tripartite coupling strength $g/(2\pi)\simeq 0.3$ MHz. At the end of the protocol the mechanical resonator is cooled down to the ground state with a phonon occupancy of 0.05, assuming an environment temperature of $T=10$ mK ($n_{\text{th}}\simeq 20$). The qubit ($\omega_{1}/(2\pi)=7$ GHz) is excited with a 200 ns Gaussian pulse, while the reset and cooling times are set to $\Delta t_{\text{reset}}=\Delta t_{\text{cool}}=200$ ns. The cooling time would be best optimised by choosing a short time for large thermal occupation and then increase the time as the resonator cools, since phonons swap faster for higher occupations. For simplicity, here, we considered a fixed time and found that for these parameters, a 200 ns time is sufficient. Furthermore, we have assumed an in-plane magnetic field $B=10$ mT which is well-below the critical field for thin Al beams meservey1971properties and does not compromise the qubit coherence schneider2019transmon . All system parameters used in the simulations are listed in Table 1.

II.4 Mechanical Fock states and qubit-phonon entanglement

Following ground-state preparation, we present a protocol which employs the tripartite coupling to deterministically generate mechanical Fock states and maximally entangled states. As schematically depicted in the inset of Fig. 3a, it consists of preparing $q_{2}$ in the excited state and tuning it to $\omega_{2}=\omega_{1}+\omega_{\text{m}}$, such that the interaction $(\hat{c}_{1}^{\dagger}\hat{b}^{\dagger}\hat{c}_{2}+\text{H.c.})$ is resonant. Ideally, assuming unitary evolution $U(t)$ under the tripartite interaction, the state would evolve as

[TABLE]

The evolution of the average number of phonons and qubit excitations, starting from an attainable mechanical state of 0.05 thermal phonons and including system dissipation (see Methods) is plotted in Fig. 3a. We consider the full interaction Hamiltonian presented in Eq. (5), including next-to-leading order non-linear correction (see Supplementary Sec. S2), for the simulation parameters shown in Table 1. In Fig. 3b we plot the Wigner function of the final state in the mechanical resonator, after time $t=\pi/(2g)$, revealing a single-phonon Fock state with $97\%$ ($99\%$) fidelity, starting from an attainable (ideal) ground state. Higher phonon states could also be prepared by resetting and repeating the protocol with modified transfer times $t_{n}=n_{\text{m}}^{-1/2}\pi/(2g)$.

The quantum state preparation scheme described above, can also be used to generate bipartite and tripartite maximally entangled states between the mechanical resonator and the qubits. In particular, in the middle of the above protocol, at $t=\pi/(4g)$, the system is in a Greenberger–Horne–Zeilinger (GHZ) state

[TABLE]

with $98~{}\%$ fidelity, with the corresponding density matrix shown in Fig. 3c. Such states are particularly interesting for applications in quantum information cleve1997substituting ; bruss1998optimal and fundamental tests of quantum theory greenberger1989going . Using the same protocol for $q_{2}$ in a superposition state $\frac{1}{\sqrt{2}}(|0_{2}\rangle+|1_{2}\rangle)$, the Bell state

[TABLE]

is generated after time $\pi/(2g)$ with $96\%$ fidelity ($98\%$ for ideal ground state), as depicted in Fig. 3d. The prepared state is a maximally entangled pair of a phonon and a qubit excitation, which could be utilised as a testbed for checking the validity of quantum mechanics at macroscopic scales without requiring tomography of the mechanical state hofer2016proposal ; vivoli2016proposal . Such states might also be suitable for integrating transmon qubits into other platforms such as spins, cold atoms, or even optical photons, via the mechanical resonator rabl2010quantum ; camerer2011realization ; arcizet2011single ; hill2012coherent . They could also provide possibilities for entangling the mechanical resonator with other physical systems via the transmon.

II.5 Multi-phonon quantum superpositions

We now extend the protocol described above to create multi-phonon quantum superposition states in the mechanical resonator, simply by flux-pulsing the qubits. In the protocols discussed previously, the qubit frequencies are tuned at $\omega_{2}=\omega_{1}+\omega_{\text{m}}$ such that the states $|0_{1}n_{\text{m}}1_{2}\rangle$ and $|1_{1}(n+1)_{\text{m}}0_{2}\rangle$ are coupled. However, when tuned at $\omega_{2}=\omega_{1}-\omega_{\text{m}}$ the interaction term $(\hat{c}_{1}\hat{b}^{\dagger}\hat{c}_{2}^{\dagger}+\text{H.c.})$ becomes resonant, which couples $|0_{1}(n+1)_{\text{m}}1_{2}\rangle$ and $|1_{1}n_{\text{m}}0_{2}\rangle$. Therefore, by interchanging the qubit frequencies with flux-tuning pulses during each cycle it could be possible to create higher phonon Fock states and multi-phonon quantum superposition states, as depicted schematically in Fig. 4a.

As a proof-of-concept, using the same simulation parameters as above (Table 1), we demonstrate the creation of superposition states $|\psi_{\text{m}}\rangle^{{}^{\prime}}=\frac{1}{2}(|0\rangle+\sqrt{2}|2\rangle+|4\rangle)$ and $|\psi_{\text{m}}\rangle^{{}^{\prime\prime}}=\frac{1}{\sqrt{2}}(|0\rangle+|4\rangle)$, after exciting qubit 2 and applying three flux pulses that interchange the qubit frequencies at variable times $t_{1},t_{2}$ and $t_{3}$. Figs. 4b, d show the evolution of the qubit and resonator occupancy, starting from an attainable mechanical state of 0.05 phonons. The dashed lines indicate the times that a flux-tuning pulse is applied. The corresponding Wigner functions at the end of each protocol, following post-selection on $|0_{1}1_{2}\rangle$, are shown in Figs. 4c, e, with preparation fidelities $98~{}\%$ and $97~{}\%$, respectively. After preparation, the states evolve naturally as $U(t)|\psi_{\text{m}}\rangle^{{}^{\prime}}=(\alpha|0\rangle+\beta e^{-in\omega_{\text{m}}t}|n\rangle+\gamma e^{-im\omega_{\text{m}}t}|m\rangle)$ including dissipation, which is however not a limiting factor because of the long lifetimes of these mechanical resonators rodrigues2019coupling . Readout of the prepared states including Wigner tomography could be performed using similar techniques to the ones developed in Ref. hofheinz2009synthesizing .

The scheme described above enables the generation of interesting classes of multi-phonon superposition states, such as the ones shown in Fig. 4, requiring only flux-tuning pulses and a projective measurement at the end. As we show in Supplementary Sec. S3, by including a projective measurement after each step of the protocol, it is possible to generate states with arbitrary phonon number probability distributions, although constrained in the relative phases of the superpositions. Furthermore, we find that quantum superpositions with arbitrary complex coefficients can also be generated with this platform by additionally employing the qubit-qubit interaction in a controllable fashion to perform exchange-type and C-Phase gates between the two qubits (see Supplementary Sec. S3). This would enable the creation of truly arbitrary states, similar to those produced in resonantly coupled qubit-resonator systems law1996arbitrary ; hofheinz2009synthesizing , with the trade-off of increased complexity in the protocol. Alternatively, using the radiation-pressure coupling with one qubit in combination with a sequence of driving pulses, could enable the creation of mechanical Schrödinger cat states as discussed in Ref. mehdi2015quantum .

III Discussion

A reconciliation of quantum mechanics and general relativity remains elusive at a theoretical level, however, there exist several proposals for testing the quantum-classical boundary with mechanical resonators offering an ideal testbed. More specifically, it has been theorised that a massive object in quantum superposition results in two coexisting space-time geometries, leading to issues with the unitary evolution, which eventually causes it to collapse penrose1996gravity ; kleckner2008creating . Importantly, this relies on the zero-point motion $X_{\text{ZPF}}$ being much larger than the approximate size of the nucleus ($\sim 1$ fm), which is the case in our system ($X_{\text{ZPF}}=33$ fm). The collapse timescale $t_{\text{G}}$ is inversely proportional to the mass of the object, resulting in $t_{\text{G}}\sim 1-10$ s for the parameters considered here ($m\sim 1$ pg), therefore the resonator coherence time should be larger than that. Recent advances in strain engineering techniques can enable the enhancement of beam quality factors up to $Q_{\text{m}}\sim 10^{9}$ ghadimi2018elastic , leading to relaxation times of hundreds of seconds, which would be sufficient for observing gravitational effects. Moreover, the ability to prepare a large variety of superposition states could offer an additional tool in testing such theories. The proposed architecture provides a very versatile platform in this regard, enabling not only generic quantum state preparation, but also with high fidelity, which has so far been a very challenging task.

Our approach combines the best of both worlds of two very versatile systems, namely the exquisite level of quantum control of qubits in circuit quantum electrodynamics blais2004cavity , with the long lifetimes and flexibility of mechanical elements in coupling to electromagnetic radiation. The high-fidelity generation of hybrid entangled states of phonons and qubit excitations, which have no classical analogue, may provide alternative routes for testing the limits of quantum theory at macroscopic scales greenberger1989going ; hofer2016proposal ; vivoli2016proposal . Additionally, such states are of particular importance in enabling quantum technologies with hybrid quantum systems, from quantum simulation to quantum computing and communication kurizki2015quantum , and could also be used for coupling qubit excitations with other systems such as optical photons, cold atoms, or spin systems rabl2010quantum ; camerer2011realization ; arcizet2011single ; hill2012coherent .

Furthermore, we have tested the robustness of our proposal against several imperfections that may occur in a realistic experimental scenario (see Supplementary Sec. S4). The most important limitation would be the presence of a considerable amount of flux noise, resulting in stray qubit-qubit coupling; for example, adding a fluctuation of $\delta\Phi_{\text{b}}=1-10~{}\mu\Phi_{0}$ results in $0.1-1$ MHz added qubit-qubit coupling $J$, respectively. We find that for $J<1$ MHz the fidelity of the cooling and quantum state preparation protocols is not compromised (Supplementary Fig. S4), therefore $\delta\Phi_{\text{b}}<10~{}\mu\Phi_{0}$ is required, which is compatible with observations in similar devices kumar2016origin ; hutchings2017tunable . We note that despite the steep slope of the tripartite coupling $g$ versus flux bias for $\Phi_{\text{b}}/\Phi_{0}>0.495$ (Fig. 1b), the tripartite coupling never changes by more than $1\%$ for $\delta\Phi_{\text{b}}<10~{}\mu\Phi_{0}$. Another possible experimental limitation is the deviation from the target qubit frequencies due to imperfect flux tuning pulses. We have studied the effect of this imperfection and find that targeting the qubit frequencies within 100 kHz is sufficient for high-fidelity quantum state preparation (see Supplementary Fig. S5). Additionally, we have studied the robustness of the protocol against qubit coherence and we find that high-fidelity quantum state preparation can be obtained for relaxation and dephasing times $T_{1,2}\gtrsim 10~{}\mu$s, which are typical in the superconducting qubit community and compatible with 10 mT magnetic fields schneider2019transmon .

In conclusion, we have analysed a hybrid circuit architecture featuring strong and tuneable flux-mediated electromechanical interactions between a mechanical resonator and two superconducting transmon qubits. Using state-of-the-art parameters, we find that the coupled system can operate in the single-photon strong-coupling regime, which has been a long-standing goal in the field of optomechanics. We have proposed and numerically demonstrated several protocols for achieving ground-state cooling and preparing multi-phonon quantum superposition states as well as hybrid entanglement with high fidelities, which has been a tremendous challenge so far. Moreover, the proposed schemes for quantum manipulation are applicable to a wider range of tripartite quantum systems where a lower frequency mode, that is not directly accessible, is within the tuning range of the two other controllable modes. Our work significantly extends the quantum control toolbox of parametrically coupled radio-frequency mechanical resonators and provides a versatile on-chip interface with transmon-based processors, offering rich opportunities for technological applications as well as fundamental tests of quantum mechanics.

IV Methods

IV.1 Numerical modelling

We model the dynamical evolution of the system, including environmental dissipation, with the Lindblad master equation

[TABLE]

which is numerically solved using QuTiP johansson2012qutip . Here, $\mathcal{L}[\hat{o}]\rho\doteq(2\hat{o}\rho\hat{o}^{\dagger}-\hat{o}^{\dagger}\hat{o}\rho-\rho\hat{o}^{\dagger}\hat{o})/2$ are superoperators describing each dissipation process, and $n_{\text{th}}~{}=~{}1/[\exp(\hbar\omega_{m}/(k_{\text{B}}T))-1]$ is the thermal phonon number at temperature $T$. More specifically, we consider qubit decay and dephasing times $T_{1}=T_{2}=30~{}\mu$s, which are consistent with measured values in a similar tuneable coupling transmon architecture kounalakis2018tuneable . The coupling of the mechanical mode to the environment is determined by $\gamma_{m}=\omega_{m}/Q$, where the quality factor $Q=10^{6}$ is chosen in agreement with experimental observations in recently fabricated SQUID-embedded beams rodrigues2019coupling . For completeness, we additionally include $\mathcal{O}(\phi^{4}X)$ terms in the interaction Hamiltonian (see Supplementary Sec S2), which nevertheless cause insubstantial corrections to the system dynamics. We model the mechanical resonator using forty levels and each transmon as a three-level system (including an anharmonicity of $E_{\text{C}_{i}}/h\simeq 320$ MHz). The same parameters, shown in Table 1, were considered in all the simulations.

V Acknowledgements

We thank M.F. Gely, D. Bothner and I.C. Rodrigues for useful discussions. This research was supported by the Dutch Foundation for Scientific Research (NWO) through the Casimir Research School.

VI Author contributions

M.K. conceptualized the work and developed the theory and numerical simulations with supervision from Y.M.B. and G.A.S.; M.K. wrote the manuscript with input from all coauthors.

VII Data and code availability statement

The simulation code and data sets generated and analysed during the current study supporting the main and supplementary figures are publicly available in Zenodo with the identifier 10.5281/zenodo.3469853.

VIII Lagrangian-Hamiltonian description of the circuit

The Lagrangian describing the electromechanical system in Fig.1(a) is

[TABLE]

where $X,~{}\phi_{\text{i}}$ are variables representing the beam displacement and the flux on circuit node $i$, respectively, and $\phi_{0}=\hbar/2e$ is the reduced flux quantum. $C_{i},~{}C_{\text{c}}$ and $E_{\text{J}_{i}},~{}E_{\text{J}}^{\text{c}}$ denote the capacitances and Josephson energies of each transmon and the coupler, respectively, and $m,~{}\omega_{0}$ are the mass and frequency of the beam. Following a Legendre transformation $H=\sum_{i}\phi_{i}Q_{i}-\mathcal{L}$ we obtain the system Hamiltonian

[TABLE]

where $P~{}=~{}\frac{\partial\mathcal{L}}{\partial\dot{X}}$ is the mechanical conjugate momentum and $Q_{i}~{}=~{}\frac{\partial\mathcal{L}}{\partial\dot{\phi}_{i}}$ are the electrical conjugate momenta representing charges on each circuit node. $\widetilde{C}_{2}=C_{2}+\frac{C_{1}C_{\text{c}}}{C_{1}+C_{\text{c}}}$ and $\widetilde{C}_{1}=C_{1}+\frac{C_{2}C_{\text{c}}}{C_{2}+C_{\text{c}}}$ denote the modified transmon capacitances due to the coupling capacitance $C_{\text{c}}$.

VIII.1 Motion-dependent flux-tuneable Josephson energy

If the two junctions of the SQUID are identical, the total Josephson energy is $E_{\text{J}}^{\text{c}}(\Phi_{\text{b}})=E_{\text{J},\Sigma}^{\text{c}}|\cos(\pi\Phi_{\text{b}}/\Phi_{0})|$, where $E_{\text{J},\Sigma}^{\text{c}}$ is the sum of the Josephson energies of each junction. In the presence of an in-plane magnetic field, $B$, the SQUID loop can pick up an additional flux, $\beta_{0}BlX$, due the beam displacement, where $l$ is the length of the suspended arm and $\beta_{0}$ is a geometric constant associated with the mode shape of the beam (for the first mechanical mode $\beta_{0}\sim 1$) etaki2008motion ; rodrigues2019coupling . Therefore, the flux- and motion-dependent Josephson energy reads

[TABLE]

where $\alpha~{}=~{}\pi\beta_{0}Bl/\Phi_{0}$ and $\Phi_{0}~{}=~{}h/2e$ is the magnetic flux quantum. Using basic trigonometry and expanding to lowest order in $X$, assuming $\pi\Phi_{\text{b}}/\Phi_{0}\gg\alpha X$, we find

[TABLE]

The above approximation holds for $\alpha X_{\text{ZPF}}\ll 1$, which is valid for the beams considered in similar experiments and small magnetic fields ($B<1$ T) rodrigues2019coupling .

In the above, the case of a symmetric SQUID $E_{\text{J,1}}^{\text{c}}=E_{\text{J,2}}^{\text{c}}$ is considered, however in realistic devices a finite asymmetry needs to be taken into account. The Josephson energy is therefore more accurately described by

[TABLE]

where $a_{\text{J}}=|({E_{\text{J,1}}^{\text{c}}-E_{\text{J,2}}^{\text{c}}})/({E_{\text{J,1}}^{\text{c}}+E_{\text{J,2}}^{\text{c}}})|$ is the SQUID asymmetry koch2007charge . For $\pi\Phi_{\text{b}}/\Phi_{0}\gg\alpha X$, we have

[TABLE]

Substituting into Eq. (S16), yields

[TABLE]

which can further be simplified, assuming $\alpha X\frac{(1-a_{\text{J}}^{2})\tan(\pi\Phi_{\text{b}}/\Phi_{0})}{1+a_{\text{J}}^{2}\tan^{2}{(\pi\Phi_{\text{b}}/\Phi_{0})}}\ll 1$, into

[TABLE]

where $c_{\text{J}}=\sqrt{1+a_{\text{J}}^{2}\tan^{2}{(\pi\Phi_{\text{b}}/\Phi_{0})}}$ and $s_{\text{J}}=(1-a_{\text{J}}^{2})/c_{\text{J}}$. The higher order terms

[TABLE]

are considered in the analysis below, however, as we will see later, they do not affect the system dynamics for the parameters used in this work.

VIII.2 Flux-mediated interactions

Expanding the last term in Eq. (S13) up to $\mathcal{O}[\phi^{4}X^{2}]$, in combination with Eq. (S19), yields the following flux-mediated interaction terms

[TABLE]

The first term describes a three-body interaction, between the two qubits and the beam

[TABLE]

while the second term

[TABLE]

describes a radiation-pressure type coupling of the mechanical mode with each qubit.

The next two terms describe flux-mediated qubit-qubit interactions, where

[TABLE]

is a linear dipole coupling term, and

[TABLE]

is a nonlinear interaction including cross-Kerr as well as pair- and correlated-hopping terms.

The last three terms

[TABLE]

[TABLE]

and

[TABLE]

describe nonlinear tripartite interactions, which are much weaker and do not contribute significantly to the dynamics compared to the leading order electromechanical terms of Eq. (S22) and (S23).

In addition, the last term in Eq. (S13) leads to corrections in the bare subsystem Hamiltonians. More specifically, the inductive energy term

[TABLE]

results in an effective qubit Josephson energy

[TABLE]

while the potential energy term

[TABLE]

leads to a displaced rest position

[TABLE]

The latter does not affect the dynamics and can be absorbed in a redefinition of $X\rightarrow(X-X_{0})$.

IX Circuit quantisation

We now switch to a quantum mechanical description of the circuit, promoting all canonical variables to quantum operators

[TABLE]

where $\hat{c}_{i}^{(\dagger)}$ and $\hat{b}^{(\dagger)}$ are ladder operators describing the annihilation (creation) of photons and phonons, respectively, and satisfy bosonic commutation relations $[\hat{c}_{i},~{}\hat{c}_{j}^{\dagger}]~{}=~{}\delta_{ij}$ and $[\hat{b},~{}\hat{b}^{\dagger}]~{}=~{}1$. The zero-point fluctuations in the mechanical displacement and momentum are given by $X_{\text{ZPF}}~{}=~{}\sqrt{{\hbar}/{(2m\omega_{\text{m}})}}$ and $P_{\text{ZPF}}~{}=~{}\sqrt{{\hbar m\omega_{\text{m}}}/{2}}$, respectively. $Z_{i}=\frac{\hbar}{e^{2}}\sqrt{E_{\text{C}_{i}}/2\widetilde{E}_{\text{J}_{i}}}$ denotes the impedance of each transmon, where $E_{\text{C}_{i}}~{}=~{}\frac{e^{2}}{2\widetilde{C}_{i}}$ is its charging energy. Since the qubits are in the transmon regime koch2007charge , $\widetilde{E}_{\text{J}_{i}}\gg E_{\text{C}_{i}}$, the uncoupled electromechanical system is well-described by

[TABLE]

with qubit frequencies $\omega_{i}~{}=~{}\frac{1}{\hbar}\left(\sqrt{8\widetilde{E}_{\text{J}_{i}}E_{\text{C}_{i}}}-E_{\text{C}_{i}}\right)$.

IX.1 Tripartite coupling

Following Eq. (S22) we can now write the Hamiltonian operator describing the three-body interaction

[TABLE]

where, in the second step, we have made a rotating-wave approximation (RWA) to neglect fast rotating terms $\hat{c}_{\text{1}}\hat{c}_{\text{2}}\hat{b}^{(\dagger)}$ and $\hat{c}_{\text{1}}^{\dagger}\hat{c}_{\text{2}}^{\dagger}\hat{b}^{(\dagger)}$, which do not contribute to the dynamics since $\omega_{1},\omega_{2}\gg\omega_{\text{m}},g$. The tripartite coupling strength is given by

[TABLE]

IX.2 Radiation-pressure couplings

The radiation-pressure interaction between each qubit and the beam, in Eq. (S23), can be expressed in second quantisation form as

[TABLE]

following a RWA. The radiation-pressure coupling strengths are given by

[TABLE]

IX.3 Qubit-qubit couplings

Following circuit quantisation and a RWA, the linear qubit-qubit interaction Hamiltonian in Eq. (S24) becomes

[TABLE]

where

[TABLE]

The qubits also couple via their charge degrees of freedom ($Q_{1}Q_{2}$ term in Eq. (S13)), which results in the same type of linear interaction with coupling strength

[TABLE]

Combining these together leads to an overall exchange-type interaction

[TABLE]

which can be suppressed at the desired operating point ($\Phi_{\text{b}}$) with the right choice of coupling capacitor $C_{\text{c}}$.

The nonlinear interaction between the two transmons in Eq. (S25) is given by

[TABLE]

The first term, following a RWA, results in a cross-Kerr interaction, $V\hat{c}_{1}^{\dagger}\hat{c}_{1}\hat{c}^{\dagger}_{2}\hat{c}_{2}$, with coupling strength

[TABLE]

Additionally the same term yields a pair-hopping interaction $\frac{V}{4}(\hat{c}_{1}^{\dagger}\hat{c}_{1}^{\dagger}\hat{c}_{2}\hat{c}_{2}+\text{H.c.})$, which however is not contributing to the dynamics for single-excitation levels. Finally, the other two terms in Eq. (S43) result in correlated hopping interactions $[\hat{c}_{1}^{\dagger}(J_{n_{1}}\hat{c}_{1}^{\dagger}\hat{c}_{1}+J_{n_{2}}\hat{c}_{2}^{\dagger}\hat{c}_{2})\hat{c}_{2}+\text{H.c.}]$, with

[TABLE]

which contribute to the linear coupling in Eq. (S42) as $J\rightarrow J+J_{n_{1}}+J_{n_{2}}$.

IX.4 Higher-order tripartite interactions

The next-to-leading order electromechanical interactions are given by

[TABLE]

with nonlinear coupling strengths

[TABLE]

[TABLE]

and

[TABLE]

The last two terms also lead to a small correction of the tripartite coupling strength $g\rightarrow g+3g_{\phi_{1}^{3}\phi_{2}x}+3g_{\phi_{1}\phi_{2}^{3}x}$. These terms, although weaker, are not negligible (since they are of the same order or larger than dissipation rates) and they are include in simulating the system dynamics.

The $\mathcal{O}[\phi^{2}X^{2}]$ terms in Eq. (S27) can be written as,

[TABLE]

where

[TABLE]

For the parameters used in this work ($\alpha X_{\text{ZPF}}\sim 10^{-6}$) the coupling is negligible and does not affect the system dynamics (the same holds for the even weaker $\mathcal{O}[\phi^{4}X^{2}]$ terms). Interestingly, however, for $i=j$ we find a dispersive (cross-Kerr) interaction $\hat{c}_{i}^{\dagger}\hat{c}_{i}\hat{b}^{\dagger}\hat{b}$ of each qubit with the resonator, which is active even when the qubits are far detuned and could potentially be employed for phonon-sensitive measurements of the mechanical state viennot2018phonon .

X Protocols for arbitrary quantum state generation

X.1 States with arbitrary complex coefficients

In the main text we demonstrated a protocol for creating multi-phonon quantum superposition states by alternating the qubit frequencies such that different parts of the tripartite interaction become resonant. By controlling the interaction times and post-selecting on the qubit state at the end of the protocol, we showed the possibility of creating interesting classes of superposition states with high fidelity. However, this protocol alone is not sufficient for generating mechanical superposition states $|\psi\rangle~{}=~{}\sum_{n=0}^{N}c_{n}|n_{\text{m}}\rangle$ with arbitrary complex coefficients $c_{n}$, since there is no phase degree of freedom to control in the protocol. This is due to the fact that the coupling constant $g$ is not an adjustable complex number and in practical implementations the interaction time is the only parameter that can be varied.

However, there is another tuning knob that can be employed by taking advantage of the qubit-qubit exchange-type interaction $J~{}(\hat{c}_{1}^{\dagger}\hat{c}_{2}+\hat{c}_{1}\hat{c}_{2}^{\dagger})$. Detuning the qubits by $\Delta=|\omega_{1}-\omega_{2}|\gg\omega_{\text{m}}$, while tuning the coupling SQUID such that there is a finite exchange-type coupling strength $J\geq\Delta$, results in a resonant qubit-qubit interaction that couples $|0_{1}1_{2}\rangle$ and $|1_{1}0_{2}\rangle$. For example, starting from $|0_{1}1_{2}\rangle$ one would end up with state $[\cos(Jt)|0_{1}1_{2}\rangle-i\sin(Jt)|1_{1}0_{2}\rangle]$ after time $t$ (at $t=\pi/(2J)$ this realises a SWAP gate). Furthermore, by detuning the qubit frequencies such that $\Delta>J$ it is also possible to introduce a relative phase between the qubit states $|0_{1}1_{2}\rangle$ and $|1_{1}0_{2}\rangle$ (C-Phase gate). Combining this gate with the resonant qubit-qubit interaction, acting on $|0_{1}1_{2}\rangle$ would result in

[TABLE]

We note that our system, is analogous to the one studied theoretically by Law and Eberly law1996arbitrary and experimentally by Hofheinz et al. hofheinz2009synthesizing . The system studied in these references considers a resonator that is controllably coupled to a qubit (with local qubit driving) via a resonant exchange-type interaction (Jaynes-Cummings). It is shown that an arbitrary resonator state can be generated by interleaving the Jaynes-Cummings evolution with qubit driving. The two-qubit states $|0_{1}1_{2}\rangle$ and $|1_{1}0_{2}\rangle$ in our system can be mapped to the single qubit basis $|0\rangle$ and $|1\rangle$ considered in Law and Eberly. Furthermore, the tripartite interaction can be mapped to the qubit-resonator Jaynes-Cummings interaction, where in our system we additionally have the possibility of realising the equivalent of a counter-rotating Jaynes-Cummings interaction by exchanging the qubit frequencies. The equivalent of qubit driving can then be performed by controlling the evolution under the exchange qubit-qubit interaction $U_{\text{J}}$ described in Eq. (S52).

We will consider the case where an arbitrary state $|\psi\rangle_{\text{m}}=\sum_{n=0}^{N}a_{n}|n_{\text{m}}\rangle$ is generated following post-selection from the arbitrary entangled state $|\psi\rangle=\sum_{n=0}^{N-1}(a_{n}|0_{1}n_{\text{m}}1_{2}\rangle+b_{n}|1_{1}n_{\text{m}}0_{2}\rangle)+a_{N}|0_{1}N_{\text{m}}1_{2}\rangle$. That is, if we can create the above state $|\psi\rangle$ with arbitrary complex coefficients, it is then straightforward to collapse it to $|\psi\rangle_{\text{m}}$ following post-selection on the $|0_{1}1_{2}\rangle$ qubit state. Inspired by Refs. law1996arbitrary ; hofheinz2009synthesizing , let us now consider the problem of generating the arbitrary state $|\psi\rangle$ from the inverse point of view, i.e. by proving that it is always possible to empty it. Suppose we start from $|\psi\rangle=\sum_{n=0}^{N-1}(a_{n}|0_{1}n_{\text{m}}1_{2}\rangle+b_{n}|1_{1}n_{\text{m}}0_{2}\rangle)+a_{N}|0_{1}N_{\text{m}}1_{2}\rangle$ (Fig. S5a). Then by applying the (inverse of the) C-Phase gate and the tripartite interaction $\hat{U}^{\dagger}(t)=e^{i(\hat{c}_{1}^{\dagger}\hat{b}\hat{c}_{2}+\text{H.c})t}$ for phase $\theta_{\text{A}}$ and time $t_{\text{A}}$, the probability amplitude of state $|0_{1}N_{\text{m}}1_{2}\rangle$ becomes

[TABLE]

By appropriately choosing $\theta_{\text{A}}$ and $t_{\text{A}}$, it is possible to make the above probability amplitude zero, such that all phonons occupying the state $|0_{1}N_{\text{m}}1_{2}\rangle$ are transferred to $|1_{1}(N-1)_{\text{m}}0_{2}\rangle$ (see Fig. S5b).

The step above also results in incomplete transfer of phonons between states $|0_{1}n_{\text{m}}1_{2}\rangle$ and $|1_{1}(n-1)_{\text{m}}0_{2}\rangle$ (for $n<N$), after which we are left with the state $|\psi\rangle^{\prime}=\sum_{n=0}^{N-1}(a^{\prime}_{n}|0_{1}n_{\text{m}}1_{2}\rangle+b^{\prime}_{n}|1_{1}n_{\text{m}}0_{2}\rangle)$. Now by combining the (inverse of the) C-Phase gate and resonant qubit-qubit interaction $U_{\text{J}}^{\dagger}(t_{\text{J}},\theta_{\text{J}})|\psi\rangle^{\prime}$ the probability amplitude of state $|1_{1}(N-1)_{\text{m}}0_{2}\rangle$ becomes

[TABLE]

which can be made zero by appropriately choosing $\theta_{\text{J}}$ and $t_{\text{J}}$ such that all phonons occupying $|1_{1}(N-1)_{\text{m}}0_{2}\rangle$ are transferred to $|0_{1}(N-1)_{\text{m}}1_{2}\rangle$ (see Fig. S5b). By applying the above two steps $N$ times it is possible to completely empty the initial state, leading the system to $|0_{1}0_{\text{m}}1_{2}\rangle$ (Fig. S5c). Therefore, reversing the problem, it is possible to generate any arbitrary mechanical state $|\psi\rangle_{\text{m}}=\sum_{n=0}^{N}a_{n}|n_{\text{m}}\rangle$ by creating the arbitrary entangled state

[TABLE]

and post-selecting on $|0_{1}1_{2}\rangle$.

X.2 States with arbitrary phonon number probability distributions

The protocol presented in the previous section can enable the creation of any mechanical quantum state with arbitrary coefficients. Although this protocol is experimentally feasible, it can become complex as it requires a lot of additional tuning to realise the C-Phase and exchange-type gates between the qubits. Here we describe an alternative protocol that only employs the tripartite interactions, therefore requiring only alternating between the qubit frequencies, and post-selective measurements at each step. Reducing the complexity of tuning pulses comes at the cost of not being able to create states with arbitrary complex coefficients, although it is possible to generate states with arbitrary phonon number probability distributions as we see below.

The protocol relies on initially preparing the two qubits in an arbitrary entangled state (Fig.S6a). Assuming the resonator is in the ground state, the tripartite system is initially described by the following wavefunction

[TABLE]

where $\alpha$, $\beta$ are complex numbers that can be chosen arbitrarily and are related by $|\alpha|^{2}+|\beta|^{2}=1$. As described in the previous section, one can prepare this state by activating the exchange-type qubit-qubit interaction $U_{\text{J}}(t,\theta)$. This can be done by detuning the qubits sufficiently such that their frequency difference is much greater than the mechanical frequency, and at the same time smaller than the direct qubit-qubit coupling, i.e. $J\geq|\omega_{1}-\omega_{2}|\gg\omega_{\text{m}}$, which can be adjusted by changing the flux $\Phi_{\text{b}}$ on the coupling SQUID.

The next step, schematically depicted in Fig. S6b, consists of the following three substeps:

1. Flux pulse into resonance condition $\omega_{2}^{+}=\omega_{1}+\omega_{\text{m}}$ for variable time $t_{\text{1}}$

2. Flux pulse into resonance condition $\omega_{2}^{-}=\omega_{1}-\omega_{\text{m}}$ for variable time $t_{\text{2}}$,

3. Projective measurement on one of the qubits and post-selection on $|0_{1}1_{2}\rangle$.

The initial wavefunction is transformed into

[TABLE]

where $\hat{P}~{}=~{}|0_{1}1_{2}\rangle\langle 0_{1}1_{2}|$ and $\hat{U}_{\text{A,B}}(t_{1,2})~{}=~{}e^{-i\hat{H}_{\text{A,B}}t_{1,2}}$, with interaction Hamiltonians $\hat{H}_{\text{A}}~{}=~{}\hat{c}_{1}^{\dagger}\hat{b}^{\dagger}\hat{c}_{2}+\text{H.c.}$ and $\hat{H}_{\text{B}}~{}=~{}\hat{c}_{1}\hat{b}^{\dagger}\hat{c}_{2}^{\dagger}+\text{H.c.}$, respectively.

Below we derive the resulting wavefunction after each substep:

[TABLE]

where $P$ is a normalisation factor determined by the Born rule.

Evidently, following this first step one has enough degrees of freedom ($\alpha,t_{1},t_{2}$) to create states with arbitrary phonon number probability distribution up to three levels, $|\psi\rangle_{1}~{}=~{}\sum_{n=0}^{2}c_{n}|n_{\text{m}}\rangle|0_{1}1_{2}\rangle$. Following a second step (Fig. S6c), we find

[TABLE]

The addition of two more variables $t_{3}$, $t_{4}$ enables the creation of states with arbitrary phonon number probability distribution up to $|4\rangle$.

Suppose we start from an arbitrary mechanical state $|\psi\rangle~{}=~{}\sum_{n=0}^{N-2}c_{n}|n_{\text{m}}\rangle|0_{1}1_{2}\rangle$, which can be created by applying the above protocol $(N/2-1)$ times. Then, following another step we have

[TABLE]

where the new coefficients are determined by the previous ones according to the following relations:

[TABLE]

Therefore, applying this protocol $N/2$ times can lead to the creation of states with arbitrary phonon number probability distribution up to $|N\rangle$.

XI Validity for non-ideal system parameters

Although in the simulations we have considered realistic parameters taken from recent experiments, we realise that in an experimental scenario system parameters such as the qubit coherence may significantly deviate from the ones considered in Table 1. Although bad qubit coherence may not affect the success of the cooling protocol, it can pose limits on the quantum state preparation protocols. We have therefore examined the dependence of the fidelity of the prepared quantum state in Fig. 3 to the ideal evolution $|\psi_{\text{ideal}}\rangle=U(t)|0_{1}0_{\text{m}}1_{2}\rangle=\cos(gt)|0_{1}0_{\text{m}}1_{2}\rangle-i\sin(gt)|1_{1}1_{\text{m}}0_{2}\rangle$ on the qubit coherence. In Fig. S7 we plot the evolution of the fidelity of the prepared density matrix $\rho$ to the ideal one $\rho_{\text{ideal}}=|\psi_{\text{ideal}}\rangle\langle\psi_{\text{ideal}}|$, for different values of $T_{1}$ (assuming $T_{2}=T_{1}$). The dashed curves correspond to $e^{-t/T_{1}}$. We find that high-fidelity state preparation ($>90~{}\%$) can be achieved with qubit coherence times of $T_{1},~{}T_{2}\sim 10~{}\mu$s which are standard in the superconducting qubit community and achievable in the presence of 10 mT magnetic fields schneider2019transmon .

One effect that is not considered in the main text is that of fluctuations on the flux bias channel $\delta\Phi_{\text{b}}$. These can occur as a result of environmental magnetic field noise or noise of the current source used for biasing the flux line. Given the stability of our current sources and assuming flux noise levels reported in similar devices kumar2016origin ; hutchings2017tunable , we estimate such fluctuations to be around $\delta\Phi_{\text{b}}/\Phi_{0}\sim 10^{-6}$ during the course of the preparation protocols. The most important effect of flux fluctuations would be to introduce a stray qubit-qubit exchange-type coupling as $J_{\text{L}}\propto E_{\text{J},\Sigma}^{\text{c}}\cos(\pi(\Phi_{\text{b}}+\delta\Phi_{\text{b}})/\Phi_{0})$. More specifically, adding a fluctuation of $\delta\Phi_{\text{b}}/\Phi_{0}\sim 10^{-6}$ would result in $J_{\text{L}}\rightarrow J_{\text{L}}+100$ kHz, while $\delta\Phi_{\text{b}}/\Phi_{0}\sim 10^{-5}$ and $\delta\Phi_{\text{b}}/\Phi_{0}\sim 10^{-4}$ translate to additional coupling of 1 and 10 MHz, respectively. The latter would be detrimental for the state preparation protocols as it is of the order of the mechanical frequency. Note that the effect of flux noise on the qubit-qubit coupling is also amplified by our choice of a large coupling Josephson energy amplitude $E_{\text{J},\Sigma}^{\text{c}}/h=200$ GHz. In Fig. S8 we plot the effect of a finite qubit-qubit coupling on the cooling and quantum state preparation protocols. We find that for $J<1$ MHz (therefore $\delta\Phi_{\text{b}}<10~{}\mu\Phi_{0}$) it is possible to maintain high-fidelity quantum state preparation. Note that despite the sharp dependence of the tripartite coupling on the flux bias (Fig. 1), this change is less than $0.1~{}\%$ in the case of $\delta\Phi_{\text{b}}<10~{}\mu\Phi_{0}$.

Additionally, in Fig. S9 we examine the effect of imperfect qubit flux pulsing on the fidelity of the prepared quantum state. We find that the fidelity of the protocol is sensitive to this parameter and that targeting the qubit frequencies within $\sim 100$ kHz is required for high-fidelity quantum state preparation.