Protocol for generating an arbitrary quantum state of the magnetization in cavity magnonics

Magnets have found commercial applications in magnetic field sensing and information storage, and are a promising building block for long-range information transfer [1, 2] and low-power logic devices [3]. Recently, there has been an interest in bringing such applications into the quantum domain, known as 'quantum magnonics' [4–7], in which the unique properties of the collective excitation of the magnetization, called magnons, are harnessed. Such properties include a high degree of frequency tunability, and low dissipation. The latter is a particular characteristic to the material of choice for quantum magnonics, the ferrimagnetic insulator yttrium iron garnet (YIG) [8], in which evidence for macroscopic (mm-long) coherence lengths were observed [9, 10]. Furthermore, magnons can couple coherently to other systems, such as optical light [11–14], microwaves (MWs) [15–17] and mechanical vibrations [18–20]. In particular, the strong coupling regime between cavity MW photons and magnons, where information can be coherently exchanged between the two excitations at rates $\gt$10 MHz much faster than each of their decay rates $\lt$1 MHz, can be routinely achieved in experiments. In contrast, the strong coupling regime in the optical domain is notoriously more difficult to achieve, however there are promising developments in this direction [21–26]. Possible applications of such hybrid magnonic quantum systems include magnon-based quantum transducers between light and MWs, as well as transduction involving other degrees of freedom such as phonons [20, 27] and electrons [28–30]. Moreover, platform specific applications such as enhanced magnetic field sensing by using non-classical states of the magnetization can be envisioned [31], potentially reaching the Heisenberg limit [32].

A quantum platform involving magnets requires generation, manipulation, and detection of non-classical states of the magnetization. While single-magnon detection has been demonstrated [33–35] by exploiting the coupling of magnets to superconducting circuits mediated by MW cavities [4, 36–38], experiments so far have probed coherent or thermal magnon states whereas 'true' non-classical states (defined e.g. by a Wigner function with negative regions) remain to be realized. Theoretical proposals in this direction include the generation of cat states of the magnetization using MW cavity photons [39], the generation of single magnon Fock states [40] and cat states [41] via optical means, and the generation of cat states using qubit [42]. In thin films, the generation of entangled pairs of travelling magnons was proposed [43]. All of these proposals are specialized to their respective target states, while there is no known method to deterministically generate an arbitrary magnon state.

In this manuscript we fill in this gap by proposing a protocol to this end, tailored for magnets coupled to a superconducting transmon via a MW cavity, as depicted in figure 1. Our method is based on the protocol proposed in [44], which shows that the deterministic generation of an arbitrary quantum state of a harmonic oscillator can be accomplished by coupling the oscillator to a qubit. A superconducting transmon is a weakly anharmonic oscillator whose first two levels can be treated as a qubit. We show how the protocol of [44] can be modified to take into account the anharmonicity of the superconducting transmon, in order to generate arbitrary quantum magnon states with high fidelity.

Zoom In Zoom Out Reset image size

The paper is organized as follows. We discuss the system and effective Hamiltonian in section 1. In section 2, we approximate the superconducting transmon by a two-level system and review the protocol from [44] for the cavity mangnonic system. In section 3, we discuss how to correct for errors due to deviations from the two-level approximation, in particular due to interferences from higher energy levels present in the transmon. In section 4, we discuss the fidelity of state generation for paradigmatic examples of non-classical states. We discuss possible improvements to our protocol that can increase the fidelity of generation further in section 5. We review the main findings qualitatively in section 6.

The hybrid system we considered is depicted in figure 1: inside a MW cavity, magnons can couple to a superconducting transmon [4, 36]. The effective coupling Hamiltonian can be obtained by tracing out the MW field. We assume a strong coupling regime, such that the dissipation can be ignored while calculating the effective Hamiltonian. In this section, we present the final results while the detailed calculations are provided in appendix A.

When the magnet is much smaller than the wavelength of the MWs (∼cm), only the total magnetization M couples to the MWs, so we need to consider only the uniform magnon mode known as the Kittel mode. It can be quantized via a Holstein–Primakoff transformation [45, 46] that, in the limit of $M_{x,y} \ll M_z$, reduces to:

$$\begin{align} M_x - iM_y\rightarrow2\mathcal{M}_\textrm{ZPF}\skew{2.4}\hat{M},\ M_z\rightarrow M_s - \frac{2\mathcal{M}_\textrm{ZPF}^2}{M_s}\skew{2.4}\hat{M}^{\dagger}\skew{2.4}\hat{M}, \end{align} \tag{ 1 }$$

where $\skew{2.4}\hat{M}$ is the magnon annihilation operator satisfying the bosonic commutation relations $[\skew{2.4}\hat{M},\skew{2.4}\hat{M}^{\dagger}] = 1$, M_s is the saturation magnetization, and the zero-point fluctuations of the magnetization are given by:

$$\begin{align} \mathcal{M}_\textrm{ZPF} = \sqrt{\frac{\gamma\hbar M_s}{2V_m}}, \end{align} \tag{ 2 }$$

with γ > 0 being the absolute gyromagnetic ratio and V_m being the volume of the magnet. Assuming a spherical magnet and an external magnetic field H_app, the Hamiltonian is $\skew{2.4}\hat{H}_\textrm{m} = \hbar\Omega_m\skew{2.4}\hat{M}^{\dagger}\skew{2.4}\hat{M}$ with $\Omega_m = \gamma\mu_0H_{\textrm{app}}$ which can be tuned to be in the range of ∼2π × 5–10 GHz. For a YIG sample with volume $V_m \gt 1~\mu$m³, $M_s = 140$ kA m⁻¹, and $\gamma = 2\pi\times 28$ GHz T⁻¹, we have $\mathcal{M}_\textrm{ZPF}/M_s\lt8\times10^{-6}$, implying that the deviations from the ground state are small.

A flux-tunable transmon consists of two Josephson junctions in parallel with Josephson energy $E_{J1,J2}$ respectively, forming a SQUID loop along with a capacitor with capacitance C also in parallel [47–50]. When the charging energy $E_C = e^2/2C$ is such that $E_C\ll E_{J1,J2}$ , the energy levels of the circuit resemble that of a weakly anharmonic oscillator,

$$\begin{align} E_n \approx \hbar\left(\Omega_q - \frac{A}{2}(n - 1)\right)n, \end{align} \tag{ 3 }$$

where $A = E_C/4$ and $\Omega_q = \sqrt{E_CE_{J,\textrm{eff}}}$ with the effective Josephson energy $E_{J,\textrm{eff}}$ being:

$$\begin{align} E_{J,\textrm{eff}}(\Phi) = \sqrt{\left(E_{J1} + E_{J2}\right)^2\cos^2\frac{\pi\Phi}{\Phi_0} + \left(E_{J1} - E_{J2}\right)^2\sin^2\frac{\pi\Phi}{\Phi_0}}. \end{align} \tag{ 4 }$$

Here Φ is the flux inside the SQUID loop and $\Phi_0 = h/2e$ is the superconducting flux quantum. Such a set of energy levels can be labelled using an annihilation operator $\skew{2.4}\hat{Q}$ and its corresponding number operator $\skew{2.4}\hat{Q}^{\dagger}\skew{2.4}\hat{Q}$. The number of quanta exchanged between the two Josephson junctions is given by:

$$\begin{align} \hat{N}_{\textrm{sc}} = \left(\frac{E_{J,\textrm{eff}}}{2E_C}\right)^{1/4}\frac{\skew{2.4}\hat{Q} - \skew{2.4}\hat{Q}^{\dagger}}{\sqrt{2}i}. \end{align} \tag{ 5 }$$

The Hamiltonian becomes:

$$\begin{align} \skew{2.4}\hat{H}_\textrm{tr} = \hbar\left[ \Omega_q - \frac{A}{2}\left(\skew{2.4}\hat{Q}^{\dagger}\skew{2.4}\hat{Q} - 1\right) \right] \skew{2.4}\hat{Q}^{\dagger}\skew{2.4}\hat{Q}. \end{align} \tag{ 6 }$$

The frequency $\Omega_q$ can be tuned dynamically via the external flux in SQUID Φ, to be in a typical range of $2\pi\times5$–20 GHz. Tuning via Φ can be achieved at a fast timescale $\ll$10 ns [51, 52], but it comes at the cost of increased decoherence proportional to the allowed tunability $\left|E_{J1} - E_{J2}\right|$, as a higher tunability implies larger susceptibility to flux noises [49, 53]. For our case, we assume that $\Omega_q$ is tunable within a range of $2\pi\times1$ GHz.

We assume that the magnet and the transmon are both placed inside a MW cavity which can be modelled as a set of harmonic modes with frequencies $\Omega_{a,r}$ and annihilation operator $\skew{2.4}\hat{A}_r$. The details of the modes depend on the shape and size of the cavity. The Hamiltonian for the unloaded cavity is $\skew{2.4}\hat{H}_\textrm{mw} = \sum_r\hbar\Omega_{a,r}\skew{2.4}\hat{A}_r^{\dagger}\skew{2.4}\hat{A}_r$. Its coupling to both the magnet and the transmon corresponds to exchange of quanta and can be modelled as:

$$\begin{align} \frac{\skew{2.4}\hat{H}_\textrm{int}}{\hbar} = \sum_r\left(G_{m,r}\skew{2.4}\hat{M}^{\dagger}\skew{2.4}\hat{A}_r + G_{m,r}^*\skew{2.4}\hat{M}\skew{2.4}\hat{A}_r^{\dagger} + G_{q,r}\skew{2.4}\hat{Q}^{\dagger}\skew{2.4}\hat{A}_r + G_{q,r}^*\skew{2.4}\hat{Q}\skew{2.4}\hat{A}_r^{\dagger}\right), \end{align} \tag{ 7 }$$

where the coupling $G_{m,r}\propto\left|\boldsymbol{B}(\boldsymbol{r}_{\textrm{magnet}})\right|$ with $\boldsymbol{B}(\boldsymbol{r}_{\textrm{magnet}})$ the magnetic field at the location of the magnet, and $G_{q,r}\propto\left|\boldsymbol{E}(\boldsymbol{r}_{\textrm{transmon}})\right|$ where $\boldsymbol{E}(\boldsymbol{r}_{\textrm{transmon}})$ is the electric field at the location of the transmon. In deriving the Hamiltonian, we have used a rotating wave approximation to drop the terms $\propto \skew{2.4}\hat{M} \skew{2.4}\hat{A}_r$, and $\propto \skew{2.4}\hat{Q} \skew{2.4}\hat{A}_r$. This is justified since $G_{m(q),r} \ll \Omega_{m(q)} + \Omega_a$ for the frequencies we are considering. If for all r modes, $\Omega_{a,r} - \Omega_m\gt G_{m,r}$ and $\Omega_{a,r} - \Omega_q\gt G_{q,r}$, we can define cavity-dressed magnon and transmon modes with annihilation operators $\skew{2.4}\hat{m}$ and $\skew{2.4}\hat{q}$ respectively. Ignoring contributions that are quartic or higher order in $G_{q,r}/(\omega_a-\omega_q)$ and $G_{m,r}/(\omega_a-\omega_m)$,

$$\begin{align} \frac{\skew{2.4}\hat{H}}{\hbar} = \omega_m\skew{2.4}\hat{m}^{\dagger}\skew{2.4}\hat{m} + \left(\omega_q - \hat{\alpha}\left[\skew{2.4}\hat{q}^{\dagger}\skew{2.4}\hat{q}\right]\right)\skew{2.4}\hat{q}^{\dagger}\skew{2.4}\hat{q} + \skew{2.4}\hat{m}\skew{2.4}\hat{q}^{\dagger}\skew{2.4}\hat{g}\left[\skew{2.4}\hat{q}^{\dagger}\skew{2.4}\hat{q}\right] + \skew{2.4}\hat{g}\left[\skew{2.4}\hat{q}^{\dagger}\skew{2.4}\hat{q}\right]\skew{2.4}\hat{m}^{\dagger}\skew{2.4}\hat{q}. \end{align} \tag{ 8 }$$

Details about this Hamiltonian are provided in appendix A, while here we list down the features relevant for the rest of the article. The renormalized magnon's frequency ω_m and transmon's frequency ω_q stay close to their bare values $\Omega_m$ and $\Omega_q$ respectively. A schematic representation of each term of the Hamiltonian is shown in figure 1.

The function $\hat{\alpha}$ gives the anharmonicity of the transmon and $\skew{2.4}\hat{g}$ gives the magnon-transmon coupling. Below, we use the same characters $\alpha(n) = \langle n|\hat{\alpha}|n\rangle$ and $g(n) = \langle n|\skew{2.4}\hat{g}|n\rangle$ for the values of anharmonicity and coupling. Typically $\alpha(2)\sim2\pi\times150$–300 MHz increasing with n and the coupling g can achieve ∼2π × 20–30 MHz depending on the system parameters.

Restricting the analysis to only the first two levels of the transmon, we get the Jaynes–Cummings Hamiltonian:

$$\begin{align} \frac{\skew{2.4}\hat{H}^{(2)}}{\hbar} = \omega_m\skew{2.4}\hat{m}^{\dagger}\skew{2.4}\hat{m} + \omega_q\left|e\rangle\langle e\right| + g\left(\skew{2.4}\hat{m}\hat{\sigma}_+ + \skew{2.4}\hat{m}^{\dagger}\hat{\sigma}_-\right), \end{align} \tag{ 9 }$$

where $g\equiv g(0)$, $\{|g\rangle,|e\rangle\}$ are the first two excited states of the transmon, $\hat{\sigma}_+ = |g\rangle\langle e|$ and $\hat{\sigma}_- = |e\rangle\langle g|$. As mentioned before, the frequency ω_q is tunable. When $\left|\omega_q - \omega_m\right|\gg g$, the magnons and transmon are decoupled besides a negligible dispersive coupling, while at resonance, $\omega_q = \omega_m$, the two subsystems can exchange quanta.

In this section, we review the protocol to generate an arbitrary state of the magnetization via its coupling to a two-level system. This is discussed in [44] theoretically for a general harmonic oscillator. The protocol consists of: (a) exciting the qubit; (b) partially transferring the excitation to the oscillator and (c) repeating the process multiple times to construct the Fock-state superposition that describes the target state. In this setup, the transmon can be excited via a time-dependent drive, while the transmon and the magnons can be put in and out of resonance by tuning the transmon frequency (via an applied flux), effectively changing the coupling rate (see figure 2). The procedure was experimentally implemented to generate quantum states of MW photons coupled to a superconducting transmon [52].

Zoom In Zoom Out Reset image size

If we start from the qubit's state $c_0|e\rangle + c_1|g\rangle$ while the magnon is in the ground state $|0\rangle$, switching the coupling on for a period $t = \pi/(2g)$ and in the absence of dissipation evolves the state as:

$$\begin{align} \left(c_0|e\rangle + c_1|g\rangle\right)|0\rangle\rightarrow|g\rangle\left(c_0|0\rangle + c_1|1\rangle\right), \end{align} \tag{ 10 }$$

where $|n\rangle$ is n-magnon state. Thus, we can create an arbitrary superposition of $|0\rangle$ and $|1\rangle$. If we repeat this process N times, we can, in principle, generate an arbitrary magnon state containing up to a maximum of N magnons. The protocol then consists of a series of operations in which first the transmon is excited, and then the qubit's frequency is changed to bring it in resonance with the magnons for a time of the order $\pi/g$.

To make the process concrete, consider the Hamiltonian restricted to the first two levels of the transmon $\skew{2.4}\hat{H}^{(2)}(t) = \skew{2.4}\hat{H}_\textrm{s}^{(2)} + \skew{2.4}\hat{H}_\textrm{d}^{(2)}(t)$, where the static part is:

$$\begin{align} \frac{\skew{2.4}\hat{H}_\textrm{s}^{(2)}}{\hbar} = \omega_m\skew{2.4}\hat{m}^{\dagger}\skew{2.4}\hat{m} + \omega_q|e\rangle\langle e| + g\left(\skew{2.4}\hat{m}\hat{\sigma}_+ + \skew{2.4}\hat{m}^{\dagger}\hat{\sigma}_-\right), \end{align} \tag{ 11 }$$

and the dynamic part is:

$$\begin{align} \frac{\skew{2.4}\hat{H}_\textrm{d}^{(2)}}{\hbar} = \Delta(t)|e\rangle\langle e| + \skew{2.4}\tilde{\varepsilon}(t)\hat{\sigma}_+ + \skew{2.4}\tilde{\varepsilon}^*(t)\hat{\sigma}_-, \end{align} \tag{ 12 }$$

here and in what follows, the superscript (2) indicates the approximation of the transmon as a two-level system. In this expression, $\Delta(t)$ is the externally induced change in the transmon's frequency, and $\skew{2.4}\tilde{\varepsilon}(t)$ is the transmon drive amplitude. We assume that the qubit and the magnons are far detuned, $\left|\omega_m - \omega_q\right|\gg g$, implying that magnons and qubit are decoupled at $\Delta = 0$. We review the main steps of the derivation here leaving the detailed calculations to appendices B and C.

A non-zero $\skew{2.4}\tilde{\varepsilon}$ can excite a qubit. The time evolution operator $\hat{\mathcal{Q}}^{(2)}$ for $\Delta = 0$ and $\skew{2.4}\tilde{\varepsilon}(t) = \varepsilon(t)e^{- i\omega_qt + i\phi}$ reads:

$$\begin{align} \hat{\mathcal{Q}}^{(2)}\left(\theta,\phi\right) = e^{- i\omega_mt_g\skew{2.4}\hat{m}^{\dagger}\skew{2.4}\hat{m} - i\omega_qt_g|e\rangle\langle e|} \left({\hat{I}\cos\theta + i\sin\theta\left(e^{i\phi}\hat{\sigma}_+ + e^{- i\phi}\hat{\sigma}_-\right)}\right), \end{align} \tag{ 13 }$$

where $\theta = \int_0^{t_g} \varepsilon(\tau) d\tau$ with t_g a pre-determined gate time. This applies a Bloch rotation of magnitude 2θ around the axis $(\cos\phi, - \sin\phi,0)$.

After the qubit is excited, we want to (partially) transfer the quanta to magnons. If the coupling is switched on for a time $t_\textrm{on}$ ($\Delta(t) = \Delta_0$ where the detuning is $\Delta_0 = \omega_m - \omega_q$) and then switched off for a time $t_\textrm{off}$, ($\Delta(t) = 0$), keeping the qubit driving off $\skew{2.4}\tilde{\varepsilon} = 0$, then the time evolution operator can be written as $\hat{\mathcal{T}}^{(2)} = \hat{\mathcal{T}}_{\textrm{ph}}^{(2)}\hat{\mathcal{T}}_{\textrm{ex}}^{(2)}$ where the phases evolve as:

$$\begin{align} \hat{\mathcal{T}}_{\textrm{ph}}^{(2)} = \exp\left( - i\omega_m\left(t_{\textrm{off}} + t_{\textrm{on}}\right)\skew{2.4}\hat{m}^{\dagger}\skew{2.4}\hat{m} - i\left(\omega_mt_{\textrm{on}} + \omega_qt_{\textrm{off}}\right)|e\rangle\langle e|\right), \end{align} \tag{ 14 }$$

and the quanta are exchanged via:

$$\begin{align} \hat{\mathcal{T}}_{\textrm{ex}}^{(2)} = \cos\left(gt_{\textrm{on}}\sqrt{\hat{N}}\right) - \frac{i}{\sqrt{\hat{N}}}\sin\left(gt_{\textrm{on}}\sqrt{\hat{N}}\right)\left(\skew{2.4}\hat{m}\hat{\sigma}_+ + \skew{2.4}\hat{m}^{\dagger}\hat{\sigma}_-\right), \end{align} \tag{ 15 }$$

with the total number operator $\hat{N} = \skew{2.4}\hat{m}^{\dagger}\skew{2.4}\hat{m} + |e\rangle\langle e|$.

A given target state $|\psi_T\rangle$ containing a maximum of N_m magnons, can be systematically achieved by adding magnons one by one as described above. For that, we need a set of time evolution operators $\hat{\mathcal{Q}}_j^{(2)}\equiv\hat{\mathcal{Q}}^{(2)}(\theta_j,\phi_j)$ and $\hat{\mathcal{T}}_j^{(2)}\equiv\hat{\mathcal{T}}^{(2)}(t_{\textrm{on},j},t_{\textrm{off},j})$ (see figure 2) such that:

$$\begin{align} |\Psi_0\rangle = |0,g\rangle,\ |\Psi_{j + 1}\rangle = \hat{\mathcal{T}}_j^{(2)}\hat{\mathcal{Q}}_j^{(2)}|\Psi_j\rangle,\ |\Psi_{N_m}\rangle = |\psi_T,g\rangle, \end{align} \tag{ 16 }$$

and each $|\Psi_j\rangle$ has a maximum of j excitations. Explicitly, for $|\Psi_j\rangle = \left|\psi_j^e,e\right\rangle + \left|\psi_j^g,g\right\rangle$ there are a maximum of $(j - 1)$ and j magnons in $|\psi_j^e\rangle$ and $|\psi_j^g\rangle$ respectively. The coefficients $\{\theta_j,\phi_j,t_{\textrm{off},j},t_{\textrm{on},j}\}$ are found by reversing the problem in order to systematically remove magnons one by one from $|\Psi_{N_m}\rangle$. Inductively, we assume that $|\Psi_{j + 1}\rangle$ is found. Explicitly, we find $\{t_{\textrm{on},j},t_{\textrm{off},j}\}$ to ensure $\left\langle\, j + 1,g\left|\hat{{\cal T}}_j^{(2),\dagger}\right|\Psi_{j + 1}\right\rangle = 0$ achieved via:

$$\begin{align} \left\langle\, j + 1,g\left|\Psi_{j + 1}\right.\right\rangle\cos\left(gt_{\textrm{on},j}\sqrt{j + 1}\right) + i\left\langle\, j,e\left|\Psi_{j + 1}\right.\right\rangle\sin\left(gt_{\textrm{on},j}\sqrt{j + 1}\right)e^{- i\Delta_0t_{\textrm{off},j}} = 0. \end{align} \tag{ 17 }$$

Next, we find $\{\theta_j,\phi_j\}$ to ensure $\left\langle\, j,e\left|\hat{{\cal Q}}_j^{(2),\dagger}\hat{{\cal T}}_j^{(2),\dagger}\right|\Psi_{j + 1}\right\rangle = 0$ giving:

$$\begin{align} \left\langle\, j,e\left|\hat{{\cal T}}_j^{(2),\dagger}\right|\Psi_{j + 1}\right\rangle\cos\theta_j + i\left\langle\, j,g\left|\hat{{\cal T}}_j^{(2),\dagger}\right|\Psi_{j + 1}\right\rangle\sin\theta_je^{i\phi_j - i\omega_qt_g} = 0. \end{align} \tag{ 18 }$$

Then, we find $|\Psi_j\rangle = \left(\hat{\mathcal{T}}_j^{(2)}\hat{\mathcal{Q}}_j^{(2)}\right)^{\dagger}|\Psi_{j + 1}\rangle$ and the induction continues until all quanta are removed. The procedure to inductively obtain the protocol's parameters is schematically shown in figure 3.

Zoom In Zoom Out Reset image size

The protocol discussed in the previous section is suitable to generate states with only a small number of excitations. States with higher number occupation require longer protocols and hence are spoiled by dissipation. The protocol duration can be decreased via two ways. Firstly, the gate time to excite the qubit can be decreased. Nevertheless, for a transmon qubit, a shorter gate time, $t_g\sim1/\alpha(n)$ where α is the level-dependent anharmonicity (see below equation (8)), implies a broader frequency spectrum allowing for significant non-resonant transitions to higher transmon levels. Secondly, a stronger coupling between magnons and the transmon speeds up the excitation transfer step. However, a very strong coupling, $g(n)\sqrt{n}\sim\alpha(n)$ where g is the level-dependent magnon-transmon coupling (see below equation (8)), implies a leakage of magnon excitations into higher transmon states. Such interferences limit the protocol time and eventually the largest size of the harmonic oscillator state that can be created. In this section, we show that transitions to higher levels can be incorporated by adding extra pulses for cancelling out such transition amplitudes. The extra pulses are analogous to the ones used for correcting qubit gates [54, 55] and can significantly decrease the protocol time allowing for states containing a large magnon number.

For clarity in presentation, we consider the case of $g(n)\equiv g$ and:

$$\begin{align} \alpha(n)\equiv\frac{A}{2}\left(n - 1\right), \end{align} \tag{ 19 }$$

while the calculations for the general case are presented in appendices B and C. The Hamiltonian, equation (8), becomes $\skew{2.4}\hat{H} = \skew{2.4}\hat{H}_\textrm{s} + \skew{2.4}\hat{H}_\textrm{d}(t)$, where the static part is:

$$\begin{align} \frac{\skew{2.4}\hat{H}_\textrm{s}}{\hbar} = \omega_m\skew{2.4}\hat{m}^{\dagger}\skew{2.4}\hat{m} + \left(\omega_q - \frac{A}{2}\left(\skew{2.4}\hat{q}^{\dagger}\skew{2.4}\hat{q} - 1\right)\right)\skew{2.4}\hat{q}^{\dagger}\skew{2.4}\hat{q} + g\left(\skew{2.4}\hat{m}\skew{2.4}\hat{q}^{\dagger} + \skew{2.4}\hat{m}^{\dagger}\skew{2.4}\hat{q}\right), \end{align} \tag{ 20 }$$

and the dynamic part is:

$$\begin{align} \frac{\skew{2.4}\hat{H}_\textrm{d}(t)}{\hbar} = \Delta(t)\skew{2.4}\hat{q}^{\dagger}\skew{2.4}\hat{q} + \skew{2.4}\tilde{\varepsilon}(t)\skew{2.4}\hat{q}^{\dagger} + \skew{2.4}\tilde{\varepsilon}^*(t)\skew{2.4}\hat{q}.\end{align} \tag{ 21 }$$

For a target magnon state $|\psi\rangle$, we want to find a set of operations $\hat{{\cal Q}}_j$ and $\hat{\mathcal{T}}_j$ satisfying (see figure 2):

$$\begin{align} |\Psi_0\rangle = |0,0\rangle,\ |\skew{2.4}\tilde{\Psi}_j\rangle = \hat{\mathcal{Q}}_j|\Psi_j\rangle,\ |\Psi_{j + 1}\rangle = \hat{\mathcal{T}}_j|\skew{2.4}\tilde{\Psi}_j\rangle,\ |\Psi_{N_m}\rangle = |\psi,0\rangle, \end{align} \tag{ 22 }$$

where $|\psi_1,\psi_2\rangle$ refers to the magnon state $|\psi_1\rangle$ and transmon state $|\psi_2\rangle$. As before, starting from $|\Psi_{N_m}\rangle$, we inductively find $\{\hat{\mathcal{T}}_j,\hat{\mathcal{Q}}_j,|\Psi_j\rangle\}$ where the transfer operator removes a magnon, i.e. ensures:

$$\begin{align} \left\langle\, j + 1,0\left|\hat{\mathcal{T}}_j^{\dagger}\right|\Psi_{j + 1}\right\rangle = 0, \end{align} \tag{ 23 }$$

and the gate operator removes transmon excitations, i.e. ensures:

$$\begin{align} \left\langle\, j + 1 - s,s\left|\hat{\mathcal{Q}}_j^{\dagger}\right|\skew{2.4}\tilde{\Psi}_j\right\rangle = 0, \end{align} \tag{ 24 }$$

for $s\geqslant1$.

3.1. Transfer of excitation

We discuss how to find $\hat{\cal T}_j$ satisfying equation (23). Here, we write down the salient features of the derivation in a self-contained manner while the calculations are provided in detail in appendix B.

We start with the same ansatz as in the simpler two-level case, i.e. $\Delta(t) = \Delta_0$ for $t\lt t_\textrm{on}$ and = 0 for $t_\textrm{on} \lt t \lt t_\textrm{on} + t_\textrm{off}$, along with no transmon driving $\skew{2.4}\tilde{\varepsilon} = 0$, for which:

$$\begin{align} \hat{{\cal T}}^{\textrm{ans}}\left(t_{\textrm{on}},t_{\textrm{off}}\right) = \exp\left[\frac{- i\skew{2.4}\hat{H}_\textrm{s}t_{\textrm{off}}}{\hbar}\right]\exp\left[ - \frac{i}{\hbar}\left(\skew{2.4}\hat{H}_\textrm{s} + \Delta_0\skew{2.4}\hat{q}^{\dagger}\skew{2.4}\hat{q}\right)t_{\textrm{on}}\right]. \end{align} \tag{ 25 }$$

This ansatz can be simplified by noticing that the number of excitations $\hat{N} = \skew{2.4}\hat{m}^{\dagger}\skew{2.4}\hat{m} + \skew{2.4}\hat{q}^{\dagger}\skew{2.4}\hat{q}$ are conserved during the evolution, i.e. $\left[\hat{N},\hat{\cal T}^{\textrm{ans}}\right] = \hat{0}$. An expression for $\hat{{\cal T}}^{\textrm{ans}}$ restricted to the finite-dimensional space of $\{|j + 1 - s,s\rangle\}_{s\geqslant0}$ can be found by numerically diagonalizing the two Hamiltonians, $\skew{2.4}\hat{H}_\textrm{s}$ and $\skew{2.4}\hat{H}_\textrm{s}+\Delta_0\skew{2.4}\hat{q}^{\dagger}\skew{2.4}\hat{q}$, again restricted to this space. Finally, we get a condition on $\{t_{\textrm{on}},t_{\textrm{off}}\}$ of the form:

$$\begin{align} \sum_{pqr}e^{i\lambda_{\textrm{on}}^{(p)}t_{\textrm{on}}}e^{i\lambda_{\textrm{off}}^{(q)}t_{\textrm{off}}}a_{pq}^{(r)}\left\langle\, j + 1 - r,r\left|\Psi_{j + 1}\right.\right\rangle = 0, \end{align} \tag{ 26 }$$

where $\{\lambda_{\textrm{on}}^{(p)},\lambda_{\textrm{off}}^{(q)},a_{pq}^{(r)}\}$ can be found in terms of the eigenvectors and eigenvalues of the two Hamiltonians (see equation (B.8)). This equation can be numerically solved for $\{t_{\textrm{on}},t_{\textrm{off}}\}$.

The above ansatz yields a solution in most cases, nevertheless, it can fail for some states, for example, when $|\Psi_{j + 1}\rangle = |j + 1,0\rangle$. Intuitively, when the coupling is on, the excitation from $|j + 1,0\rangle$ transfers to all $|j + 1 - s,s\rangle$ but with different oscillation frequencies $\sim\sqrt{g^2 + \frac{A^2s^2(s - 1)^2}{4}}$. Unless these oscillation frequencies are commensurate, the amplitude of the wave-function in $|j + 1,0\rangle$ will stay non-zero. In cases where this ansatz fails, we can apply a second operation,

$$\begin{align} \hat{{\cal T}}_j = \hat{{\cal T}}^{\textrm{ans}}\left(t_{\textrm{on},j}^{(2)},t_{\textrm{off},j}^{(2)}\right)\hat{{\cal T}}^{\textrm{ans}}\left(t_{\textrm{on},j}^{(1)},t_{\textrm{off},j}^{(1)}\right), \end{align} \tag{ 27 }$$

where we choose $\{t_{\textrm{on},j}^{(1)},t_{\textrm{off},j}^{(1)}\}$ as the values found from a two-level approximation (see equation (18)), and $\{t_{\textrm{on},j}^{(2)},t_{\textrm{off},j}^{(2)}\}$ are found by applying the above optimization to $\hat{{\cal T}}^{\textrm{ans},\dagger}\left(t_{\textrm{on},j}^{(1)},t_{\textrm{off},j}^{(1)}\right)|\Psi_{j + 1}\rangle$. In principle, this process can be continued until a solution is found. However, in our simulations, we never had to go beyond a second step.

3.2. Transmon gates

We discuss the main steps of deriving $\hat{\cal Q}_j$ that satisfies equation (24), while the detailed calculations are provided in appendix C.

As noted before, the application of $\hat{\mathcal{Q}}_j$ ensures that $\hat{\mathcal{Q}}_j^{\dagger}|\skew{2.4}\tilde{\Psi}_j\rangle$ has no excitations above the level j. We systematically apply gates transferring amplitudes from higher excitation states to lower ones, i.e. (suppressing the index j):

$$\begin{align} \hat{\mathcal{Q}}_j^{\dagger} = \dots \hat{\mathcal{Q}}_{12}^{\dagger} \hat{\mathcal{Q}}_{01}^{\dagger}, \end{align} \tag{ 28 }$$

where we want to ensure (see condition equation (24)) $\left\langle\, j,1\left|\hat{\mathcal{Q}}_{01}^{\dagger}\right|\skew{2.4}\tilde{\Psi}_j\right\rangle = 0$, $\left\langle\, j-1,2\left|\hat{\mathcal{Q}}_{12}^{\dagger}\hat{\mathcal{Q}}_{01}^{\dagger}\right|\skew{2.4}\tilde{\Psi}_j\right\rangle = 0$ and so on. First, we discuss how to implement $\hat{\mathcal{Q}}_{01}$.

We try the ansatz $\Delta = 0$ and $\skew{2.4}\tilde{\varepsilon}(t) = -\varepsilon(t)e^{-i\omega_q t + i\phi_\textrm{app}}$ where $\varepsilon(t)$ and φ_app are real, giving:

$$\begin{align} \hat{\mathcal{Q}}_{01}(t) = \mathrm{Texp}\left[ - \frac{i}{\hbar}\int_0^tdt\left(\skew{2.4}\hat{H}_\textrm{s} - \hbar\varepsilon(\tau)e^{- i\omega_q\tau + i\phi_{\textrm{app}}}\skew{2.4}\hat{q}^{\dagger} +\textrm{h.c.}\right) \right], \end{align} \tag{ 29 }$$

We can separate the phases by:

$$\begin{align} \hat{{\cal Q}}_{01}(t) \approx \exp\left[ - i\left(\omega_m\skew{2.4}\hat{m}^{\dagger}\skew{2.4}\hat{m} + \omega_q\skew{2.4}\hat{q}^{\dagger}\skew{2.4}\hat{q}\right)t\right]e^{i\phi_{\textrm{app}}\skew{2.4}\hat{q}^{\dagger}\skew{2.4}\hat{q}}\hat{{\cal U}}_{01}(t)e^{- i\phi_{\textrm{app}}\skew{2.4}\hat{q}^{\dagger}\skew{2.4}\hat{q}}, \end{align} \tag{ 30 }$$

where we ignore the interaction between magnons and the transmon, since $g\ll |\omega_m - \omega_q|$, and:

$$\begin{align} \hat{{\cal U}}_{01}(t) = \mathrm{Texp}\left[ i\int_0^td\tau\ \varepsilon(\tau)\left(\sum_{n = 0}^{\infty}\sqrt{n + 1}\ e^{- inA\tau}|n + 1\rangle\langle n| + \textrm{h.c.}\right) \right]. \end{align} \tag{ 31 }$$

For $At_g\gg 1$, where t_g is the gate time, the solution is the same as that in the two-level approximation, choosing $\varepsilon^{(0)}$ and $\phi_\textrm{app}^{(0)}$ satisfying:

$$\begin{align} \left\langle\, j,1\left|\skew{2.4}\tilde{\Psi}_j\right.\right\rangle \cos\theta + i \left\langle\, j,0\left|\skew{2.4}\tilde{\Psi}_j\right.\right\rangle \sin\theta e^{i\left(\phi_\textrm{app}^{(0)} - \omega_q t_g\right)} = 0, \end{align} \tag{ 32 }$$

where $\theta = \int_0^{t_g} \varepsilon^{(0)}(\tau)d\tau$. For concreteness, we choose:

$$\begin{align} \varepsilon^{(0)}(t) = \frac{\pi\theta}{2t_g}\sin\left(\frac{\pi t}{t_g}\right). \end{align} \tag{ 33 }$$

The approximate time-dependence is given by:

$$\begin{align} \hat{{\cal U}}_{01}^{(0)}(t) = \cos\Theta(t)\left(|0\rangle\langle 0| + |1\rangle\langle 1|\right) + i\sin\Theta(t)\left(|1\rangle\langle 0| + |0\rangle\langle 1|\right) + \sum_{n = 2}^{\infty}|n\rangle\langle n|, \end{align} \tag{ 34 }$$

where $\Theta(t) = \int_0^t\varepsilon^{(0)}(\tau)d\tau$.

However, a drive amplitude with a general form $\varepsilon(t)$ can induce transitions to higher transmon levels, e.g. transitions between $|j,1\rangle$ and $|j,2\rangle$. This is problematic as a significant transition amplitude $\left\langle\, j,2\left|\hat{{\cal Q}}_j^{\dagger}\right|\skew{2.4}\tilde{\Psi}_j\right\rangle$ implies excitation number j + 2, while we want $\hat{\mathcal{Q}}_j^{\dagger}|\skew{2.4}\tilde{\Psi}_j\rangle$ to have no excitations above j. Thus, we consider $\varepsilon(t) = \varepsilon^{(0)}(t) + \varepsilon^{(1)}(t)$, where $\varepsilon^{(1)}(t)$ are higher frequency components that cancel out higher order transitions, making $\hat{{\cal U}}_{01}(t_g) \approx \hat{{\cal U}}_{01}^{(0)}(t_g)$.

The error in the transmon gate is given by:

$$\begin{align} \hat{{\cal U}}_{\textrm{err}}(t) = \hat{{\cal U}}_{01}^{(0),\dagger}(t)\hat{{\cal U}}_{01}(t) = \mathrm{Texp}\left[ i\int_0^td\tau\ \delta\hat{h}(\tau) - \textrm{h.c}. \right], \end{align} \tag{ 35 }$$

where,

$$\begin{align} \delta\hat{h}(\tau) = \hat{{\cal U}}_{01}^{(0),\dagger}(\tau)\left(\varepsilon(\tau)\sum_{n = 0}^{\infty}\sqrt{n + 1}\ e^{- inA\tau}|n + 1\rangle\langle n| - \varepsilon^{(0)}(\tau)|1\rangle\langle 0|\right)\hat{{\cal U}}_{01}^{(0)}(\tau). \end{align} \tag{ 36 }$$

$\hat{{\cal U}}_{\textrm{err}}$ cannot be found analytically. To linear order in the perturbation, we ignore the corrections due to non-commutativity, and impose:

$$\begin{align} \int_0^{t_g}d\tau\ \delta\hat{h}(\tau) = \hat{0}. \end{align} \tag{ 37 }$$

This in turn implies the following set of equations that the drive amplitude $\varepsilon(t)$ must satisfy:

$$\begin{align} \int_0^{t_g}\left(\varepsilon(\tau) - \varepsilon^{(0)}(\tau)\right)d\tau & = 0, \end{align} \tag{ 38 }$$

$$\begin{align} \int_0^{t_g}\varepsilon(\tau)e^{iA\tau\pm i\Theta(t)}d\tau & = 0, \end{align} \tag{ 39 }$$

$$\begin{align} n\geqslant2:\ \int_0^{t_g}\varepsilon(\tau)e^{inA\tau}d\tau & = 0. \end{align} \tag{ 40 }$$

The first condition is a no-bias condition ensuring that $\int_0^{t_g}\varepsilon(\tau)d\tau = \Theta$. The second condition, with $\Theta(t) = \int_0^td\tau\varepsilon^{(0)}(\tau)$, ensures that the transitions $|1\rangle\leftrightarrow|2\rangle$ cancel out at the end of the gate. Similarly, the transitions $|n\rangle\leftrightarrow|n + 1\rangle$ are cancelled by the third condition. In our simulations, we find that for typical parameters, the population beyond level 2 is negligible, so we can ignore the conditions for $n\geqslant3$.

The above conditions leave sufficient freedom in choosing $\varepsilon(t)$. In what follows, we choose $t_g = 3\pi/A$ and:

$$\begin{align} \varepsilon(t) = \sum_{\mu = 0}^4\varepsilon_{\mu}\sin\left(\frac{\left(2\mu + 1\right)\pi t}{t_g}\right), \end{align} \tag{ 41 }$$

where the parameters ε_µ are governed by an underdetermined set of linear equations that can be found from equations (38)–(40). Given $\varepsilon(t)$, we can find φ_app by minimizing $\left|\left\langle\, j,1\left|\hat{{\cal Q}}_{01}^{\dagger}\right|\skew{2.4}\tilde{\Psi}_j\right\rangle\right|$.

We point out that although we choose a linear approximation, the corrections are not always small. However, in our simulations we find that the error matrix $\hat{{\cal U}}_{\textrm{err}}$ is dominated by extraneous phases for all transmon levels. One might think that these can be removed by adding higher order corrections, however, we observed that such higher order corrections typically shifts the errors from $|2\rangle$ to $|n\geqslant3\rangle$ instead of removing them. Thus, we account for these extra phases by modifying the next steps of the induction.

If $\left|\left\langle\, j - 1,2\left|\hat{{\cal Q}}_{01}^{\dagger}\right|\skew{2.4}\tilde{\Psi}_j\right\rangle\right|$ is significant (we take the threshold as 0.1 in our simulations), a second gate between $|1\rangle$ and $|2\rangle$, say $\hat{{\cal Q}}_{12}$, needs to be applied. The correction method is analogous to above. The levels $|n\geqslant3\rangle$ are found to be approximately unoccupied, so we do not need to apply more gates. The time required for $\hat{\cal Q}_j$ is either $3\pi/A$ or $6\pi/A$ depending on whether a second gate is required. If we want to treat the transmon as a two-level system, we would need gate times $\gt$7$\pi/A$ (found empirically) to avoid leakage, which is even longer than applying two gates as we do here.

To evaluate the above algorithm, we numerically solve the time evolution that implements the protocol for the density matrix, using QuTiP [56, 57], satisfying:

$$\begin{align} \frac{d\hat{\rho}}{dt} = - \frac{i}{\hbar}\left[\skew{2.4}\hat{H}_\textrm{s} + \skew{2.4}\hat{H}_\textrm{d}(t),\hat{\rho}(t)]\right] + \sum_n\frac{1}{2}\left[2\hat{C}_n\hat{\rho}(t)\hat{C}_n^{\dagger} - \hat{\rho}(t)\hat{C}_n^{\dagger}\hat{C}_n - \hat{C}_n^{\dagger}\hat{C}_n\hat{\rho}(t)\right], \end{align} \tag{ 42 }$$

where the Hamiltonian is given by equations (20) and (21). Assuming zero temperature, we take the Lindblad operators pertaining to magnon decay $\hat{C}_1 = \sqrt{\kappa_m}\skew{2.4}\hat{m}$, transmon decay $\hat{C}_2 = \sqrt{\kappa_q}\skew{2.4}\hat{q}$ and transmon dephasing $\hat{C}_3 = \sqrt{\gamma_q}\skew{2.4}\hat{q}^{\dagger}\skew{2.4}\hat{q}$. We did not include magnon dephasing as it has not yet been experimentally observed, although there are theoretical reasons for pure dephasing of magnons [58]. We choose typical parameters $\omega_m = 2\pi\times6$ GHz, $\omega_q = 2\pi\times 5$ GHz, $\kappa_m = 10^{-4}\omega_m$, $\kappa_q = \gamma_q = 10^{-4}\omega_q$, $A = 2\pi\times 300$ MHz, and coupling $g = 2\pi\times25$ MHz. For a given target state $|\psi_T\rangle$, the algorithm described in the previous sections gives the signals $\Delta[\psi_T](t)$ and $\skew{2.4}\tilde{\varepsilon}[\psi_T](t)$, along with the protocol time $T_{\textrm{tot}}[\psi_T]$. Starting from the ground state, $|0,0\rangle$, we numerically solve the above Lindblad equation with the given signals, Δ and $\skew{2.4}\tilde{\varepsilon}$, and find the reduced density matrix for the magnon, $\hat{\rho}_{\textrm{mag}}[\psi_T] = \textrm{Tr}_{\skew{2.4}\hat{q}}\left[\hat{\rho}(T_\textrm{tot}[\psi_T])\right]$:

$$\begin{align} |\psi_T\rangle \xrightarrow{\mathrm{protocol}} \{\Delta(t),\skew{2.4}\tilde{\varepsilon}(t)\} \xrightarrow{\mathrm{simulation}} \hat{\rho}_{\textrm{mag}}. \end{align} \tag{ 43 }$$

Finally, we calculate the fidelity of state generation $F[\psi_T] = \left\langle\psi_T\left|\hat{\rho}_{\textrm{mag}}\left[\psi_T\right]\right|\psi_T\right\rangle$ which can be interpreted as the probability of finding the magnetization in the target state $|\psi_T\rangle$.

We compare the Wigner function of the target state $|\psi_T\rangle$ and the achieved state $\hat{\rho}_{\textrm{mag}}[\psi_T]$ for some special cases. The Wigner function is defined as:

$$\begin{align} W(\beta,\hat{\rho}) = \int\frac{d^2\gamma}{\pi}\textrm{Tr}\left[\hat{\rho}\exp\left[\left(\skew{2.4}\hat{m} - \beta\right)^{\dagger}\gamma - \gamma^*\left(\skew{2.4}\hat{m} - \beta\right)\right]\right]. \end{align} \tag{ 44 }$$

The Wigner function can be interpreted as follows. Define a symmetric product $\left\{\skew{2.4}\hat{m}^{\dagger,p}\skew{2.4}\hat{m}^q\right\}_S$ as an average over all permutations of p $\skew{2.4}\hat{m}^{\dagger}$'s and q $\skew{2.4}\hat{m}$'s. For example,

$$\begin{align} \left\{\skew{2.4}\hat{m}^{\dagger,2}\skew{2.4}\hat{m}\right\}_S = \frac{\skew{2.4}\hat{m}^{\dagger,2}\skew{2.4}\hat{m} + \skew{2.4}\hat{m}^{\dagger}\skew{2.4}\hat{m}\skew{2.4}\hat{m}^{\dagger} + \skew{2.4}\hat{m}\skew{2.4}\hat{m}^{\dagger,2}}{3}. \end{align} \tag{ 45 }$$

In a given state $\hat{\rho}$, the average of the above operator is given by [59, 60]:

$$\begin{align} \textrm{Tr}\left[\hat{\rho}\left\{\skew{2.4}\hat{m}^{\dagger,p}\skew{2.4}\hat{m}^q\right\}_S\right] = \int\frac{d^2\beta}{\pi}\beta^{*,p}\beta^qW(\beta,\hat{\rho}). \end{align} \tag{ 46 }$$

We can convert this into expectation values of the magnetization operators $\skew{2.4}\hat{M}_x$ and $\skew{2.4}\hat{M}_y$ via the substitution $2\mathcal{M}_\textrm{ZPF}\skew{2.4}\hat{m}\rightarrow\skew{2.4}\hat{M}_x - i\skew{2.4}\hat{M}_y$ (note that we are using the dressed operator $\skew{2.4}\hat{m}$ instead of the bare magnetization operator $\skew{2.4}\hat{M}$ as the distinction between them is small) getting:

$$\begin{align} \textrm{Tr}\left[\hat{\rho}\left\{\skew{2.4}\hat{M}_x^p\skew{2.4}\hat{M}_y^q\right\}_S\right] = \int\frac{dM_xdM_y}{\pi}M_x^pM_y^q\ \skew{2.4}\tilde{W}\left(M_x,M_y,\hat{\rho}\right), \end{align} \tag{ 47 }$$

where $\skew{2.4}\tilde{W}\left(M_x,M_y . \hat{\rho}\right) = W\left(\frac{M_x - iM_y}{2\mathcal{M}_\textrm{ZPF}},\hat{\rho}\right)$. The above equation suggests an interpretation of $\skew{2.4}\tilde{W}$ as a quasi-classical probability distribution. However, note that $\skew{2.4}\tilde{W}$ is not always positive and a negative $\skew{2.4}\tilde{W}$ implies a lack of a probabilistic interpretation considered to be a signature of non-classicality. The unit $\mathcal{M}_\textrm{ZPF}$, defined in equation (2), is given by:

$$\begin{align} \frac{\mathcal{M}_\textrm{ZPF}}{M_s} = 8\times10^{- 6}\sqrt{\frac{1~\mu\mathrm{m}^{3}}{V_m}}, \end{align} \tag{ 48 }$$

where V_m is the volume of the magnet. This is a very small angle but it can be measured via the magnet's coupling to the transmon [34, 35].

In figure 4, we consider the case of Fock states $|\psi_{\textrm{Fock}}(n)\rangle = |n\rangle$ and their superposition with the ground state $|\psi_{\textrm{Sup}}(n)\rangle\propto|0\rangle + |n\rangle$ where $n\in\{1,2,\dots,14\}$. The fidelity decreases with increasing n. This is caused by two factors. The dominant factor is an increase in the protocol time with the number of magnons leading to higher decay. Secondly, each gate introduces a small amount of error, including residual transfer to higher levels, that compounds with the number of gates, linearly increasing with n. It appears that Fock states are more difficult to generate than the superposition states, which could be due to a larger average number of magnons in $|\psi_{\textrm{Fock}}\rangle$, compared to $|\psi_{\textrm{Sup}}\rangle$, leading to a faster decay. For $n\in\{6,10\}$ and $T\in\{\textrm{Fock},\textrm{Sup}\}$, we also show Wigner functions of the target state, $\skew{2.4}\tilde{W}\left(M_x,M_y,|\psi_T(n)\rangle\langle \psi_T(n)|\right)$ and the achieved state $\skew{2.4}\tilde{W}\left(M_x,M_y,\hat{\rho}_{\textrm{mag}}\left[\psi_T(n)\right]\right)$. The essential feature of $|\psi_{\textrm{Sup}}(n)\rangle$ are the azimuthal fringes of order n, which is well captured for n = 6 and weakly reproduced for n = 10. The essential feature of $|\psi_{\textrm{Fock}}(n)\rangle$ are the n circular fringes, that are destroyed by dephasing.

Zoom In Zoom Out Reset image size

A superposition of two semi-classical states, known as a cat state, is useful for error correction protocols because of their insensitivity to particle loss noise [61], particularly as carriers of information [62–64] and in metrology [65–68]. In case of a magnet, we can define semi-classical states as the minimum-uncertainty states with isotropic quantum fluctuations, given by coherent states [69]:

$$\begin{align} |\xi\rangle = e^{- \left|\xi\right|^2/2}\sum_{n = 0}^{\infty}\frac{\xi^n}{\sqrt{n!}}|n\rangle, \end{align} \tag{ 49 }$$

for an arbitrary complex number ξ. A coherent state has average $\left\langle \skew{2.4}\hat{m} \right\rangle = \xi$ and can be interpreted as the quantum state 'closest' to the classical state with magnetization $M_x = 2\mathcal{M}_\textrm{ZPF}\mathrm{Re}\left(\xi\right)$ and $M_y = - 2\mathcal{M}_\textrm{ZPF}\mathrm{Im}\left(\xi\right)$. An 'odd' cat state given by $|\psi_{\textrm{odd}}(\xi)\rangle\propto|\xi/2\rangle - |- \xi/2\rangle$ has two semi-classical components separated by the magnetization $2\mathcal{M}_\textrm{ZPF}\left|\xi\right|$; the same holds for an 'even' cat state $|\psi_{\textrm{even}}(\xi)\propto|\xi/2\rangle + |- \xi/2\rangle\rangle$. Figure 5 shows the fidelity of generation as a function of ξ. The Wigner function for an ideal cat state has two blobs at the expected peaks along with fringes in the centre. The presence of the fringes signify the quantum coherence between the two semi-classical components. In contrast, the Wigner function for a classical distribution, $\hat{\rho}_{\textrm{cl}}(\xi)\propto|\xi/2\rangle\langle \xi/2| + |- \xi/2\rangle\langle - \xi/2|$ would still have the blobs but no fringes.

Zoom In Zoom Out Reset image size

In the achieved states, we find such fringes for cat sizes up to 4 until which the fidelities are not too low. At cat size 5, we see that the achieved state are closer to a classical probability distribution (they have $\gt$0.85 fidelity with the corresponding classical states $\hat{\rho}_{\textrm{cl}}(\xi = 5)$). Differences in the odd and cat states for low sizes can be explained on the basis that odd cats have larger number of magnons, $\frac{\xi^2}{4}\frac{1 + e^{- \xi^2/2}}{1 - e^{- \xi^2/2}}$, than even cats, $\frac{\xi^2}{4}\frac{1 - e^{- \xi^2/2}}{1 + e^{- \xi^2/2}}$, an effect that is negligible for ξ > 2. Fidelity for generation of any other state can be found using the publicly available codes [70].

Here, we focussed on adapting and characterizing the protocol for a system with a magnet and a transmon. There are several methods by which the protocol can be improved, as we now elaborate.

5.1. Pre-processing

We considered only a specific type of operations, where the transmon is excited and the excitation is transferred to the magnons. With the help of the MW cavity, one can also apply other operations, in particular, a displacement given by $\skew{2.4}\hat{D}(\beta) = e^{\beta\skew{2.4}\hat{m}^{\dagger} - \beta^*\skew{2.4}\hat{m}}$ caused by a coherent excitation of the magnet by MWs. This can be used to decrease the time of the protocol as one can choose an appropriate β, such that $\skew{2.4}\hat{D}^{\dagger}(\beta)|\psi_T\rangle$ has fewer magnons than $|\psi_T\rangle$. For example, if $\left\langle\psi_T\left|\skew{2.4}\hat{m}\right|\psi_T\right\rangle = C$, then:

$$\begin{align} \left\langle\psi_T\left|\skew{2.4}\hat{D}(C)\skew{2.4}\hat{m}^{\dagger}\skew{2.4}\hat{m}\skew{2.4}\hat{D}^{\dagger}(C)\right|\psi_T\right\rangle = \left\langle\psi_T\left|\skew{2.4}\hat{m}^{\dagger}\skew{2.4}\hat{m}\right|\psi_T\right\rangle - \left|C\right|^2, \end{align} \tag{ 50 }$$

implying that $\skew{2.4}\hat{D}^{\dagger}(C)|\psi_T\rangle$ has fewer average magnons than $|\psi_T\rangle$. A better β can be found numerically that decreases the highest amplitude in $|\psi_T\rangle$ to within some given tolerance. Then, we run the protocol to create $\skew{2.4}\hat{D}^{\dagger}(\beta)|\psi_T\rangle$ and apply a MW pulse to displace keeping a high fidelity. Pre-processing with more complicated analog signals can decrease the duration further [42].

5.2. High fidelity gates

5.3. Combining steps

We proposed and numerically demonstrated a protocol to generate an arbitrary quantum state of the magnetization using its coupling to a superconducting transmon, which is a weakly anharmonic oscillator. The protocol involves repeatedly applying transmon gates and transferring an excitation to the magnons. As longer protocol times lead to higher dissipation, we can improve fidelity of generation by reducing the gate time, t_g , and the transfer time, t_on, in the notation of sections 2 and 3. Fast gate pulses with duration ∼2$\pi/A$ (A is defined in equation (19)) would imply leakage of excitation into higher transmon levels that can be reduced by modifying the input pulses [51, 71, 73, 74]. We use correction pulses by removing errors at first order in a Magnus expansion, similar to [54, 55], bringing the gate time down to $3\pi/A$. Correcting for second order term in the expansion does not seem to decrease errors, but rather transfer the population to higher levels. The transfer time, t_on, can be reduced with a higher magnon-transmon coupling. As the magnon-MW coupling can reach very high values, it is likely to be experimentally feasible to have high magnon-transmon couplings. However, a coupling $g\sim A$ causes leakage of magnons into higher transmon levels leading to loss of fidelities, so we choose $g\lt0.1A$. For such couplings, the leakage is small and can be incorporated by renormalizing the duration of the coupling, $t_{\textrm{on},j}$.

The main source of errors is the magnon dissipation putting a limit on the size of the states. We find that the protocol faithfully produces states up to a largest occupation of ∼9 magnons and an average of ∼4 magnons. We note that our protocol creates states in dressed magnons, given by $\skew{2.4}\hat{m}$ instead of bare magnons $\skew{2.4}\hat{M}$. As discussed in appendix A, for a sufficiently large detuning between MWs and magnons, the distinction is small. We demonstrated our protocol for cases of Fock states $|n\rangle$, Fock superpositions $|0\rangle + |n\rangle$, and cat states $|\beta\rangle\pm|-\beta\rangle$. While our protocol gives good fidelities for states of moderate sizes, we expect that advanced numerical techniques such as optimal quantum control [74, 75, 78] can improve the state generation (likely by reducing the protocol time) by using our protocol as an initial guiding guess.

With a setup to generate an arbitrary magnetization, there still lies a question of how to experimentally measure it, which we leave to future work. The figures can be generated via the publicly available codes [70].

The authors thank H Ribeiro for useful discussions. We acknowledge funding from the Max Planck Society and from the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) through Project-ID 429529648-TRR 306 QuCoLiMa ('Quantum Cooperativity of Light and Matter').

The codes to generate the plots in the paper are available at the following URL/DOI: https://github.com/sancharsharma/Magnon-State-Generation.

In this appendix, we derive the effective Hamiltonian for a magnet coupled to a transmon indirectly via a MW cavity.

A.1. Microwaves

The MW cavity is modelled by the classical Hamiltonian density:

$$\begin{align} \mathcal{H}_\textrm{mw} = \frac{\epsilon_0\left|\boldsymbol{E}\right|^2 + \mu_0\left|\boldsymbol{H}\right|^2}{2}. \end{align} \tag{ A.1 }$$

It contains several modes labelled by r giving the quantization:

$$\begin{align} \boldsymbol{H}(\boldsymbol{r})\rightarrow\sum_r\boldsymbol{H}_r(\boldsymbol{r})\skew{2.4}\hat{A}_r + \boldsymbol{H}_r^*(\boldsymbol{r})\skew{2.4}\hat{A}_r^{\dagger}, \end{align} \tag{ A.2 }$$

and similar holds for E . Here, H _r and $\skew{2.4}\hat{A}_r$ are respectively the mode profile and the annihilation operator of mode r. With appropriate normalization, this gives the Hamiltonian,

$$\begin{align} \skew{2.4}\hat{H}_\textrm{mw} = \hbar\sum_r\Omega_{a,r}\skew{2.4}\hat{A}_r^{\dagger}\skew{2.4}\hat{A}_r, \end{align} \tag{ A.3 }$$

where $\Omega_{a,r}$ is the resonance frequency of mode r.

A.2. Magnet

If the magnet is much smaller than MW wavelengths (∼cm), then only the total magnetization M couples to the MWs. Under a DC field along $-\boldsymbol{z}$ with magnitude H_app, its classical Hamiltonian density is:

$$\begin{align} \mathcal{H}_{\textrm{m}, \text{mw}} = \mu_0H_{\textrm{app}}M_z - \mu_0\boldsymbol{M}\cdot\boldsymbol{H}, \end{align} \tag{ A.4 }$$

where H is the magnetic field due to the cavity. The magnetization is quantized by the Holstein-Primakoff transformation [45, 46], given by a form:

$$\begin{align} M_x - iM_y\rightarrow2\mathcal{M}_\textrm{ZPF}\skew{2.4}\hat{M}\sqrt{\frac{1 - \beta\skew{2.4}\hat{M}^{\dagger}\skew{2.4}\hat{M}}{1 - \beta}}, \end{align} \tag{ A.5 }$$

where $\mathcal{M}_\textrm{ZPF}$ and β are to be found. We know that a spin in an applied field undergoes a Larmor precession with frequency $\Omega_m = \gamma\mu_0H_{\textrm{app}}$ where γ is the absolute value of the gyromagnetic ratio. This implies that the Kittel mode should have the resonance frequency, i.e. after quantization:

$$\begin{align} \mu_0V_mH_{\textrm{app}}M_z\rightarrow\hbar\Omega_m\skew{2.4}\hat{M}^{\dagger}\skew{2.4}\hat{M}. \end{align} \tag{ A.6 }$$

This implies $M_z = -\skew{2.4}\tilde{M}_s + M_{\textrm{sing}}M^{\dagger}M$ where $\skew{2.4}\tilde{M}_s$ (to be found) is the saturation magnetization reduced by zero point fluctuations and $M_{\textrm{sing}} = \gamma\hbar/V_m$ is the magnetization reduction because of one (delocalized) spin flip. Using $M_x^2 + M_y^2 + M_z^2 = M_s^{\,2}$, we find the relations:

$$\begin{align} \skew{2.4}\tilde{M}_s^2 + 2\mathcal{M}_\textrm{ZPF}^2 = M_s^{\,2},\ \mathcal{M}_\textrm{ZPF} = \sqrt{\frac{\skew{2.4}\tilde{M}_sM_{\textrm{sing}}}{2}},\ M_{\textrm{sing}}^{\,2} = \frac{4\mathcal{M}_\textrm{ZPF}^2\beta}{1 - \beta}. \end{align} \tag{ A.7 }$$

This is solved explicitly by:

$$\begin{align} \skew{2.4}\tilde{M}_s = \sqrt{M_s^{\,2} + \frac{M_{\textrm{sing}}^{\,2}}{4}} - \frac{M_{\textrm{sing}}}{2},\ \beta = \frac{M_{\textrm{sing}}}{2\skew{2.4}\tilde{M}_s + M_{\textrm{sing}}} \end{align} \tag{ A.8 }$$

Typically $M_{\textrm{sing}}/M_s\lt10^{- 10}$, so we can resort to the leading order corrections,

$$\begin{align} \skew{2.4}\tilde{M}_s \approx M_s - \frac{M_{\textrm{sing}}}{2},\ \mathcal{M}_\textrm{ZPF} \approx \sqrt{\frac{M_sM_{\textrm{sing}}}{2}},\ \beta \approx \frac{M_{\textrm{sing}}}{2M_s}. \end{align} \tag{ A.9 }$$

For $\left\langle \skew{2.4}\hat{M}^{\dagger}\skew{2.4}\hat{M} \right\rangle\ll1/\beta\sim10^{10}$, we can simply substitute:

$$\begin{align} M_x - iM_y\rightarrow2\mathcal{M}_\textrm{ZPF}\skew{2.4}\hat{M},\ M_z\rightarrow -M_s + \frac{2\mathcal{M}_\textrm{ZPF}^2}{M_s}\skew{2.4}\hat{M}^{\dagger}\skew{2.4}\hat{M}. \end{align} \tag{ A.10 }$$

This gives the magnon Hamiltonian (coupled to the cavity):

$$\begin{align} \skew{2.4}\hat{H}_{\textrm{m}, \text{mw}} = \hbar\Omega_m\skew{2.4}\hat{M}^{\dagger}\skew{2.4}\hat{M} + \sum_rG_{m,r}\skew{2.4}\hat{M}^{\dagger}\skew{2.4}\hat{A}_r + G_{m,r}^*\skew{2.4}\hat{M}\skew{2.4}\hat{A}_r^{\dagger}, \end{align} \tag{ A.11 }$$

where $\Omega_m = \gamma\mu_0H_{\textrm{app}}$ and:

$$\begin{align} \hbar G_{m,r} = - \mu_0\mathcal{M}_\textrm{ZPF}\left(H_{x,r} - iH_{y,r}\right). \end{align} \tag{ A.12 }$$

We assumed rotating wave approximation which is valid when $G\ll\Omega_m + \Omega_a$.

A.3. Transmon

We consider a flux-tunable transmon that consists of two Josephson junctions in parallel forming a SQUID loop. Its Hamiltonian is [48, 49]:

$$\begin{align} \skew{2.4}\hat{H}_{\textrm{tr}, \text{mw}} = E_C\left(\hat{n} + \hat{n}_{\textrm{mw}}\right)^2 - E_{J1}\cos\hat{\phi}_1 - E_{J2}\cos\hat{\phi}_2. \end{align} \tag{ A.13 }$$

Here, $\hat{n}$ is the number of quanta exchanged between the two junctions, $\hat{n}_{\textrm{mw}}$ is the number of quanta excited by MWs, and $\hat{\phi}_{1,2}$ are the fluxes through each junction. Under weak coupling, we have $\hat{n}_{\textrm{mw}}\propto\left|\boldsymbol{E}\right|$. Flux quantization constraints $\hat{\phi}_1 - \hat{\phi}_2 = 2\pi m + 2\pi\Phi/\Phi_0$ where m is an arbitrary integer, Φ is the flux through the SQUID loop, and $\Phi_0 = h/2e$ is the superconducting flux quantum. Defining $2\hat{\phi} = \hat{\phi}_1 + \hat{\phi_2}$,

$$\begin{align} \skew{2.4}\hat{H}_\textrm{tr} = E_C\hat{n}^2 - E_{J,\textrm{eff}}\cos\left(\hat{\phi} + \phi_0\right). \end{align} \tag{ A.14 }$$

where the phase is:

$$\begin{align} \phi_0 = \tan^{- 1}\left(\frac{E_{J1} - E_{J2}}{E_{J1} + E_{J2}}\tan\frac{\pi\Phi}{\Phi_0}\right) + n\pi \end{align} \tag{ A.15 }$$

and the effective junction energy is:

$$\begin{align} E_{J,\textrm{eff}}(\Phi) = \sqrt{\left(E_{J1} + E_{J2}\right)^2\cos^2\frac{\pi\Phi}{\Phi_0} + \left(E_{J1} - E_{J2}\right)^2\sin^2\frac{\pi\Phi}{\Phi_0}}. \end{align} \tag{ A.16 }$$

In the transmon limit $E_{J,\textrm{eff}}\gg E_C$, the fluctuations in $\hat{\phi} + \phi_0$ are small, so we can expand the cosine and get an approximately quadratic Hamiltonian. This suggests the quantization,

$$\begin{align} \hat{n} = \left(\frac{E_{J,\textrm{eff}}}{2E_C}\right)^{1/4}\frac{\skew{2.4}\hat{Q} - \skew{2.4}\hat{Q}^{\dagger}}{\sqrt{2}i},\ \hat{\phi} + \phi_0 = \left(\frac{2E_C}{E_{J,\textrm{eff}}}\right)^{1/4}\frac{\skew{2.4}\hat{Q} + \skew{2.4}\hat{Q}^{\dagger}}{\sqrt{2}}, \end{align} \tag{ A.17 }$$

that reduces the Hamiltonian to:

$$\begin{align} \skew{2.4}\hat{H}_\textrm{q,mw} = \hbar\Omega_q\skew{2.4}\hat{Q}^{\dagger}\skew{2.4}\hat{Q} - \hbar\frac{A}{2}\skew{2.4}\hat{Q}^{\dagger}\skew{2.4}\hat{Q}\left(\skew{2.4}\hat{Q}^{\dagger}\skew{2.4}\hat{Q} - 1\right) + \hbar\sum_r\left(G_{q,r}\skew{2.4}\hat{Q}^{\dagger}\skew{2.4}\hat{A}_r + G_{q,r}^*\skew{2.4}\hat{Q}\skew{2.4}\hat{A}_r^{\dagger}\right), \end{align} \tag{ A.18 }$$

ignoring higher order terms. Here $\hbar\Omega_q = \sqrt{2E_CE_{J,\textrm{eff}}} - E_C/4$ and $\hbar A = E_C/4$.

A.4. Dressed magnons and transmons

The dressed states are found by diagonalizing the total Hamiltonian $\skew{2.4}\hat{H} = \skew{2.4}\hat{H}_\textrm{0} + \skew{2.4}\hat{H}_\textrm{int}$ where:

$$\begin{align} \frac{\skew{2.4}\hat{H}_\textrm{0}}{\hbar} = \sum_r \left( \Omega_{a,r}\skew{2.4}\hat{A}_r^{\dagger}\skew{2.4}\hat{A}_r + \Omega_m\skew{2.4}\hat{M}^{\dagger}\skew{2.4}\hat{M} + \Omega_q\skew{2.4}\hat{Q}^{\dagger}\skew{2.4}\hat{Q} - \frac{A}{2}\skew{2.4}\hat{Q}^{\dagger}\skew{2.4}\hat{Q}\left(\skew{2.4}\hat{Q}^{\dagger}\skew{2.4}\hat{Q} - 1\right) \right), \end{align} \tag{ A.19 }$$

and the interaction:

$$\begin{align} \frac{\skew{2.4}\hat{H}_\textrm{int}}{\hbar} = \sum_r\left(G_{m,r}\skew{2.4}\hat{M}^{\dagger}\skew{2.4}\hat{A}_r + G_{m,r}^*\skew{2.4}\hat{M}\skew{2.4}\hat{A}_r^{\dagger} + G_{q,r}\skew{2.4}\hat{Q}^{\dagger}\skew{2.4}\hat{A}_r + G_{q,r}^*\skew{2.4}\hat{Q}\skew{2.4}\hat{A}_r^{\dagger}\right). \end{align} \tag{ A.20 }$$

However, exact diagonalization is not possible analytically, so we resort to a perturbative treatment, in particular a Schrieffer–Wolff (SW) transformation [79]. The SW transformation perturbatively gives the rotation operator $\skew{2.4}\hat{S} = - \skew{2.4}\hat{S}^{\dagger}$ such that $e^{- \skew{2.4}\hat{S}}\skew{2.4}\hat{H}e^{\skew{2.4}\hat{S}}$ is diagonal to a given order. We find $\skew{2.4}\hat{S}$ to linear order in coupling constants G that gives the rotated Hamiltonian,

$$\begin{align} e^{- \skew{2.4}\hat{S}}\skew{2.4}\hat{H}e^{\skew{2.4}\hat{S}} = \skew{2.4}\hat{H}_\textrm{0} + \skew{2.4}\hat{H}_\textrm{int} + \left[\skew{2.4}\hat{H}_\textrm{0},\skew{2.4}\hat{S}\right] + \left[\skew{2.4}\hat{H}_\textrm{int},\skew{2.4}\hat{S}\right] + \frac{1}{2}\left[\left[\skew{2.4}\hat{H}_\textrm{0},\skew{2.4}\hat{S}\right],\skew{2.4}\hat{S}\right] + O\left(G^3\right). \end{align} \tag{ A.21 }$$

To remove first order corrections, we choose:

$$\begin{align} \skew{2.4}\hat{S} = \sum_r\left(\frac{G_{m,r}}{\Omega_{a,r} - \Omega_m}\skew{2.4}\hat{A}_r\skew{2.4}\hat{M}^{\dagger} + G_{q,r}\chi_r\left[\skew{2.4}\hat{Q}^{\dagger}\skew{2.4}\hat{Q}\right]\skew{2.4}\hat{A}_r\skew{2.4}\hat{Q}^{\dagger}\right) - \textrm{h.c.}, \end{align} \tag{ A.22 }$$

where the transmon susceptibility is:

$$\begin{align} \chi_r\left[x\right] = \frac{1}{\Omega_{a,r} - \Omega_q + A\left(x - 1\right)}, \end{align} \tag{ A.23 }$$

ensuring $\skew{2.4}\hat{H}_\textrm{int} + \left[\skew{2.4}\hat{H}_\textrm{0},\skew{2.4}\hat{S}\right] = \hat{0}$. Note that we need $\left\Vert \skew{2.4}\hat{S}\right\Vert \ll1$ that implies a limit on the number of excitations $\{N_{a,r},N_m,N_q\}$ s.t. for all r,

$$\begin{align} G_{m,r}\sqrt{N_{a,r}N_m}\ll\Omega_{a,r} - \Omega_m,\ G_{q,r}\sqrt{N_{a,r}N_q}\ll\Omega_{a,r} - \Omega_q, \end{align} \tag{ A.24 }$$

are satisfied. This gives the rotated Hamiltonian,

$$\begin{align} e^{- \skew{2.4}\hat{S}}\skew{2.4}\hat{H}e^{\skew{2.4}\hat{S}} = \skew{2.4}\hat{H}_\textrm{0} + \frac{1}{2}\left[\skew{2.4}\hat{H}_\textrm{int},\skew{2.4}\hat{S}\right] + O\left(G^3\right). \end{align} \tag{ A.25 }$$

As there are four terms in both $\skew{2.4}\hat{H}_\textrm{int}$ and $\skew{2.4}\hat{S}$, a total of 16 commutators are required. If we ignore all non-resonant terms, such as $\skew{2.4}\hat{A}_r^2\skew{2.4}\hat{Q}^{\dagger,2}$, and use the Hermiticity of $\skew{2.4}\hat{H}_\textrm{int}$ we are left with:

$$\begin{align} \frac{1}{2\hbar}\left[\skew{2.4}\hat{H}_\textrm{int},\skew{2.4}\hat{S}\right] = \sum_{rs}\left[\left(G_{m,r}^*\skew{2.4}\hat{M} + G_{q,r}^*\skew{2.4}\hat{Q}\right)\skew{2.4}\hat{A}_r^{\dagger},\left(\frac{G_{m,s}}{\Omega_{a,s} - \Omega_m}\skew{2.4}\hat{M}^{\dagger} + G_{q,s}\chi_s\left(\skew{2.4}\hat{Q}^{\dagger}\skew{2.4}\hat{Q}\right)\skew{2.4}\hat{Q}^{\dagger}\right)\skew{2.4}\hat{A}_s\right] + h . c. \end{align} \tag{ A.26 }$$

For $\skew{2.4}\hat{f},\skew{2.4}\hat{g}$ commuting with all photon operators,

$$\begin{align} \left[\skew{2.4}\hat{f}\skew{2.4}\hat{A}_r^{\dagger},\skew{2.4}\hat{g}\skew{2.4}\hat{A}_s\right] = - \skew{2.4}\hat{f}\,\skew{2.4}\hat{g}\delta_{rs} + \left[\skew{2.4}\hat{f},\skew{2.4}\hat{g}\right]\skew{2.4}\hat{A}_r^{\dagger}\skew{2.4}\hat{A}_s. \end{align} \tag{ A.27 }$$

This leads to

$$\begin{align} \frac{1}{2\hbar}\left[\skew{2.4}\hat{H}_\textrm{int},\skew{2.4}\hat{S}\right] & = - \sum_r\left(G_{m,r}^*\skew{2.4}\hat{M} + G_{q,r}^*\skew{2.4}\hat{Q}\right)\left(\frac{G_{m,r}}{\Omega_{a,r} - \Omega_m}\skew{2.4}\hat{M}^{\dagger} + G_{q,r}\chi_r\left(\skew{2.4}\hat{Q}^{\dagger}\skew{2.4}\hat{Q}\right)\skew{2.4}\hat{Q}^{\dagger}\right) \nonumber\\ & \quad + \, \delta\Omega_{a,r}\left[\skew{2.4}\hat{Q}^{\dagger}\skew{2.4}\hat{Q}\right]\skew{2.4}\hat{A}_r^{\dagger}\skew{2.4}\hat{A}_r + h . c . \end{align} \tag{ A.28 }$$

where

$$\begin{align} \delta\Omega_{a,r}\left[x\right] = \frac{\left|G_{m,r}\right|^2}{\Omega_{a,r} - \Omega_m} + \left|G_{q,r}\right|^2\left(\left(x + 1\right)\chi_r\left(x + 1\right) - x\chi_r\left(x\right)\right). \end{align} \tag{ A.29 }$$

This procedure yields the rotated Hamiltonian. To return back to the original one, we define the dressed operators: $\skew{2.4}\hat{m} = e^{\skew{2.4}\hat{S}}\skew{2.4}\hat{M}e^{- \skew{2.4}\hat{S}}$, $\skew{2.4}\hat{q} = e^{\skew{2.4}\hat{S}}\skew{2.4}\hat{Q}e^{- \skew{2.4}\hat{S}}$, and $\hat{a}_r = e^{\skew{2.4}\hat{S}}\skew{2.4}\hat{A}_re^{- \skew{2.4}\hat{S}}$. Up to linear order in G, the dressed magnon operator is:

$$\begin{align} \skew{2.4}\hat{m} \approx \skew{2.4}\hat{M} - \sum\frac{G_{m,r}}{\Omega_{a,r} - \Omega_m}\skew{2.4}\hat{A}_r + O\left(G^2\right), \end{align} \tag{ A.30 }$$

and the dressed transmon operator is:

$$\begin{align} \skew{2.4}\hat{q} \approx \skew{2.4}\hat{Q} - \sum G_{q,r}\skew{2.4}\hat{A}_r\left(\chi_r\left[\skew{2.4}\hat{Q}\skew{2.4}\hat{Q}^{\dagger}\right] - \chi_r\left[\skew{2.4}\hat{Q}^{\dagger}\skew{2.4}\hat{Q}\right]\right) + O\left(G^2\right). \end{align} \tag{ A.31 }$$

This gives the final Hamiltonian:

$$\begin{align} \frac{\skew{2.4}\hat{H}}{\hbar} = \omega_m\skew{2.4}\hat{m}^{\dagger}\skew{2.4}\hat{m} + \left(\omega_q + \hat{\alpha}\left[\skew{2.4}\hat{q}^{\dagger}\skew{2.4}\hat{q}\right]\right)\skew{2.4}\hat{q}^{\dagger}\skew{2.4}\hat{q} + \skew{2.4}\hat{m}\skew{2.4}\hat{q}^{\dagger}\skew{2.4}\hat{g}\left[\skew{2.4}\hat{q}^{\dagger}\skew{2.4}\hat{q}\right] + \skew{2.4}\hat{g}\left[\skew{2.4}\hat{q}^{\dagger}\skew{2.4}\hat{q}\right]\skew{2.4}\hat{m}^{\dagger}\skew{2.4}\hat{q} + \frac{\skew{2.4}\hat{H}_\textrm{mw,eff}}{\hbar}. \end{align} \tag{ A.32 }$$

Here, $\skew{2.4}\hat{H}_\textrm{mw,eff}$ is the effective MW Hamiltonian that is not relevant for our discussion. The renormalized magnon's frequency is Stark shifted by the MWs:

$$\begin{align} \omega_m = \Omega_m - \sum_r\frac{\left|G_{m,r}\right|^2}{\Omega_{a,r} - \Omega_m}. \end{align} \tag{ A.33 }$$

Similarly, the transmon's frequency is:

$$\begin{align} \omega_q = \Omega_q - \sum_r\frac{\left|G_{q,r}\right|^2}{\Omega_{a,r} - \Omega_q}. \end{align} \tag{ A.34 }$$

The Stark shift of higher transmon levels renormalizes the anharmonicity to:

$$\begin{align} \hat{\alpha}\left[\skew{2.4}\hat{q}^{\dagger}\skew{2.4}\hat{q}\right] = \frac{A}{2}\left(\skew{2.4}\hat{q}^{\dagger}\skew{2.4}\hat{q} - 1\right) - \sum_r\left|G_{q,r}\right|^2\left[\frac{1}{\Omega_{a,r} - \Omega_q + A\left(\skew{2.4}\hat{q}^{\dagger}\skew{2.4}\hat{q} - 1\right)} - \frac{1}{\Omega_{a,r} - \Omega_q}\right]. \end{align} \tag{ A.35 }$$

Finally, the coupling of each transmon transition to the magnons is given by:

$$\begin{align} \skew{2.4}\hat{g}\left[\skew{2.4}\hat{q}^{\dagger}\skew{2.4}\hat{q}\right] = - \frac{1}{2}\sum_rG_{m,r}^*G_{q,r}\left(\frac{1}{\Omega_{a,r} - \Omega_q + A\skew{2.4}\hat{q}^{\dagger}\skew{2.4}\hat{q}} + \frac{1}{\Omega_{a,r} - \Omega_m}\right). \end{align} \tag{ A.36 }$$

In this appendix, we show how to obtain the periods required for transferring quanta from transmon to magnons. In what follows, $|m,q\rangle$ denotes a state with m magnons and transmon at qth level. Given a state $|\Psi\rangle$ having a maximum of j excitations, i.e. $\left\langle m,q\left|\Psi\right.\right\rangle = 0$ for $m + q \gt j$ and a tunable Hamiltonian:

$$\begin{align} \frac{\skew{2.4}\hat{H}(t)}{\hbar} = \omega_m\skew{2.4}\hat{m}^{\dagger}\skew{2.4}\hat{m} + \left(\omega_q + \Delta(t)\right)\skew{2.4}\hat{q}^{\dagger}\skew{2.4}\hat{q} - \sum_{n = 2}^{\infty}\alpha_n|n\rangle\langle n| + \sum_{n = 0}^{\infty}g_n\skew{2.4}\hat{m}^{\dagger}|n\rangle\langle n + 1| + g_n^*\skew{2.4}\hat{m}|n + 1\rangle\langle n|, \end{align} \tag{ B.1 }$$

we want to find $\{\Delta(t),t_{\textrm{tot}}\}$ s.t. $\left\langle\, j,0\left|\hat{{\cal T}}^{\dagger}\right|\Psi\right\rangle = 0$, where the time-evolution operator is:

$$\begin{align} \hat{{\cal T}} = \mathrm{Texp}\left[ - \frac{i}{\hbar}\int_0^{t_{\textrm{tot}}}\skew{2.4}\hat{H}(\tau)d\tau \right], \end{align} \tag{ B.2 }$$

with Texp being the time-ordered exponential operator. Here, $\alpha_n\equiv\left\langle n\left|\hat{\alpha}\right|n\right\rangle$ and $g_n\equiv \left\langle n\left|\skew{2.4}\hat{g}\right|n\right\rangle$ as derived in the previous section.

First, consider the ansatz $\Delta(t) = \Delta_0\left[u(t - t_{\textrm{on}}) - u(t)\right]$ with $\Delta_0 = \omega_m - \omega_q$, s.t. the magnons and transmons are resonant for a time-period t_on, effectively switching on the coupling. Adding an off-period of t_off afterwards, we get $t_{\textrm{tot}} = t_{\textrm{on}} + t_{\textrm{off}}$ and we want to find $\{t_{\textrm{on}},t_{\textrm{off}}\}$ to satisfy the above condition. We can derive the time-evolution unitary explicitly using $[\skew{2.4}\hat{H}(t),\hat{N}] = 0$ where the total number of excitations are given by:

$$\begin{align} \hat{N} = \skew{2.4}\hat{m}^{\dagger}\skew{2.4}\hat{m} + \skew{2.4}\hat{q}^{\dagger}\skew{2.4}\hat{q}. \end{align} \tag{ B.3 }$$

Thus, we divide the wave-function as:

$$\begin{align} |\Psi(t)\rangle = \sum_{N = 0}^{\infty}\sum_{s = 0}^N\xi_{Ns}(t)|N - s,s\rangle. \end{align} \tag{ B.4 }$$

During the off-time, we have the equations of motion, (using the Schrödinger equation $|\dot{\Psi}\rangle = - i\skew{2.4}\hat{H}|\Psi\rangle$):

$$\begin{align} \dot{\xi}_{Ns} = - i\left(\omega_m\left(N - s\right) + \omega_qs - \alpha_s\right)\xi_{Ns} - i\left[g_s\sqrt{N - s}\ \xi_{N,s + 1} + g_{s - 1}^*\sqrt{N - s + 1}\ \xi_{N,s - 1}\right]. \end{align} \tag{ B.5 }$$

This is a set of linear differential equations which can be solved by diagonalizing the dynamical matrix. For the case N = j, we can write:

$$\begin{align} \boldsymbol{\xi}_j(t_{\textrm{on}} + t_{\textrm{off}}) = V_{\textrm{off}}e^{- i\Lambda_{\textrm{off}}t_{\textrm{off}}}V_{\textrm{off}}^{\dagger}\boldsymbol{\xi}_j(t_{\textrm{on}}), \end{align} \tag{ B.6 }$$

where the vector $\boldsymbol{\xi}_j = \left(\xi_{j0}\ \xi_{j1}\ \dots\ \xi_{jj}\right)^T$, and the dynamical matrix is diagonalized by $\{- i\Lambda_{\textrm{off}},V_{\textrm{off}}\}$. Similar analysis with $\omega_q\rightarrow\omega_m$ works for the on-time,

$$\begin{align} \boldsymbol{\xi}_j\left(t_{\textrm{on}} + t_{\textrm{off}}\right) = V_{\textrm{off}}e^{- i\Lambda_{\textrm{off}}t_{\textrm{off}}}V_{\textrm{off}}^{\dagger}V_{\textrm{on}}e^{- i\Lambda_{\textrm{on}}t_{\textrm{on}}}V_{\textrm{on}}^{\dagger}\boldsymbol{\xi}_j(0). \end{align} \tag{ B.7 }$$

We want to ensure $\left\langle \Psi\left|\hat{{\cal T}}\right|j,0\right\rangle = 0$. We find $\hat{{\cal T}}|j,0\rangle$ by putting $\boldsymbol{\xi}_j(0) = \left(1\ 0\ \dots\ 0\right)^T\equiv\boldsymbol{u}_j$ and then, the condition $\left\langle \Psi\left|\hat{{\cal T}}\right|j,0\right\rangle = 0$ becomes:

$$\begin{align} \boldsymbol{\Psi}_j^{\dagger}\left(V_{\textrm{off}}e^{- i\Lambda_{\textrm{off}}t_{\textrm{off}}}V_{\textrm{off}}^{\dagger}V_{\textrm{on}}e^{- i\Lambda_{\textrm{on}}t_{\textrm{on}}}V_{\textrm{on}}^{\dagger}\right)\boldsymbol{u}_j = 0, \end{align} \tag{ B.8 }$$

where the vector $\boldsymbol{\Psi}_j = \left(\left \langle\, j,0\left|\Psi\right.\right\rangle\ \left \langle\, j - 1,1\left|\Psi\right.\right\rangle\ \dots\ \left \langle0,j\left|\Psi\right.\right\rangle\right)^T$. This can be numerically solved for $\{t_{\textrm{on}},t_{\textrm{off}}\}$.

In the special case of a two-level approximation (ignoring higher levels) and for sufficiently large detuning $g_0\ll\omega_m - \omega_q$, the above condition reduces to:

$$\begin{align} \cos\left(g_0t_{\textrm{on}}^{(0)}\sqrt{j}\right)\left \langle\, j,0\left|\Psi\right.\right\rangle + ie^{- i\Delta_0t_{\textrm{off}}^{(0)}}\sin\left(g_0t_{\textrm{on}}^{(0)}\sqrt{j}\right)\left \langle\, j - 1,1\left|\Psi\right.\right\rangle = 0. \end{align} \tag{ B.9 }$$

This gives:

$$\begin{align} g_0\sqrt{j}t_{\textrm{on}}^{(0)} & = \tan^{- 1}\left|\frac{\left \langle\, j,0\left|\Psi\right.\right\rangle}{\left \langle\, j - 1,1\left|\Psi\right.\right\rangle}\right|, \end{align} \tag{ B.10 }$$

$$\begin{align} \Delta_0t_{\textrm{off}}^{(0)} & = \textrm{arg}\left[\frac{\left \langle\, j,0\left|\Psi\right.\right\rangle}{\left \langle\, j - 1,1\left|\Psi\right.\right\rangle}\right] + \frac{\pi}{2}. \end{align} \tag{ B.11 }$$

In some special cases, equation (B.8) will have no solution. Then, we try the ansatz with a pulse of width $t_{\textrm{on}}^{(0)}$, followed by an off-time $t_{\textrm{off}}^{(0)}$ and a second pulse of an undetermined time t_on,

$$\begin{align} \Delta(t) = \Delta_0\left[u\left(t - t_{\textrm{on}} - t_{\textrm{tot}}^{(0)}\right) - u\left(t - t_{\textrm{tot}}^{(0)}\right) + u\left(t - t_{\textrm{on}}^{(0)}\right) - u\left(t\right)\right], \end{align} \tag{ B.12 }$$

where $t_{\textrm{tot}}^{(0)} = t_{\textrm{on}}^{(0)} + t_{\textrm{off}}^{(0)}$. After another off-period of time t_off, the optimization condition becomes (following the same analysis as before):

$$\begin{align} \boldsymbol{\Psi}_j^{\dagger}\left(V_{\textrm{off}}e^{- i\Lambda_{\textrm{off}}t_{\textrm{off}}}V_{\textrm{off}}^{\dagger}V_{\textrm{on}}e^{- i\Lambda_{\textrm{on}}t_{\textrm{on}}}V_{\textrm{on}}^{\dagger}\right)\left(V_{\textrm{off}}e^{- i\Lambda_{\textrm{off}}t_{\textrm{off}}^{(0)}}V_{\textrm{off}}^{\dagger}V_{\textrm{on}}e^{- i\Lambda_{\textrm{on}}t_{\textrm{on}}^{(0)}}V_{\textrm{on}}^{\dagger}\right)\boldsymbol{u}_j = 0, \end{align} \tag{ B.13 }$$

which can be solved for $\{t_{\textrm{on}},t_{\textrm{off}}\}$. This process can be repeated indefinitely, however, in practice, we did not observe the need to go beyond a second iteration.

We now describe how the transmon gates can be implemented in our scheme. In what follows, $|m,q\rangle$ denotes a state with m magnons and transmon at qth level. Given an initial state $|\Psi\rangle$ with a maximum of j magnons and $(j + 1)$ excitations, i.e. $\left \langle m,q\left|\Psi\right.\right\rangle = 0$ for $m + q\gt j + 1$ and $\left \langle\, j + 1,0\left|\Psi\right.\right\rangle = 0$, and a tunable Hamiltonian:

$$\begin{align} \frac{\skew{2.4}\hat{H}(t)}{\hbar} = \omega_m\skew{2.4}\hat{m}^{\dagger}\skew{2.4}\hat{m} + \omega_q\skew{2.4}\hat{q}^{\dagger}\skew{2.4}\hat{q} - \sum_{n = 2}^{\infty}\alpha_n|n\rangle\langle n| + \skew{2.4}\tilde{\varepsilon}(t)\skew{2.4}\hat{q}^{\dagger} + \skew{2.4}\tilde{\varepsilon}^*(t)\skew{2.4}\hat{q}, \end{align} \tag{ C.1 }$$

we want to find $\skew{2.4}\tilde{\varepsilon}(t)$ s.t. $\left \langle\, j,1\left|\hat{{\cal Q}}^{\dagger}\right|\Psi\right\rangle = 0$, where the time-evolution operator:

$$\begin{align} \hat{{\cal Q}} = \mathrm{Texp}\left[ - \frac{i}{\hbar}\int_0^{t_g}\skew{2.4}\hat{H}(\tau)d\tau \right], \end{align} \tag{ C.2 }$$

with Texp being the time ordered exponential operator and t_g being a pre-determined gate time.

Note that we have ignored the coupling between magnons and transmons assuming large enough detuning $\omega_m - \omega_q$. Then, the evolution is given by:

$$\begin{align} \hat{{\cal Q}}(t) = e^{- i\omega_mt\skew{2.4}\hat{m}^{\dagger}\skew{2.4}\hat{m}}\otimes\hat{{\cal U}}_g(t). \end{align} \tag{ C.3 }$$

Here,

$$\begin{align} \frac{d}{dt}\hat{{\cal U}_g} = - i\left(\omega_q\skew{2.4}\hat{q}^{\dagger}\skew{2.4}\hat{q} - \sum_{n = 2}^{\infty}\alpha_n|n\rangle\langle n| + \skew{2.4}\tilde{\varepsilon}(t)\skew{2.4}\hat{q}^{\dagger} + \skew{2.4}\tilde{\varepsilon}^*(t)\skew{2.4}\hat{q}\right)\hat{{\cal U}}_g. \end{align} \tag{ C.4 }$$

An ideal gate between $|0\rangle$ and $|1\rangle$, with angular parameters $\{\theta,\phi\}$ is:

$$\begin{align} \hat{{\cal U}}_{\textrm{ideal}} = e^{i\Phi_{\textrm{global}}}\left[\cos\theta\left(|0\rangle\langle 0| + |1\rangle\langle 1|\right) + i\sin\theta\left(e^{i\phi}|1\rangle\langle 0| + e^{- i\phi}|0\rangle\langle 1|\right) + \sum_{n = 2}^{\infty}|n\rangle\langle n|\right], \end{align} \tag{ C.5 }$$

barring an unimportant global phase $\Phi_{\textrm{global}}$. The parameters for this gate can be found by ensuring $\left \langle\, j,1\left|\hat{{\cal Q}}_{\textrm{ideal}}^{\dagger}\right|\Psi\right\rangle = 0$, where $\hat{{\cal Q}}_{\textrm{ideal}} = e^{- i\omega_mt_g\skew{2.4}\hat{m}^{\dagger}\skew{2.4}\hat{m}}\otimes\hat{{\cal U}}_{\textrm{ideal}}$ giving:

$$\begin{align} \theta_{\textrm{ideal}} = \tan^{- 1}\left|\frac{\left \langle\, j,1\left|\Psi\right.\right\rangle}{\left \langle\, j,0\left|\Psi\right.\right\rangle}\right|, \end{align} \tag{ C.6 }$$

and,

$$\begin{align} \phi_{\textrm{ideal}} = - \frac{\pi}{2} + \textrm{arg}\left[\frac{\left \langle\, j,1\left|\Psi\right.\right\rangle}{\left \langle\, j,0\left|\Psi\right.\right\rangle}\right]. \end{align} \tag{ C.7 }$$

To implement this, we try the ansatz $\skew{2.4}\tilde{\varepsilon}(t) = - \varepsilon(t)e^{- i\omega_qt + i\phi_{\textrm{app}}}$ with real $\{\varepsilon(t),\phi_\textrm{app}\}$ to be found. We separate out the phases by the transformation,

$$\begin{align} \hat{{\cal U}}_g(t) = e^{- i\left(\omega_q\skew{2.4}\hat{q}^{\dagger}\skew{2.4}\hat{q} - \sum_{n = 2}^{\infty}\alpha_n|n\rangle\langle n|\right)t}\hat{{\cal U}}_{\textrm{ph}}\hat{{\cal U}}_{\theta}(t)\hat{{\cal U}}_{\textrm{ph}}^{\dagger}, \end{align} \tag{ C.8 }$$

where

$$\begin{align} \hat{{\cal U}}_{\textrm{ph}} = e^{i\phi_{\textrm{app}}\skew{2.4}\hat{q}^{\dagger}\skew{2.4}\hat{q}}, \end{align} \tag{ C.9 }$$

and $\hat{{\cal U}}_{\theta}$ satisfies:

$$\begin{align} \frac{d}{dt}\hat{{\cal U}}_{\theta} = i\varepsilon(t)\left(\sum_{n = 0}^{\infty}\sqrt{n + 1}\ e^{- i\left(\alpha_{n + 1} - \alpha_n\right)t}|n + 1\rangle\langle n| + h . c . \right)\hat{{\cal U}}_{\theta}, \end{align} \tag{ C.10 }$$

where by definition $\alpha_0 = \alpha_1 = 0$.

First, consider the two-level approximation, where we want to solve:

$$\begin{align} \frac{d}{dt}\hat{{\cal U}}_{\theta}^{(0)} = i\varepsilon^{(0)}(t)\left(|0\rangle\langle 1| + |1\rangle\langle 0|\right)\hat{{\cal U}}_{\theta}^{(0)}, \end{align} \tag{ C.11 }$$

where $\varepsilon^{(0)}(t)$ is the 'bare' solution to be chosen. The unitary becomes:

$$\begin{align} \hat{{\cal U}}_{\theta}^{(0)}(t) = \cos\Theta(t)\left(|0\rangle\langle 0| + |1\rangle\langle 1|\right) + i\sin\Theta(t)\left(|1\rangle\langle 0| + |0\rangle\langle 1|\right) + \sum_{n = 2}^{\infty}|n\rangle\langle n|, \end{align} \tag{ C.12 }$$

where

$$\begin{align} \Theta(t) = \int_0^t\varepsilon^{(0)}(\tau)d\tau. \end{align} \tag{ C.13 }$$

For such a case, we have:

$$\begin{align*} \hat{{\cal U}}_{g}^{(0)}(t)& = \left[\cos\Theta(t)\left(|0\rangle\langle 0|+e^{-i\omega_{q}t_{g}}|1\rangle\langle 1|\right)\right.\\ & \qquad + \left. i\sin\Theta(t)\left(e^{i\phi_{\textrm{app}}^{(0)}}e^{-i\omega_{q}t_{g}}|1\rangle\langle 0|+e^{-i\phi_{\textrm{app}}^{(0)}}|0\rangle\langle 1|\right)+\sum_{n = 2}^{\infty}e^{-i\left(\omega_{q}-\alpha_{n}\right)t}|n\rangle\langle n|\right].\end{align*}$$

This resembles the ideal gate $\hat{{\cal U}}_{\textrm{ideal}}$ except for the extra phases due to the free evolution of the states. For this two-level approximation, we can ensure $\left \langle\, j,1\left|\hat{{\cal Q}}^{(0),\dagger}\right|\skew{2.4}\tilde{\Psi}_{j + 1}\right\rangle = 0$, where $\hat{{\cal Q}}^{(0)} = e^{- i\omega_mt_g\skew{2.4}\hat{m}^{\dagger}\skew{2.4}\hat{m}}\otimes\hat{{\cal U}}_g^{(0)}(t_g)$ by choosing:

$$\begin{align} \Theta(t_g) = \theta_{\textrm{ideal}},\ \phi_{\textrm{app}}^{(0)} = \phi_{\textrm{ideal}} + \omega_qt_g . \end{align} \tag{ C.14 }$$

Notice that we can choose any $\varepsilon^{(0)}(t)$ satisfying $\int_0^{t_g}d\tau\varepsilon^{(0)}(\tau) = \theta_{\textrm{ideal}}$. For concreteness, we choose a sinusoidal:

$$\begin{align} \varepsilon^{(0)}(t) = \frac{\pi\theta_{\textrm{ideal}}}{2t_g}\sin\left(\frac{\pi t}{t_g}\right).\end{align} \tag{ C.15 }$$

C.1. Including all levels

Transition between $|1\rangle$ and $|2\rangle$ are non-resonant by the frequency α₂, so the above two-level approximation works when $t_g\gg\alpha_2$ [assuming that the anharmonicity increases with levels, i.e. $\alpha_2\lt\alpha_3\lt\alpha_4\lt\dots$]. For a general t_g , we modify the pulse $\varepsilon(t) = \varepsilon^{(0)}(t) + \varepsilon^{(1)}(t)$ and define the error as $\hat{{\cal U}}_{\theta} = \hat{{\cal U}}_{\theta}^{(0)}\hat{{\cal U}}_{\textrm{err}}$. The error satisfies:

$$\begin{align} \frac{d\hat{{\cal U}}_{\textrm{err}}}{dt} = i\delta\skew{2.4}\hat{H}\ \hat{{\cal U}}_{\textrm{err}}, \end{align} \tag{ C.16 }$$

where the error Hamiltonian is:

$$\begin{align} \delta\skew{2.4}\hat{H} = \hat{{\cal U}}_{\theta}^{(0)\dagger}\left(\varepsilon(t)\sum_{n = 0}^{\infty}\sqrt{n + 1}\ e^{- i\left(\alpha_{n + 1} - \alpha_n\right)t}|n + 1\rangle\langle n| - \varepsilon^{(0)}(t)\left(|0\rangle\langle 1| + |1\rangle\langle 0|\right) + h . c . \right)\hat{{\cal U}}_{\theta}^{(0)}. \end{align} \tag{ C.17 }$$

This can be divided as:

$$\begin{align} \delta\skew{2.4}\hat{H} = \left[\varepsilon(t)\left(\delta\hat{h}_{12} + \delta\hat{h}_{\textrm{rest}}\right) + \left(\varepsilon(t) - \varepsilon^{(0)}(t)\right)\delta\hat{h}_{01}\right], \end{align} \tag{ C.18 }$$

with the reduced Hamiltonians defined as:

$$\begin{align*} \delta\hat{h}_{01} & = |0\rangle\langle 1|+|1\rangle\langle 0|\\ \delta\hat{h}_{12} & = \sqrt{2}\ e^{-i\alpha_{2}t}\left[\cos\Theta(t)|2\rangle\langle 1|+i\sin\Theta(t)|2\rangle\langle 0|\right]+h.c.\\ \delta\hat{h}_{\textrm{rest}} & = \sum_{n = 2}^{\infty}\sqrt{n+1}\ e^{-i\left(\alpha_{n+1}-\alpha_{n}\right)t}|n+1\rangle\langle n|+h.c. \end{align*}$$

We perform a Magnus expansion to find [54, 55]:

$$\begin{align} \hat{{\cal U}}_{\textrm{err}}(t) = \exp\left[i\int_0^td\tau\delta\skew{2.4}\hat{H}(\tau) - \frac{1}{2}\int_0^td\tau_1\int_0^{\tau_1}d\tau_2\left[\delta\skew{2.4}\hat{H}(\tau_1),\delta\skew{2.4}\hat{H}(\tau_2)\right] + \dots\right], \end{align} \tag{ C.19 }$$

so, to ensure $\hat{{\cal U}}_{\textrm{err}}(t_g) \approx \hat{I}$ up to first order, we require $\int_0^{t_g}d\tau\delta\skew{2.4}\hat{H}(\tau) = \hat{0}$ which translates to:

$$\begin{align} \int_0^{t_g}\varepsilon(\tau)d\tau & = \theta_{\textrm{ideal}}, \end{align} \tag{ C.20 }$$

$$\begin{align} \int_0^{t_g}\varepsilon(\tau)e^{i\alpha_2\tau\pm i\Theta(t)}d\tau & = 0, \end{align} \tag{ C.21 }$$

$$\begin{align} \int_0^{t_g}\varepsilon(\tau)e^{i\left(\alpha_3 - \alpha_2\right)\tau}d\tau & = 0. \end{align} \tag{ C.22 }$$

We used $\int_0^{t_g}\varepsilon^{(0)}d\tau = \theta_{\textrm{ideal}}$. The last condition ensures that there are no $|2\rangle\leftrightarrow|3\rangle$ transitions, and we ignored all higher transitions as they are approximately unoccupied throughout the protocol.

C.2. Satisfying the conditions

The conditions allow for a lot of freedom in choosing $\varepsilon(t)$. Here, we take an ansatz of a sum of sinusoids, ensuring that the pulse starts and ends at 0,

$$\begin{align} \varepsilon(t) = \sum_{\nu}\varepsilon_{\nu}\sin\left(\frac{\nu\pi t}{t_g}\right), \end{align} \tag{ C.23 }$$

where the set of indices ν can be freely chosen. Equation (C.20) reads:

$$\begin{align} \sum_{\nu\:\textrm{odd}}\frac{2\varepsilon_{\nu}t_g}{\pi\nu} = \theta_{\textrm{ideal}}. \end{align} \tag{ C.24 }$$

Equation (C.22) is:

$$\begin{align} \sum_{\nu}\frac{\pi\nu\varepsilon_{\mu}t_g}{\left(\alpha_3 - \alpha_2\right)^2t_g^2 - \pi^2\nu^2}\left(1 - \left( - 1\right)^{\nu}e^{i\left(\alpha_3 - \alpha_2\right)t_g}\right) = 0. \end{align} \tag{ C.25 }$$

Equation (C.21) is:

$$\begin{align} e^{\pm i\theta/2}\sum_{\nu}\frac{\varepsilon_{\mu}t_g}{\pi}\int_0^{\pi}du\sin\left(\nu u\right)\exp\left[i\frac{\alpha_2t_g}{\pi}u\mp\frac{i\theta}{2}\cos u\right] = 0. \end{align} \tag{ C.26 }$$

To simplify the integration, we define a real integral:

$$\begin{align} G_{\zeta}(\xi) = \int_0^{\pi/2}du\cos\left[\zeta u - \xi\sin u\right] , \end{align} \tag{ C.27 }$$

giving

$$\begin{align} \int_0^{\pi}du\exp\left[i\zeta u + i\xi\cos u\right] = 2i^{\zeta}G_{\zeta}(\xi). \end{align} \tag{ C.28 }$$

For $\zeta = \alpha_2t_g/\pi$, this gives:

$$\begin{align} 0 & = \sum_{\nu\:\textrm{odd}}i^{\nu - 1}\frac{\varepsilon_{\nu}t_g}{\pi}\left[G_{\zeta + \nu}\left(\frac{\pm\theta}{2}\right) + G_{\zeta - \nu}\left(\frac{\pm\theta}{2}\right)\right] \end{align} \tag{ C.29 }$$

$$\begin{align} & = \sum_{\nu\:\textrm{even}}i^{\nu}\frac{\varepsilon_{\nu}t_g}{\pi}\left[G_{\zeta + \nu}\left(\frac{\pm\theta}{2}\right) - G_{\zeta - \nu}\left(\frac{\pm\theta}{2}\right)\right]. \end{align} \tag{ C.30 }$$

If we choose the gate time s.t. $\left(\alpha_3 - \alpha_2\right)t_g = 2N\pi$ for some integer N, we find that the odd and even frequency components can be separated, and in particular, we need only the odd components to satisfy all the above conditions.

For simulation purposes, we choose the specific case of $\alpha_n = An(n - 1)/2$, the bare transmon, and choose $At_g = N\pi$ for some natural number N giving (for some arbitrary set of indices H),

$$\begin{align} \varepsilon(t) & = \sum_{\mu\in H}\varepsilon_{\mu}\sin\left(\frac{\left(2\mu+1\right)\pi t}{t_{g}}\right), \end{align} \tag{ C.31 }$$

$$\begin{align} \frac{\pi\theta_{\textrm{ideal}}}{2t_{g}} & = \sum_{\mu\in H}\frac{\varepsilon_{\mu}}{2\mu+1}, \end{align} \tag{ C.32 }$$

$$\begin{align} 0 & = \sum_{\mu\in H}\frac{\left(2\mu+1\right)\varepsilon_{\mu}}{4N^{\,2}-\left(2\mu+1\right)^{2}}, \end{align} \tag{ C.33 }$$

$$\begin{align} 0 & = \sum_{\mu\in H}\left(-1\right)^{\mu}\varepsilon_{\mu}\frac{G_{N+2\mu+1}(\pm\theta/2)+G_{N-2\mu-1}\left(\pm\theta/2\right)}{2}. \end{align} \tag{ C.34 }$$

This is a set of linear equations which can be solved to give ε_µ . To find an expression for $G_n(\theta)$, we use the fact that it satisfies the same recursion relation as that of Bessel functions:

$$\begin{align} 2G_n^{^{\prime}}(\xi) = G_{n - 1}(\xi) - G_{n + 1}(\xi), \end{align} \tag{ C.35 }$$

with the initial conditions,

$$\begin{align} G_0(\xi) - \frac{\pi}{2}J_0(\xi) = 0,\ G_1(\xi) - \frac{\pi}{2}J_1(\xi) = \frac{\sin\xi}{\xi}. \end{align} \tag{ C.36 }$$

Then, an expression for $G_n(\xi) - \pi J_n(\xi)/2$ can be found recursively.

C.3. Phase

Given $\varepsilon(t)$, we find φ_app to minimize $\left|\left \langle\, j,1\left|\hat{{\cal Q}}^{\dagger}\right|\Psi\right\rangle\right|$. We have $\hat{{\cal Q}}(t) = e^{- i\omega_mt\skew{2.4}\hat{m}^{\dagger}\skew{2.4}\hat{m}}\otimes\hat{{\cal U}}_g(t)$ where $\hat{{\cal U}}_g(t)$ can be decomposed as in equation (C.8)]. Then,

$$\begin{align} \hat{{\cal U}_g}(t_g)|1\rangle = e^{- i\phi_{\textrm{app}}}e^{- i\left(\omega_q\skew{2.4}\hat{q}^{\dagger}\skew{2.4}\hat{q} - \sum_{n = 2}^{\infty}\alpha_n|n\rangle\langle n|\right)t_g}\hat{{\cal U}}_{\textrm{ph}}\hat{{\cal U}}_{\theta}(t_g)|1\rangle. \end{align} \tag{ C.37 }$$

As the leakage into higher levels is suppressed via an appropriate choice of $\varepsilon(t)$, we expect the 2x2 sub-block of $\hat{{\cal U}}_g(t_g)$ for the subspace spanned by $\{|0\rangle,|1\rangle\}$ to be approximately unitary. For concreteness, we find the unitary matrix closest to this sub-block, say $\hat{\cal U}_{g,01}$, in 2-norm via $\hat{\cal U}_{g,01}^{\textrm{unit}} = VW^{\dagger}$ using the singular value decomposition $\hat{\cal U}_{g,01} = V\Lambda W^{\dagger}$. This gives:

$$\begin{align} \hat{\cal U}_{g,01}^{\textrm{unit}}|1\rangle = e^{i\Phi_G}\left(i\sin\tilde{\theta}e^{i\delta\phi}|0\rangle + \cos\tilde{\theta}|1\rangle\right), \end{align} \tag{ C.38 }$$

for some constants $\{\Phi_G,\tilde{\theta},\delta\phi\}$. Then,

$$\begin{align} \left\langle \Psi\left|\hat{\cal U}_g(t_g)\right|j,1\right\rangle \approx e^{- ij\omega_mt - i\phi_{\textrm{app}} + i\Phi_G}\left(i\sin\tilde{\theta}e^{i\delta\phi}\left\langle \Psi\left|j,0\right.\right\rangle + \cos\tilde{\theta}e^{i\phi_{\textrm{app}} - i\omega_qt_g}\left\langle \Psi\left|j,1\right.\right\rangle\right), \end{align} \tag{ C.39 }$$

giving the amplitude:

$$\begin{align} \left|\left\langle\Psi\left|\hat{\cal U}_g(t_g)\right|j,1\right\rangle\right| \approx \sqrt{\left|\left\langle\Psi\left|j,0\right.\right\rangle\right|^2 + \left|\left\langle\Psi\left|j,1\right.\right\rangle\right|^2}\left|\cos\theta_{\textrm{ideal}}\sin\tilde{\theta} - \cos\tilde{\theta}\sin\theta_{\textrm{ideal}}e^{i\phi_{\textrm{app}} - i\phi_{\textrm{app}}^{(0)} - i\delta\phi}\right|. \end{align} \tag{ C.40 }$$

This is minimized at:

$$\begin{align} \phi_{\textrm{app}} = \phi_{\textrm{app}}^{(0)} + \delta\phi. \end{align} \tag{ C.41 }$$

C.4. Gates between higher levels

Typically, there would be leakage from magnons into higher levels of transmon, which needs to be corrected via applying gates between the higher levels. Consider the case when we want to apply a gate between $|n\rangle$ and $|n + 1\rangle$ with angular parameters $\{\theta,\phi\}$. We try for a pulse shape:

$$\begin{align} \skew{2.4}\tilde{\varepsilon}(t) = - \frac{\varepsilon(t)}{\sqrt{n}}e^{i\phi_{\textrm{app}}}e^{- i\left(n\omega_m - \alpha_n\right)t}. \end{align} \tag{ C.42 }$$

Again, we have the 'bare' pulse $\varepsilon^{(0)}$ which we take as half-sines. To remove the errors at first order in the Magnus expansion, we find that $\varepsilon(t)$ here satisfies the same set of conditions equations (C.20)–(C.22). Thus, the same calculations as above apply.