Introduction

Quantum computers take advantage of quantum coherence and entanglement that do not exist in the classical world to process information. A range of important problems such as factorizing large integers1,2 and simulating quantum many-body dynamics3,4 are promised to be solved by quantum computers exponentially faster than classical computers5. The basic unit of a quantum computer to encode information is quantum bits (qubits) made from discrete-variable (DV) systems such as spins and quantized levels. However, due to the ubiquitous noises in the environment, quantum states are fragile and the encoded information in qubits can be inevitably deteriorated and lost6. For practical quantum computing, it is crucial to implement quantum error correction (QEC) to protect against unwanted and uncontrolled errors, which remains one of the most challenging and urgent goals for building a fault-tolerant quantum computer7,8,9.

Representative QEC codes with DV qubits towards practical large-scale quantum computation include stabilizer codes10 and surface codes11,12,13, which encode a logical qubit in a subspace of multiple physical qubits so that different error events lead to distinguishable detectable syndromes, thus facilitating the recovery of corrupted quantum states. However, such DV-based QEC schemes typically consume substantial physical resources, as encoding a single logical qubit often demands a large number of physical qubits. Furthermore, logical gate operations in these schemes are complicated, as multiple physical qubits have to be manipulated simultaneously. Scaling up physical qubits to build a fault-tolerant quantum computer remains an exceptionally challenging task due to these significant demands on physical qubits and operation complexity14,15.

The core idea of QEC lies in protecting the encoded quantum states by introducing a redundant Hilbert space for error detection and correction. A single harmonic oscillator already provides an infinitely large Hilbert space that can be partitioned into a logical subspace and error subspace16,17,18,19,20. Continuous-variable (CV) systems, such as harmonic oscillators or other bosonic systems, naturally provide an infinite-dimensional Hilbert space that can be partitioned into logical and error subspaces, making them an attractive alternative for fault-tolerant universal quantum computation16,17,18,19,20. Unlike DV systems, CV systems leverage the continuous nature of their state space, which inherently offers redundancy without the need of a large number of physical qubits21,22.

It has been proved that a single nonlinear term in addition to linear Gaussian operations is possible to realize universal single-mode CV quantum computation23, e.g., a common universal gate set includes the cubic phase gate, displacement gate, squeezing gate, and Fourier gate (phase rotation). In superconducting circuit experiments, the cubic phase gate has been proposed and realized with an on-chip planar resonator terminated by a superconducting nonlinear asymmetric inductive element (SNAIL)24,25. Another popular universal set for single-mode bosonic quantum computation is composed of a displacement gate and a selective number-dependent arbitrary phase (SNAP) gate, which imparts an arbitrary phase to every Fock number state26 using an off-resonantly coupled qubit27. The universal CV quantum computation for multiple bosonic modes can be implemented by adding a simple two-mode linear operation28,29.

However, the above universal gate operations for CV quantum computation are not guaranteed to be fault-tolerant due to the intrinsic continuous error channels of CV mode quadratures. To achieve practical fault-tolerant quantum computation with CV systems, it is necessary to embed a finite-dimensional code space into the infinite-dimensional CV Hilbert space, similar to QEC using DV systems30. Since the most common and practically relevant Gaussian CV errors generally cannot be suppressed using solely Gaussian states and Gaussian operations31, non-Gaussian states are typically required to implement fault-tolerant encodings32,33,34. However, not all non-Gaussian states are suitable for fault-tolerant quantum computing. Depending on the error set relevant in experiments32,33,34, the constructed discretized non-Gaussian states have to satisfy the so-called Knill-Laflamme condition5,35, and are known as bosonic QEC codes, or simply bosonic codes36,37.

According to the symmetries in phase space, bosonic codes can be broadly classified into translational codes, rotational codes, and vortex codes38,39,40. A prominent example of translational codes is the Gottesman-Kitaev-Preskill (GKP) code18,41, while rotational bosonic codes include cat codes19,42 and binomial codes20. There are also multimode bosonic codes that are superposition states with the same excitation number by combining multiple modes16,43. Compared to conventional multiple-qubits QEC codes, the bosonic QEC modes are hardware-efficient as the logical qubit can be built from a single quantum system with limited error channels like photon loss and pure dephasing. Various bosonic codes have been realized in experiments44,45,46,47,48,49 and have won the break-even point with cat codes9 and binomial codes50.

As both information encoding and error correction continuously consume bosonic code states, it is crucial to prepare high-quality non-Gaussian CV states on demand38, which requires some form of nonlinearity in CV modes. In optical or mechanical systems, due to the weak accessible nonlinearities (relative to the intrinsic losses), the common strategy is to combine Gaussian operations with non-Gaussian measurements such as photon number subtraction performed by single-photon detectors51,52,53. In contrast, superconducting circuit architectures54 leverage Josephson junctions (JJs) to provide strong and controllable nonlinearities, enabling the engineering of non-Gaussian bosonic states. The JJ-based transmon acts as an artificial atomic system and the lowest two levels are usually harnessed to be a physical qubit, which is utilized to implement arbitrary nonlinear phase gates by repeating the noncommuting Rabi gates55,56 or the SNAP gates using a weak drive off-resonant drive27.

For engineering bosonic code states with gate operations, how to decompose the desired arbitrary unitary operation is another complex and challenging problem. For example, the SNAP gates inherently rely on an incremental approach to implement phase rotations, and strong gates with large rotations would require many repetitions limiting resource efficiency56. In practice, the gate sequence and gate parameters to realize a target operation rely on numerical optimization techniques49,57. Recent works have proposed to engineer code states with passive control based on Hamiltonian engineering58,59,60,61 or reservoir engineering62, which can be leveraged to facilitate fault-tolerant operations58,60,61,63.

In this work, we propose a universal gate set that contains only one type of gate operation, which we refer to as the quantum lattice gate, cf. Eq. (4). This gate set is in contrast to the universal gate set with cubic phase gate, which comprises four distinct gate types with three gate parameters in total, cf. Eq. (1), and the SNAP gate set, which includes two types of gates but requires an infinite number of gate parameters (depending on the truncated number of Fock states), cf. Eq. (3). Moreover, we provide an analytical deterministic framework to prepare and transform bosonic code states with quantum lattice gates based on Floquet Hamiltonian engineering and decompose the state preparation/transformation process with a sequence of quantum lattice gates. We apply our method to specific examples, including the preparation of single binomial code states, embedding binomial code space, and transforming binomial codes to cat codes. We also demonstrate autonomous error correction for cat states in the presence of single-photon losses.

Methods

In this main section, we will first introduce a universal set of quantum lattice gates and the bosonic codes that satisfy the Knill-Laflamme condition for this work. Then, we classify the code state engineering into three basic processes and provide an analytical framework to implement these processes based on Floquet Hamiltonian engineering. We provide an analytical framework to decompose a given code state engineering process into a sequence of quantum lattice gates. Last, we apply our method to specific examples, including preparing single binomial code states, embedding binomial code space, transforming binomial codes to cat codes, and demonstrating autonomous error correction for cat states in the presence of single-photon losses.

Universal gate sets for bosonic modes

To achieve universal CV quantum computation, it is essential to implement any unitary transformation between two bosonic states by using a set of universal quantum gates generated by Hamiltonians that are arbitrary polynomials of the CV mode quadratures \(\hat{x}\) and \(\hat{p}\)23,28.

For example, the universal gate set for a single CV mode can be chosen as follows23,24,29

$$\left\{{e}^{i\frac{\pi }{2}{\hat{a}}^{{{\dagger}} }\hat{a}},{e}^{ir\hat{p}},{e}^{is{\hat{x}}^{2}},{e}^{i\gamma {\hat{x}}^{3}}\right\},$$
(1)

where \(\hat{a}\equiv (\hat{x}+i\hat{p})/\sqrt{2\lambda }\) is the ladder operator with dimensionless Planck constant λ given via \([\hat{x},\hat{p}]=i\lambda\). This universal gate set includes three Gaussian gates, i.e., the phase rotation gate \({e}^{i\frac{\pi }{2}{\hat{a}}^{{{\dagger}} }\hat{a}}\), the displacement gate \({e}^{ir\hat{p}}\), and the squeezing gate \({e}^{is{\hat{x}}^{2}}\), allowing to perform any linear transformations between the CV modes generated by an arbitrary quadratic Hamiltonian. The cubic phase gate \({e}^{i\gamma {\hat{x}}^{3}}\) is a non-Gaussian gate operation that realizes arbitrary nonlinear transformations between the CV modes24,29,55. By using the following commutator-based Lloyd-Braunstein decomposition

$${e}^{-[\hat{A},\hat{B}]\delta t}={e}^{i\hat{A}\sqrt{\delta t}}{e}^{i\hat{B}\sqrt{\delta t}}{e}^{-i\hat{A}\sqrt{\delta t}}{e}^{-i\hat{B}\sqrt{\delta t}}+O(\delta {t}^{3}),$$
(2)

any desired Hamiltonian term as an arbitrary polynomial of the mode quadratures can be generated23,24,64. Another commonly used single-mode universal gate set contains the following two types of gates26

$$\left\{{e}^{ir\hat{p}},\,S[\overrightarrow{\theta }]\equiv {\sum }_{n=0}^{\infty }{e}^{i{\theta }_{n}}| n\left.\right\rangle \left\langle \right.n| \right\},$$
(3)

where \(S[\overrightarrow{\theta }]\) is the SNAP gate that endows arbitrary phase \(\overrightarrow{\theta }={\{{\theta }_{n}\}}_{n = 0}^{\infty }\) to the Fock basis \({\{| n\left.\right\rangle \}}_{n = 0}^{\infty }\).

For multiple CV modes, combining single-mode universal gate operations with simple two-mode linear operations, such as the two-mode CSUM gate (\({e}^{-i{\hat{x}}_{1}{\hat{p}}_{2}}\)) or CZ gate (\({e}^{i{\hat{x}}_{1}{\hat{x}}_{2}}\)), is sufficient to achieve universal quantum computation23,29. This is in stark contrast to performing universal quantum computation with DV qubits. In fact, the set of single-mode Gaussian gates, together with either the CSUM gate or the CZ gate, constitutes the Clifford set that enables the implementation of Clifford quantum computations, which can be efficiently simulated on a classical computer65,66. Non-Gaussian gates for single CV modes, such as the cubic phase gate and the SNAP gate, play the role of non-Clifford gates that are essential for achieving universal gates with quantum speed-up or quantum advantage.

Quantum lattice gate

The universal gate set given by Eq. (1) is comprised of four gate elements with three gate parameters. The universal set given by Eq. (3) includes one displacement gate and one SNAP gate with the number of gate parameters determined by the truncated photon number increasing with the size of the target bosonic state. Here, we propose a universal gate set that only contains one type of elementary gate, i.e.,

$$\left\{{{{\rm{PSL}}}}(\zeta ,\sigma ;\gamma ,\delta )\equiv {e}^{i\gamma \cos (\zeta \hat{x}+\sigma \hat{p}+\delta )}\right\}.$$
(4)

We refer to this gate as the phase-space lattice (PSL) gate, as the generator of such a gate operation is a consine-lattice function in phase space of the CV mode. We further define an X-space lattice (XSL) gate by

$${{{\rm{XSL}}}}(\rho ,\gamma ,\delta )\equiv {e}^{i\gamma \cos (\rho \hat{x}+\delta )},$$
(5)

which is a special case of the PSL(ζσγδ) gate with gate parameter σ = 0. In view of this, both the PSL and the XSL gates can be referred to as quantum lattice gates. With the phase rotation gate \({{{\rm{R}}}}(\theta )\equiv {e}^{-i\theta {\hat{a}}^{{{\dagger}} }\hat{a}}\), we can decompose the PSL gate as

$${{{\rm{PSL}}}}(\zeta ,\sigma ;\gamma ,\delta )={{{\rm{R}}}}(-\theta ){{{\rm{XSL}}}}(\rho ,\gamma ,\delta ){{{\rm{R}}}}(\theta ),$$
(6)

where the gate parameters ρ and θ are determined by \(\zeta =\rho \cos \theta\) and \(\sigma =\rho \sin \theta\). We illustrate this gate decomposition with quantum circuit representation in Fig. 1a.

Fig. 1: Quantum lattice gate.
figure 1

a Quantum circuit of a phase-space lattice (PSL) gate decomposed into one X-space lattice (XSL) gate and two-phase rotation gates R(±θ), cf. Eq. (6). The upper and lower lines of the circuit represent the bosonic state stored in the cavity and the external driving field exerted on the cavity, respectively. b Schematic illustration of implementing the PSL gate with a harmonic oscillator with bare frequency ω0 kicked by a cosine-lattice potential (left). The kick exerted at time moment t = θ/ω0 realizes the XSL gate \({e}^{i\gamma \cos (\rho x+\delta )}\); the free evolutions of oscillator before and after the kick realize the phase rotation gates R(θ) and R(−θ), respectively (right).

To elucidate how the PSL gate can be implemented, we consider a harmonic oscillator with a bare frequency ω0, as shown in Fig. 1b. The phase rotation gate R(θ) can be simply generated by the free-time evolution with θ = ω0t. Subsequently, the oscillator is kicked by a lattice potential \(V(x)\propto \cos (\rho \hat{x}+\delta )\), which results in the XSL gate operation described by Eq. (5). Finally, a following free-time evolution over one harmonic period (θ = 2π − ω0t) completes the process by implementing a phase rotation gate R(−θ). At the end of this paper, we will discuss further the experimental implementation of quantum lattice gates introduced above with JJ-based superconducting circuit architectures (cf. Fig. 10). For a typical superconducting microwave cavity operating at GHz frequencies, our proposed quantum lattice gate operation can be implemented in less than one nanosecond (<ns). In contrast, the typical cubic phase gate requires tens of nanoseconds25 while the typical SNAP gate costs several microseconds (>μs)27,49 due to the weakness of the dispersive interaction.

The nonlinearity of JJ-based superconducting circuits is a valuable resource for CV quantum information processing and universal quantum computing67. To generate a nonlinear cubic phase gate, three-wave mixing is usually employed to pick up the third-order nonlinearity and eliminate the higher-order contributions25. For the SNAP gates, only the two lowest levels of the transmon are used to interact with the cavity, and all other levels are neglected27. Both the cubic phase and the SNAP gates ignore the higher-order nonlinear terms of the JJ potential27,68, and thus do not utilize the full nonlinearity as a quantum resource. Even worse, the residual nonlinearity causes errors as a destructive channel and could become a dominant source of error for engineering large states25. In contrast, our proposed lattice gates directly utilize the entire nonlinear cosine potential offered by a JJ by coherently harnessing all the high-order nonlinear terms.

Bosonic codes

While the universal gates introduced above allow for arbitrary transformations over CV modes, there is no guarantee that errors can be corrected during gate operations. Due to the intrinsic continuous error channels, e.g., small diffusion along the two quadratures or weak shifts in phase rotations, precise control of CV states is typically more difficult than that of DV states. In fact, the no-go theorem claims that the most common and practically relevant Gaussian CV errors cannot be suppressed solely with Gaussian states and Gaussian operations31. It is thus very important to correct errors for practical quantum computations with CV states.

Similar to QEC using DV systems, quantum error correction requires some form of discretization over CVs. Bosonic codes provide a solution to this problem by embedding a finite-dimensional code space into the infinite-dimensional Hilbert space of a CV system30. Compared to QEC codes based on DV systems, bosonic codes are more hardware-efficient, as they leverage the infinite-dimensional Hilbert space of a single CV mode to achieve the necessary redundancy36. Bosonic codes also circumvent the no-go theorem of Gaussian CVs, enabling one to detect and correct small errors31. Unlike CV universal quantum computation without error correction, which is similar to noisy intermediate-scale quantum (NISQ) simulators in the context of qubit-based quantum computation6, bosonic codes provide the possibility of fault-tolerant digital quantum computation69.

The essence of bosonic codes lies in the encoding and protection of quantum information within a code space spanned by discretized non-Gaussian states (NGSs), of which two are the code states \(| \bar{0}\left.\right\rangle\) and \(| \bar{1}\left.\right\rangle\). However, not all NGSs are suitable for bosonic code states for quantum error correction. Precisely speaking, the chosen bosonic code states have to satisfy the so-called Knill-Laflamme condition5,35

$$\left\{\begin{array}{l}\langle \bar{0}| {\hat{\varepsilon }}_{i}^{{{\dagger}} }{\hat{\varepsilon }}_{j}| \bar{0}\rangle =\langle \bar{1}| {\hat{\varepsilon }}_{i}^{{{\dagger}} }{\hat{\varepsilon }}_{j}| \bar{1}\rangle , \hfill \quad \\ \langle \bar{0}| {\hat{\varepsilon }}_{i}^{{{\dagger}} }{\hat{\varepsilon }}_{j}| \bar{1}\rangle =\langle \bar{1}| {\hat{\varepsilon }}_{i}^{{{\dagger}} }{\hat{\varepsilon }}_{j}| \bar{0}\rangle =0\quad \end{array}\right.$$
(7)

with the error set \(\varepsilon =\{\hat{I},{\hat{\varepsilon }}_{1},{\hat{\varepsilon }}_{2},{\hat{\varepsilon }}_{3},\cdots \,\}\). For cavity modes, the usual dominant error channels are the single-photon loss \({\hat{\varepsilon }}_{1}=\hat{a}\) and dephasing \({\hat{\varepsilon }}_{2}={\hat{a}}^{{{\dagger}} }\hat{a}\). The Knill-Laflamme condition guarantees that any encoded quantum state (\(| \psi \left.\right\rangle ={c}_{0}| \bar{0}\left.\right\rangle +{c}_{1}| \bar{1}\left.\right\rangle\)) corrupted by error channels (\(| \psi \left.\right\rangle \to {c}_{0}{\hat{\varepsilon }}_{i}| \bar{0}\left.\right\rangle +{c}_{1}{\hat{\varepsilon }}_{j}| \bar{1}\left.\right\rangle\)) is only transformed into another orthogonal basis without losing the quantum information encoded in the superposition coefficients.

For example, the four-legged cat code states42,44,70

$$\left\{\begin{array}{l}| {\bar{0}}_{c}\left.\right\rangle \equiv \frac{1}{\sqrt{{{{{\mathcal{N}}}}}_{0}}}\left(| \alpha \left.\right\rangle +| -\alpha \left.\right\rangle +| i\alpha \left.\right\rangle +| -i\alpha \left.\right\rangle \right)\quad \\ | {\bar{1}}_{c}\left.\right\rangle \equiv \frac{1}{\sqrt{{{{{\mathcal{N}}}}}_{2}}}\left(| \alpha \left.\right\rangle +| -\alpha \left.\right\rangle -| i\alpha \left.\right\rangle -| -i\alpha \left.\right\rangle \right),\quad \end{array}\right.$$
(8)

were proposed to correct single-photon loss characterized by the error set \(\varepsilon =\{\hat{I},\hat{a}\}\), where \({{{{\mathcal{N}}}}}_{m}=8{e}^{-{\alpha }^{2}}[\cosh {\alpha }^{2}+{(-1)}^{\frac{m}{2}}\cos {\alpha }^{2}]\) for m = 0, 2 is the corresponding normalization factor. However, the four-legged cat code states satisfy the Knill-Laflamme condition only when α is chosen at specific “sweet spots” determined by refs. 70,71

$$\tan {\alpha }^{2}=-\tanh {\alpha }^{2}.$$
(9)

Otherwise, the Knill-Laflamme condition for such cat codes is only approximately satisfied in the limit of large coherent states, i.e., α 1.

In order to further correct pure dephasing errors described by the error set \(\varepsilon =\{\hat{I},\hat{a},\hat{n}\}\), where \(\hat{n}\equiv {\hat{a}}^{{{\dagger}} }\hat{a}\) is the photon number operator, one can encode the quantum information into the so-called binomial codes20

$$\left\{\begin{array}{l}| {\bar{0}}_{b}\left.\right\rangle \equiv \frac{1}{2}\left(| 0\left.\right\rangle +\sqrt{3}| 4\left.\right\rangle \right)\quad \\ | {\bar{1}}_{b}\left.\right\rangle \equiv \frac{1}{2}\left(\sqrt{3}| 2\left.\right\rangle +| 6\left.\right\rangle \right),\quad \end{array}\right.$$
(10)

where \(| n\left.\right\rangle\) is the Fock state with n photons. One can check that the above binomial code states exactly satisfy the Knill-Laflamme condition Eq. (7) with

$$\left\{\begin{array}{l}\langle {\bar{0}}_{b}| \hat{n}| {\bar{0}}_{b}\rangle =\langle {\bar{1}}_{b}| \hat{n}| {\bar{1}}_{b}\rangle =3\hfill \\ \langle {\bar{0}}_{b}| {\hat{n}}^{2}| {\bar{0}}_{b}\rangle =\langle {\bar{1}}_{b}| {\hat{n}}^{2}| {\bar{1}}_{b}\rangle =12.\quad \end{array}\right.$$
(11)

To correct errors that are polynomials up to a specific degree in bosonic creation and annihilation operators, one has to introduce more complicated binomial code states that are formed by a finite superposition of Fock states weighted with appropriate binomial coefficients20.

In Fig. 2, we visualize the cat and binomial code states by plotting their Wigner functions (see the detailed expressions in Supplementary Note 1. A). One can find that all the Wigner functions are invariant under a \(\frac{\pi }{2}\) rotation in phase space. In fact, both the cat and binomial code states are the eigenstates of the parity operator

$$\hat{{{{\mathcal{P}}}}}=\exp (i\pi \hat{n})$$
(12)

with eigenvalue “+1” (stabilizer). A single-photon loss will change their parity to “−1”, which can be detected by nondemolition parity measurements72,73,74 and then can be corrected accordingly.

Fig. 2: Wigner functions of bosonic code states W(xp).
figure 2

a, b Two logical cat code states \(| {\bar{0}}_{c}\left.\right\rangle\) and \(| {\bar{1}}_{c}\left.\right\rangle\) defined in Eq. (8) for the first “sweet spot” (α2 ≈ 2.34) determined by Eq. (9). c, d Two logical binomial code states \(| {\bar{0}}_{b}\left.\right\rangle\) and \(| {\bar{1}}_{b}\left.\right\rangle\) defined in Eq. (10).

Code states engineering

The ability to engineer bosonic code states is one of the key ingredients for fault-tolerant bosonic quantum computation. In this section, we study three basic state engineering processes, i.e., code state preparation, code space embedding, and code space transformation, as illustrated in Fig. 3. In most bosonic quantum computing schemes with nonlinear phase gates or SNAP gates, the gate sequences and parameters are determined by various numerical optimization methods57,75,76,77. Here, in contrast, we provide an analytical and deterministic framework for engineering bosonic code states based on Floquet Hamiltonian engineering and the adiabatic ramp protocol61.

Fig. 3: Code states engineering.
figure 3

a Code state preparation (CStP) and Code space embedding (CSpE) via an adiabatic ramp (AR) from the bare cavity Hamiltonian \({\hat{H}}_{0}\) to the target Hamiltonian \({\hat{H}}_{T}\) constructed from Eq. (13) and Eq. (14) respectively. b Code space transformation (CSpT) from the binomial Hamiltonian \({\hat{H}}_{{{{\rm{bin}}}}}\) to the cat Hamiltonian \({\hat{H}}_{{{{\rm{cat}}}}}\) via a transit Hamiltonian \({\hat{H}}_{t}\), cf. Eqs. (16)–(18).

First, we introduce the code state preparation. We begin by discussing the process of code state preparation, which typically starts from the ground state of a bare cavity (harmonic oscillator), cf. the middle part in Fig. 3a, with the bare Hamiltonian \({\hat{H}}_{0}=\frac{1}{2}{\omega }_{0}({\hat{x}}^{2}+{\hat{p}}^{2}).\) Our first objective is to design a target Hamiltonian HT whose “ground state” or eigenstate is the desired target state \(| {\psi }_{T}\left.\right\rangle\). For example, \({\hat{H}}_{T}=\beta (\cos \hat{x}+\sin \hat{p})\) could be used for the GKP state18,61, and \({\hat{H}}_{T}=\beta ({\hat{a}}^{{{\dagger}} q}-1)({\hat{a}}^{q}-1)\) for the q-leg cat state58,63, where β represents the amplitude of the designed potential. However, it is generally not always possible to find a compact Hamiltonian whose eigenstate is an arbitrary target state, such as the binomial code states20. To address this difficulty, we turn to construct the target Hamiltonian directly from the target state by

$${\hat{H}}_{T}=-\Delta | {\psi }_{T}\left.\right\rangle \left\langle \right.{\psi }_{T}| ,$$
(13)

where Δ > 0 is the energy gap of the target Hamiltonian. In this formulation of the target Hamiltonian, the target state \(| {\psi }_{T}\left.\right\rangle\) has the lowest eigenvalue  −Δ, while all other orthogonal states have zero eigenvalues. To prepare the target code state \(| {\psi }_{T}\left.\right\rangle\), we employ the adiabatic ramp protocol61 from the Hamiltonian \({\hat{H}}_{0}\) of the bare cavity to the target Hamiltonian \({\hat{H}}_{T}\). As a result, the cavity state is adiabatically tuned from the initial vacuum state to the final target state. The parameter Δ provides a gap protection for the prepared state. This scheme is illustrated by the “middle-to-left” process in Fig. 3a.

Second, we introduce the code space embedding. One can also embed a finite-dimensional code space, spanned by the two logical code words \(| \bar{0}\left.\right\rangle\) and \(| \bar{1}\left.\right\rangle\), into the infinite-dimensional Fock space of the cavity by constructing the target Hamiltonian as

$${\hat{H}}_{T}=-\Delta \left(| \bar{0}\left.\right\rangle \langle \bar{0}| +| \bar{1}\rangle \left\langle \right.\bar{1}| \right).$$
(14)

The corresponding Hilbert space forms a degenerate code space spanned by the two logical states, with a nonzero energy gap separating them from the other orthogonal eigenstates. As discussed above, we adiabatically ramp the system Hamiltonian from \({\hat{H}}_{0}\) to \({\hat{H}}_{T}\), and prepare the two bosonic code states from the Fock ground state \(| 0\left.\right\rangle\) and the first excited state \(| 1\left.\right\rangle\), i.e., \(| 0\left.\right\rangle \to | \bar{0}\left.\right\rangle\) and \(| 1\left.\right\rangle \to | \bar{1}\left.\right\rangle\). As a result, all the code states can be prepared from a properly superposed initial state of the Fock ground and first-excited states with a single-state engineering protocol. This scheme is illustrated by the “middle-to-right” process in Fig. 3a.

Third, we introduce the code space transformation. The binomial code states \(| {\bar{0}}_{b}\left.\right\rangle\) and \(| {\bar{1}}_{b}\left.\right\rangle\), as defined in Eq. (10), exactly satisfy the Knill-Laflamme condition for single-photon loss errors, cf. Eq. (11). In principle, binomial code states corrupted by errors can be perfectly recovered. However, this correction process requires sophisticated unitary operations that facilitate state transfers between the logical code words and the error words20. Implementing such recovery operations remains a challenge with state-of-the-art technological capabilities20,73. In contrast, the cat code words \(| {\bar{0}}_{c}\left.\right\rangle\) and \(| {\bar{1}}_{c}\left.\right\rangle\) defined by Eq. (8) satisfy the Knill-Laflamme condition for single-photon loss error on some “sweet spots” determined by Eq. (9). Such cat codes at sweet spots offer the important advantage that the corrupted state will automatically return to itself after four photon losses (see more details in Supplementary Note 4). By tracking the results of parity measurements and updating our knowledge on the code basis, we can implement the autonomous quantum error correction (AQEC) protocol against single-photon loss errors without feedback operations. Therefore, we introduce the state transformation between the binomial and cat code spaces on demand without deforming the encoded quantum information (preserving the superposition coefficients of the encoded state), i.e.,

$$| {\psi }_{b}\left.\right\rangle ={c}_{0}| {\bar{0}}_{b}\left.\right\rangle +{c}_{1}| {\bar{1}}_{b}\left.\right\rangle \to | {\psi }_{c}\left.\right\rangle ={c}_{0}| {\bar{0}}_{c}\left.\right\rangle +{c}_{1}| {\bar{1}}_{c}\left.\right\rangle .$$
(15)

Such a code space transformation process can be achieved by designing the following target transition Hamiltonian

$${\hat{H}}_{T}(t)=-\Delta \left(| {\bar{0}}_{t}\left.\right\rangle \langle {\bar{0}}_{t}| +| {\bar{1}}_{t}\rangle \left\langle \right.{\bar{1}}_{t}| \right),$$
(16)

with time-dependent eigenstates

$$\left\{\begin{array}{l}| {\bar{0}}_{t}\left.\right\rangle =\sqrt{1-h(t)}| {\bar{0}}_{b}\left.\right\rangle +\sqrt{h(t)}| {\bar{0}}_{c}\left.\right\rangle ,\quad \\ | {\bar{1}}_{t}\left.\right\rangle =\sqrt{1-h(t)}| {\bar{1}}_{b}\left.\right\rangle +\sqrt{h(t)}| {\bar{1}}_{c}\left.\right\rangle .\quad \end{array}\right.$$
(17)

Here, h(t) is a slow enough time-varying function with h(t0) = 0 and h(tf) = 1 such that

$$\left\{\begin{array}{l}{\hat{H}}_{{{{\rm{bin}}}}}\equiv {\hat{H}}_{T}({t}_{0})=-\Delta \left(| {\bar{0}}_{b}\left.\right\rangle \langle {\bar{0}}_{b}| +| {\bar{1}}_{b}\rangle \left\langle \right.{\bar{1}}_{b}| \right)\quad \\ {\hat{H}}_{{{{\rm{cat}}}}}\equiv {\hat{H}}_{T}({t}_{f})=-\Delta \left(| {\bar{0}}_{c}\left.\right\rangle \langle {\bar{0}}_{c}| +| {\bar{1}}_{c}\rangle \left\langle \right.{\bar{1}}_{c}| \right).\quad \end{array}\right.$$
(18)

Here, t0 and tf are the initial and final time moments, respectively. As illustrated in Fig. 3b, this process allows any encoded state in the binomial code space to be adiabatically transformed to the corresponding encoded state in the cat code space.

Last, we introduce the method of Floquet Hamiltonian engineering. In order to generate the target Hamiltonian \({\hat{H}}_{T}\), which is in general an arbitrary function of the quadrature operators \(\hat{x}\) and \(\hat{p}\), we drive the cavity by a periodic external potential \(V(\hat{x},t)=V(\hat{x},t+{T}_{d})\) with Td = 2π/ωd, i.e.,

$$\hat{{{{\mathcal{H}}}}}(t)=\frac{{\omega }_{0}}{2}\left({\hat{p}}^{2}+{\hat{x}}^{2}\right)+\beta V(\hat{x},t).$$
(19)

A periodically driven system is also called a Floquet system78,79. By transforming the above Hamiltonian into the rotating frame of frequency Ω = 2π/T with \(T=n{T}_{d}(n\in {{\mathbb{Z}}}^{+})\), we have \(\hat{O}(t)\hat{x}{\hat{O}}^{{{\dagger}} }(t)=\hat{x}\cos (\Omega t)+\hat{p}\sin (\Omega t)\) with time-evolution operator \(\hat{O}(t)\equiv {e}^{i{\hat{a}}^{{{\dagger}} }\hat{a}\Omega t}\). The transformed Hamiltonian in the rotating frame is given by

$$\begin{array}{rcl}\hat{H}(t)&\equiv &\hat{O}(t)\hat{{{{\mathcal{H}}}}}(t){\hat{O}}^{{{\dagger}} }(t)-i\lambda \hat{O}(t){\dot{\hat{O}}}^{{{\dagger}} }(t)\\ &=&\beta V\left[\hat{x}\cos (\Omega t)+\hat{p}\sin (\Omega t),t\right].\end{array}$$
(20)

Here, we have adapted the multi-photon resonance condition T = 2π/ω0 or equivalently Ω = ω0, i.e., the driving frequency is set to be n times the bare frequency of the harmonic oscillator.

The Flouqet theorem states that the stroboscopic time evolution of a periodic time-varying system is described by a time-independent Floquet Hamiltonian \({\hat{H}}_{F}\) determined by78,79,80,81,82,83

$$\exp \left(-i\frac{1}{\lambda }{\hat{H}}_{F}T\right)={{{\mathcal{T}}}}\exp \left[-i\frac{1}{\lambda }\int_{0}^{T}\hat{H}(t)dt\right],$$
(21)

where \({{{\mathcal{T}}}}\) is the time-ordering operator. Under the rotating wave approximation (RWA), the Floquet Hamiltonian \({\hat{H}}_{F}\) is just the time-averaged version of \(\hat{H}(t)\) over one Floquet period T63,82,84, i.e.,

$${\lim }_{{\omega }_{0}/\beta \to \infty }{\hat{H}}_{F}(\hat{x},\hat{p})=\frac{1}{T}\int_{0}^{T}dt\hat{H}(t).$$
(22)

By properly engineering the driving potential \(V(\hat{x},t)\)63, the Floquet Hamiltonian \({\hat{H}}_{F}(\hat{x},\hat{p})\) can be designed as the target Hamiltonian \({\hat{H}}_{T}(\hat{x},\hat{p})\).

For this purpose, we decompose a given target Hamiltonian \({\hat{H}}_{T}(\hat{x},\hat{p})\) as a sum of plane-wave operators in the noncommutative phase space63, i.e.,

$${\hat{H}}_{T}(\hat{x},\hat{p})=\frac{1}{2\pi }\int\,\int\,d{k}_{x}d{k}_{p}\,{f}_{T}({k}_{x},{k}_{p}){e}^{i({k}_{x}\hat{x}+{k}_{p}\hat{p})},$$
(23)

where the noncommutative Fourier transformation (NcFT) coefficient in Eq. (23) is given by63

$${f}_{T}({k}_{x},{k}_{p})=\frac{{e}^{\frac{\lambda }{4}({k}_{x}^{2}+{k}_{p}^{2})}}{2\pi }\int\int\,dxdp{H}_{T}^{Q}(x,p){e}^{-i({k}_{x}x+{k}_{p}p)}.$$
(24)

Here, the integrand \({H}_{T}^{Q}(x,p)=\langle \alpha | {\hat{H}}_{T}| \alpha \rangle\) is the Q-function of the target Hamiltonian with \(| \alpha \left.\right\rangle\) the coherent state defined via \(\hat{a}| \alpha \left.\right\rangle =\alpha | \alpha \left.\right\rangle\), where \(\alpha =(x+ip)/\sqrt{2\lambda }\) with \(x\equiv \langle \alpha | \hat{x}| \alpha \rangle\) and \(p\equiv \langle \alpha | \hat{p}| \alpha \rangle\).

With the NcFT coefficient, one can design the driving potential by superposing a series of cosine-type lattice potentials as ref. 63

$$V(x,\Omega t)=\int_{-\infty }^{+\infty }A(k,\Omega t)\cos [kx+\phi (k,\Omega t)]dk.\,\,$$
(25)

Here, the tunable time-dependent amplitude A(kΩt) and phase ϕ(kΩt) are given by

$$\left\{\begin{array}{l}A(k,t)=k\left\vert {f}_{T}(k\cos \Omega t,k\sin \Omega t)\right\vert \hfill \quad \\ \phi (k,t)=\,{{\mbox{Arg}}}\,\left[{f}_{T}(k\cos \Omega t,k\sin \Omega t)\right],\quad \end{array}\right.$$
(26)

where we have adopted \({k}_{x}=k\cos \Omega t\) and \({k}_{p}=k\sin \Omega t\). Each cosine component can be implemented with, e.g., an optical lattice that is formed by laser beams intersecting at an angle in cold-atom experiments85,86,87 or a JJ potential in superconducting circuits88,89,90. Note that, according to the definitions given in Eqs. (19), (20), (22), (23), and (25), the Floquet and target Hamiltonians actually differ by an overall prefactor, i.e., \({\hat{H}}_{F}=\beta {\hat{H}}_{T}\).

Following the adiabatic ramp protocol61, we prepare a target code state according to the Schrödinger equation

$$i\lambda \frac{\partial }{\partial t}| {\psi }_{{{{\rm{pre}}}}}(t)\left.\right\rangle ={\hat{H}}_{{{{\rm{adia}}}}}(t)| {\psi }_{{{{\rm{pre}}}}}(t)\left.\right\rangle ,$$
(27)

where the Hamiltonian \({\hat{H}}_{{{{\rm{adia}}}}}(t)\) is given by

$${\hat{H}}_{{{{\rm{adia}}}}}(t)=\frac{{\omega }_{0}}{2}\left({\hat{p}}^{2}+{\hat{x}}^{2}\right)+\beta (t)V[\hat{x},\Omega (t)t]$$
(28)

with time-dependent amplitude β(t) and frequency Ω(t) for the designed driving potential \(V(\hat{x},\Omega t)\), cf. Eq. (25). The amplitude β(t) [frequency Ω(t)] of the driving potential is tuned from β(t0) = 0 [Ω(t0) ≠ ω0] to β(tf) = βf [Ω(tf) = ω0] over a preparation time tf with sigmoidal modulations

$$\left\{\begin{array}{l}\beta (t)=\frac{{\beta }_{f}}{{Z}_{1}\left[1+{e}^{-{s}_{1}(t-{t}_{c,1})}\right]}-\frac{{\beta }_{f}}{1+{e}^{{s}_{1}{t}_{c,1}}}, \hfill \quad \\ \Omega (t)=\Omega (0)+\frac{{\omega }_{0}-\Omega (0)}{{Z}_{2}\left[1+{e}^{-{s}_{2}(t-{t}_{c,2})}\right]}-\frac{{\omega }_{0}}{1+{e}^{{s}_{2}{t}_{c,2}}},\quad \end{array}\right.$$
(29)

where \({Z}_{j = 1,2}={[1+{e}^{-{s}_{j}({t}_{f}-{t}_{c,j})}]}^{-1}-{(1+{e}^{{s}_{j}{t}_{c,j}})}^{-1}\) with tc,j and sj being the time centers and slopes of the two adiabatic processes, respectively.

The adiabaticity condition for Floquet systems is given by ref. 61,

$$\frac{| \langle {\psi }_{m}^{F}(s)| \frac{d}{ds}| {\psi }_{n}^{F}(s)\rangle | }{| 1-{e}^{-i({\epsilon }_{m}-{\epsilon }_{n})T}| }\ll \frac{{t}_{f}}{T},\quad m\ne n.$$
(30)

Here, s = t/tf is the parameterized time during the adiabatic ramp, and \(| {\psi }_{n}^{F}(s)\left.\right\rangle\) is the Floquet eigenstate of Hamiltonian (20) with quasienergies ϵn. Similar to the case of static systems, the above adiabaticity condition can be guaranteed by extending the preparation time tf. However, in Floquet systems, the adiabatic evolution may break down when the difference between two Floquet states ϵm − ϵn is modulo 2π/T. As a result, an initial detuning Ω(0) − ω0 > 0, which must be incommensurate with the cavity frequency, is required to avoid resonantly driving the linear oscillator. This ensures that the denominator in the above adiabaticity condition does not vanish during the early stage of the protocol. This detuning opens a gap in the Floquet spectrum that can suppress unwanted excitations during the adiabatic preparation process61. Once the ramped driving field introduces sufficient nonlinearity to the cavity mode, it can be safely tuned to be resonant with the cavity, i.e., Ω(tf) → ω0.

Lattice gates decomposition

For a given target Hamiltonian, each cosine component in the engineered driving potential, cf. Eqs. (25)–(26), corresponds to a JJ-based device (e.g., a SQUID or SNAIL) in a superconducting circuit architecture61 or an optical lattice in a cold-atom experimental platform87. The engineered driving potential Eq. (25) can be approached via a set of discretized cosine-lattice potentials as follows

$$V(x,t)\approx {\sum }_{n=0}^{N}A({k}_{n},\Omega t)\Delta k\cos \left[{k}_{n}x+\phi ({k}_{n},\Omega t)\right],\,\,\,$$
(31)

where kn = nΔk is the discretized wavenumber with n labeling the discretization step. In principle, one can use N cosine potentials to synthesize the desired driving potential63. However, it is challenging to precisely manipulate multiple cosine-lattice potentials simultaneously in experiments. Here, we propose to engineer code states with a single cosine-lattice potential that implements lattice gate operations, cf. Eq. (4). In this section, we provide an analytical framework for decomposing the code state engineering process into a sequence of quantum lattice gates.

We first introduce the grid lattice gate. For a given code state engineering process, the driving amplitude A(knt) and phase ϕ(knt) in Eq. (31) can be calculated according to Eqs. (24) and (26), as illustrated by a chart in Fig. 4a over a single Floquet period T = 2π/ω0. The chart is divided into small grids with the height and width given by the wavenumber interval \(\Delta k={k}_{\max }/N\) and the time step Δt = T/M (\(N,M\in {{\mathbb{Z}}}^{+}\)), respectively. Each grid in the chart, labeled by {kn = nΔktm = mΔt}, corresponds to a discrete unitary time evolution in the rotating frame, cf. Eq. (20), i.e.,

$$| \psi ({k}_{n},{t}_{m}+\Delta t)\left.\right\rangle ={{{\rm{PSL}}}}({k}_{n},{t}_{m})| \psi ({k}_{n},{t}_{m})\left.\right\rangle ,$$
(32)

where PSL(kntm) ≡ PSL(ζnmσnmγnmδnm) is the lattice gate defined in Eq. (4) with gate parameters

$$\left\{\begin{array}{ll}{\zeta }_{nm}={k}_{n}\cos \Omega {t}_{m}, \hfill &{\sigma }_{nm}={k}_{n}\sin \Omega {t}_{m},\\ {\gamma }_{nm}=-\frac{1}{\lambda }\beta A({k}_{n},\Omega {t}_{m})\Delta k\Delta t,&{\delta }_{nm}=\phi \left({k}_{n},\Omega {t}_{m}\right).\end{array}\right.$$
(33)

We point out that each discrete unitary time evolution PSL(kntm) can be viewed as the operation of a grid lattice gate, where the original time parameter just plays the role of a control parameter in this scheme. All gate parameters can be independently tuned based on the restriction condition given by Eq. (33).

Fig. 4: Quantum lattice gate decomposition.
figure 4

a Chart of the driving field modulation over a single Floquet period T = 2π/ω0, showing the amplitude A(knt) or the phase ϕ(knt) as a function of discrete momentum and time. The whole chart is divided into N × M grids. Each grid is labeled by {kn = nΔktm = mΔt}, corresponding to a discrete unitary evolution generated by the phase-space lattice (PSL) gate PSL(kntm) that is defined in Eqs. (32)–(33). The tm-PSL gate is composed by concatenating grid lattice gates along the row for a fixed time tm, cf. Eq. (34). b Quantum circuits of the lattice gate decomposition: (upper) the whole code state engineering process is decomposed into a sequence of Ξ(βΩ) lattice gates, cf. Eq. (36); (middle) each Ξ(βΩ) gate is decomposed into a sequence of tm-PSL gates, cf. Eq. (34); (lower) a single tm-PSL gate is decomposed into a sequence of grid lattice gates PSL(kntm) shown in (a).

We then discuss how to design the gate sequence of a quantum circuit. By concatenating the grid lattice gates with discretized wavenumbers for a fixed time step tm = mΔt, as shown in Fig. 4a, we define the lattice gate

$${{{{\rm{t}}}}}_{{{{\rm{m}}}}}-{{{\rm{PSL}}}} \equiv {\prod }_{n=0}^{N}{{{\rm{PSL}}}}({k}_{n},{t}_{m})\\ = {\prod }_{n=0}^{N}{e}^{-\frac{i}{\lambda }\beta A({k}_{n},\Omega {t}_{m})\Delta k\Delta t\cos [{\zeta }_{nm}\hat{x}+{\sigma }_{nm}\hat{p}+\phi ({k}_{n},\Omega {t}_{m})]}\\ \approx {e}^{-\frac{i}{\lambda }\Delta t{\sum }_{n = 0}^{N}\beta A({k}_{n},\Omega {t}_{m})\Delta k\cos [{\zeta }_{nm}\hat{x}+{\sigma }_{nm}\hat{p}+\phi ({k}_{n},\Omega {t}_{m})]}.$$
(34)

In the last step, we have adapted the Lie-Trotter product formula91, i.e., \({e}^{i(A+B)}=\mathop{\lim }_{n\to \infty }{({e}^{i\frac{1}{n}A}{e}^{i\frac{1}{n}B})}^{n}\), for sufficiently small time intervals. Compared to Eq. (31) in the rotating frame, cf. Eq. (20), the tm-PSL(βΩ) lattice gate faithfully realizes the time evolution with the desired driving potential at each time interval with a single cosine-potential device.

By repeatedly applying the tm-PSL gate operations at different time steps and concatenating them, we realize the following total gate operation over one full Floquet period

$$\Xi (\beta ,\Omega )\equiv {\prod }_{m=0}^{M}{{{{\rm{t}}}}}_{{{{\rm{m}}}}}-{{{\rm{PSL}}}}(\beta ,\Omega ),$$
(35)

which is referred to as the Ξ-type lattice gate in this paper. Here, we have explicitly written the Floquet gate with control parameters β and Ω. In this way, one can prepare the target state by concatenating Ξ(βΩ) gates

$$| {\psi }_{f}\left.\right\rangle ={\prod }_{s=0}^{{t}_{f}/T}\Xi [\beta (sT),\Omega (sT)]| {\psi }_{i}\left.\right\rangle ,$$
(36)

where the adiabatic parameters β(sT) and Ω(sT) with the stroboscopic step index s = 0, 1,  , tf/T are determined from the adiabatic ramp protocol, cf. Eq. (28).

In the upper panel of Fig. 4b, we show the quantum circuits of engineering a sequence of Ξ-type lattice gates for a cavity. The middle panel shows that each Ξ(βΩ) gate is decomposed into a sequence of tm-PSL(βΩ) lattice gates. In the lower panel, we show the quantum circuit for realizing the tm-PSL(βΩ) gate with a sequence of grid lattice gates PSL(kntm).

Results

In this section, we apply our method to concrete examples of bosonic code state engineering. We will first discuss how to prepare a single binomial code state and embed the binomial code space. Then we will show the transformation from binomial codes to cat codes. Finally, we discuss the autonomous quantum error correction of the obtained cat code states against single-photon losses.

Single code state preparation

As the first application of our method, we aim to prepare an arbitrary superposition state of binomial code words given by

$$| {\psi }_{T}\left.\right\rangle ={c}_{0}| {\bar{0}}_{b}\left.\right\rangle +{c}_{1}| {\bar{1}}_{b}\left.\right\rangle .$$
(37)

According to Eq. (13), we construct the target Hamiltonian by \({\hat{H}}_{T}=-\Delta | {\psi }_{T}\left.\right\rangle \left\langle \right.{\psi }_{T}|\). From Eq. (24), we calculate the NcFT coefficient of the target Hamiltonian as follows (see detailed derivation in Supplementary Note 1. B)

$${f}_{T}(k,\tau )= -\frac{\Delta }{4}\left[| {c}_{0}{| }^{2}{f}_{00}+3| {c}_{1}{| }^{2}{f}_{22}+3| {c}_{0}{| }^{2}{f}_{44}+| {c}_{1}{| }^{2}{f}_{66}\right.\\ +\sqrt{3}\left(| {c}_{0}{| }^{2}{f}_{04}+| {c}_{1}{| }^{2}{f}_{26}+{c}_{0}{c}_{1}^{* }{f}_{02}+{c}_{0}{c}_{1}^{* }{f}_{46}+\,{\mbox{c.c.}}\,\right)\\ +\left.\left(3{c}_{1}{c}_{0}^{* }{f}_{24}+{c}_{0}{c}_{1}^{* }{f}_{06}+\,{\mbox{c.c.}}\,\right)\right]$$
(38)

with the explicit expression of the Fourier component fnm given by Eq. (S10) in Supplementary Note 1. B.

Both the general superposition state \(| {\psi }_{T}\left.\right\rangle\) and the corresponding target Hamiltonian \({\hat{H}}_{T}\) are invariant under the two-fold discrete phase-space rotation given by the parity operator \(\hat{{{{\mathcal{P}}}}}=\exp (i\pi \hat{n})\)38,92, cf. Eq. (12). As an example, we choose the superposition coefficients of the target state as c0 = 1/2 and \({c}_{1}=\sqrt{3}/2\) in Eq. (37). In Fig. 5a, we plot the Q-function of \({\hat{H}}_{T}\) in phase space, i.e., \({H}_{T}^{Q}(x,p)=\langle \alpha | {\hat{H}}_{T}| \alpha \rangle ,\) which clearly shows a two-fold rotational symmetry in phase space. The engineered periodically driving potential can be straightforwardly calculated from Eqs. (25), (26) and (38). In Fig. 5b, we plot the time variation of the amplitude A(kt) and phase ϕ(kt) for the engineered driving potential, cf. Eqs. (25) and (26).

Fig. 5: State preparation of a single binomial code.
figure 5

a Q-function of the target Hamiltonian \({H}_{T}^{Q}(x,p)=-\Delta \langle \alpha | {\psi }_{T}\rangle \langle {\psi }_{T}| \alpha \rangle\) with \(| {\psi }_{T}\left.\right\rangle =(| {\bar{0}}_{b}\left.\right\rangle +\sqrt{3}| {\bar{1}}_{b}\left.\right\rangle )/2\). b Charts of the driving amplitude A(kt) (left) and the driving phase ϕ(kt) (right) for the engineered driving potential, cf. Eqs. (25) and (26), for the target Hamiltonian \(| {\psi }_{T}\left.\right\rangle\). c Snapshots of the Wigner functions W(xp) of the prepared state during the adiabatic ramp process at different time moments. d Time evolution of the infidelity of the prepared state with respect to the target state \(| {\psi }_{T}\left.\right\rangle\) given by \(1-{F}_{{{{\rm{pre}}}}}(t)\), where the fidelity \({F}_{{{{\rm{pre}}}}}(t)\) is calculated from Eq. (39). e Envelopes of driving amplitude β(t) and driving frequency Ω(t) during the adiabatic ramp process, cf. Eqs. (28)–(29). Parameters: λ = 0.25, Δ = 1.3ω0, βf = 0.02ω0, Ω(0) = ω0/(1 + π × 10−3), s1 = 40/tf, s2 = 30/tf, tc,1 = tf/6, tc,2 = 2tf/3.

We then adiabatically ramp the driving potential following Eq. (29) to prepare the target state. In Fig. 5c, we plot the snapshots of the Wigner functions of the prepared state during the adiabatic ramp process. To quantify the performance of state preparation, we define the fidelity93 of the prepared state with respect to the target code state by

$${F}_{{{{\rm{pre}}}}}(t)\equiv \,{{\mbox{Tr}}}\,\left[\sqrt{{\hat{\rho }}_{T}^{1/2}{\hat{\rho }}_{{{{\rm{pre}}}}}(t){\hat{\rho }}_{T}^{1/2}}\right],$$
(39)

where \({\hat{\rho }}_{{{{\rm{pre}}}}}=| {\psi }_{{{{\rm{pre}}}}}\left.\right\rangle \left\langle \right.{\psi }_{{{{\rm{pre}}}}}|\) is the density operator of the prepared state and \({\hat{\rho }}_{T}=| {\psi }_{T}\left.\right\rangle \left\langle \right.{\psi }_{T}|\) is the density operator of the target code state. This fidelity measures how closely the prepared state approximates the target state. For pure states, the above defined fidelity reduces to \({F}_{{{{\rm{pre}}}}}(t)=| \langle {\psi }_{{{{\rm{pre}}}}}(t)| {\psi }_{T}\rangle |\). In Fig. 5d, we plot the infidelity of the prepared state \(1-{F}_{{{{\rm{pre}}}}}(t)\) at stroboscopic time steps during the whole preparation process. The early oscillating behavior of infidelity arises from the initially detuned driving potential that rotates the prepared states in the phase space. This oscillation ceases once the driving potential becomes resonant with the cavity mode.

Unitary errors

In practice, unavoidable unitary errors can occur during the state preparation process. Here, we analyze the possible unitary errors in the code state engineering process. The first unitary error source comes from the discretization of implementing the driving potential, cf. in Eq. (31), with finite wavenumber interval and time step. In Fig. 6a, we plot the infidelity of the prepared state as a function of time for different numbers of wavenumber steps N and the number of time steps M. As expected, fewer wavenumber and time steps result in larger discretization errors, as shown, e.g., by the result for N = 10 and M = 10. As the number of discretized steps increases, the errors are quickly suppressed and eventually remain almost unchanged for NM > 20, see the final infidelity of prepared state for (N = 20, M = 20) and (N = 50, M = 50). Hereafter in this paper, unless otherwise specified, we set the wavenumber and time steps to NM = 20.

Fig. 6: Unitary errors.
figure 6

a Infidelity of the prepared state \(1-{F}_{{{{\rm{pre}}}}}(t)\), where the fidelity \({F}_{{{{\rm{pre}}}}}(t)\) is calculated from Eq. (39), as a function of time for different wavenumber steps N and time steps M. b Final infidelity of the prepared state, i.e., \(1-{F}_{{{{\rm{pre}}}}}({t}_{f})\), as a function of energy gap Δ for NM = 20. Other system parameters are set to be the same as those in Fig. 5.

The second unitary error source is the non-adiabatic excitation that causes leakage of code state to the excited error states during the finite-time adiabatic ramp process. This type of error can, in principle, be suppressed by using a larger energy gap Δ. To reveal the effect of the energy gap Δ, we plot the final infidelity of the prepared state as a function of the energy gap in Fig. 6b. Indeed, it shows that the infidelity is relatively high for a small gap but drops dramatically to the minimum around the optimal gap value Δ ≈ 1.3ω0. However, the infidelity increases gradually again as the gap continues increasing beyond the optimal value. This is because our method relies on the RWA (valid for βω0). As the gap of the Floquet spectrum increases with the driving strength β, i.e., \({\hat{H}}_{F}=\beta {\hat{H}}_{T}\), a large enough energy gap will break the RWA and thus affect the validity of Eq. (22). The higher-order harmonics of the rotating-frame Hamiltonian \({\hat{H}}_{l}(\hat{x},\hat{p})={T}^{-1}\int_{0}^{T}\hat{H}(t){e}^{-il\Omega t}dt\) deviate from the engineered Hamiltonian in the Floquet-Magnus expansion. Specifically, the Floquet Hamiltonian in Eq. (21) is given by \({\hat{H}}_{F}={\sum }_{n = 0}^{\infty }{(\beta /\Omega )}^{n}{\hat{H}}_{F}^{(n)}(\hat{x},\hat{p})\), where \({\hat{H}}_{F}^{(0)}(\hat{x},\hat{p})\) is the RWA term given by Eq. (22) and the non-RWA terms \({\hat{H}}_{F}^{(n\ge 1)}(\hat{x},\hat{p})\) are functions of the higher-order harmonics \({\hat{H}}_{l}(\hat{x},\hat{p})\) that can be obtained from the standard Floquet-Magnus expansion via a recursive procedure84,94,95. To suppress unitary non-RWA errors, one can introduce additional higher-order driving potentials, i.e., \(V(\hat{x},t)={\sum }_{n = 0}{\beta }^{n}{V}^{(n)}(\hat{x},t)\), to approximate the target Hamiltonian \({\hat{H}}_{T}(\hat{x},\hat{p})\) to the desired order of precision96. In this work, we restrict our analytical calculation to the lowest order of the Flouqet-Magnus expansion.

Code space embedding

As the second application, we aim to embed a finite binomial code space into the infinite Fock space of a cavity. According to Eq. (14), we set the target Hamiltonian as

$${\hat{H}}_{T}=-\Delta \left(| {\bar{0}}_{b}\left.\right\rangle \langle {\bar{0}}_{b}| +| {\bar{1}}_{b}\rangle \left\langle \right.{\bar{1}}_{b}| \right).$$
(40)

In Fig. 7(a), we calculate and plot the Q-function \({H}_{Q}^{T}\) of target Hamiltonian \({H}_{T}^{Q}(x,p)=\langle \alpha | {\hat{H}}_{T}| \alpha \rangle\). In contrast to the case of single code state preparation, the current target Hamiltonian Q-function is invariant under four-fold discrete rotations in phase space38,92. The NcFT Fourier coefficient is given by (see details in Supplementary Note 1. B)

$${f}_{T}(k,\tau )=-\frac{\Delta }{4}\left[{f}_{00}+3{f}_{22}+3{f}_{44}+{f}_{66}+\sqrt{3}\left(\,{f}_{04}+{f}_{26}+\,{\mbox{c.c.}}\,\right)\right],$$
(41)

which yields the driving field with the charts of periodically modulated amplitude and phase shown in Fig. 7b.

Fig. 7: Code space embedding.
figure 7

a Q-function of the target Hamiltonian given by Eq. (40). b Charts of the driving amplitude A(kt) and the driving phase ϕ(kt) for the target Hamiltonian. c Snapshots of the Wigner functions W(xp) of the prepared states at different time moments: (upper) binomial logical state \(| {\psi }_{T}\left.\right\rangle =| {\bar{0}}_{b}\left.\right\rangle\), (middle) binomial logical state \(| {\psi }_{T}\left.\right\rangle =| {\bar{1}}_{b}\left.\right\rangle\) and (lower) binomial superposition state \(| {\tilde{\psi }}_{T}\left.\right\rangle =(| {\bar{0}}_{b}\left.\right\rangle +{e}^{i\Delta \varphi }| {\bar{1}}_{b}\left.\right\rangle )/\sqrt{2}\), cf. Eq. (43). d Infidelities of three prepared states with respect to the corresponding target states as a function of time during the adiabatic ramp process. For the superposition state, the target state has been updated according to Eq. (43). Parameters: λ = 0.25, Δ = 1.4ω0, βf = 0.02ω0, Ω(0) = ω0/(1 + π × 10−3), s1 = 40/tf, s2 = 30/tf, tc,1 = tf/6, tc,2 = 2tf/3.

As claimed above, the code space embedding is supposed to enable preparing any superposition state with the same adiabatically ramped driving protocol. In the upper and middle rows of Fig. 7c, we show the Wigner functions of two prepared logical code states \(| {\bar{0}}_{b}\left.\right\rangle\) and \(| {\bar{1}}_{b}\left.\right\rangle\) at different moments, starting from the cavity vacuum state \(| 0\left.\right\rangle\) and the second excited Fock state \(| 2\left.\right\rangle\), respectively. In Fig. 7d, we also plot the time evolution of infidelities for both binomial code words. It is clear that the prepared state eventually shows a very high fidelity as expected. In the lower row of Fig. 7c, we show the Wigner functions of a prepared superposition state starting from the superposition Fock state \(| {\psi }_{i}\left.\right\rangle ={c}_{0}| 0\left.\right\rangle +{c}_{1}| 2\left.\right\rangle\) (with superposition coefficients \({c}_{0}={c}_{1}=1/\sqrt{2}\)) at different moments.

It is important to note that the two prepared code words (i.e., \(| {\bar{0}}_{b}\left.\right\rangle\) and \(| {\bar{1}}_{b}\left.\right\rangle\)) may accumulate different phase factors during the adiabatic preparation process, i.e.,

$$| {\psi }_{i}\left.\right\rangle \to | {\psi }_{{{{\rm{pre}}}}}\left.\right\rangle ={c}_{0}{e}^{i{\varphi }_{0}}| {\bar{0}}_{b}\left.\right\rangle +{c}_{1}{e}^{i{\varphi }_{1}}| {\bar{1}}_{b}\left.\right\rangle .$$
(42)

Each phase of the code words includes both dynamical and geometric contributions. Since the embedded code space is degenerate, the phase difference between the two code words Δφ = φ1 − φ0 is only determined by the geometric phases, which can be deterministically extracted via \(\Delta \varphi =\arg (\langle {\bar{1}}_{b}| {\psi }_{{{{\rm{pre}}}}}\rangle \langle {\psi }_{{{{\rm{pre}}}}}| {\bar{0}}_{b}\rangle )\) from numerical simulation. In our case, the extracted geometric phase difference is Δφ ≈ 0.27π, and it remains independent of the superposition coefficients c0 and c1 within our numerical precision. This phase difference also accounts for the slight asymmetry observed in the final Wigner function, cf. the lower row of Fig. 7c. By updating our knowledge of the final target state with the extracted phase difference, i.e.,

$$| {\tilde{\psi }}_{T}\left.\right\rangle \to {c}_{0}| {\bar{0}}_{b}\left.\right\rangle +{c}_{1}{e}^{i\Delta \varphi }| {\bar{1}}_{b}\left.\right\rangle ,$$
(43)

the prepared states consistently exhibit high fidelities, as shown in Fig. 7d.

Code spaces transformation

As the final application, we study the code space transformation between the binomial and cat code states by engineering the transition Hamiltonian given by Eqs. (16)–(17). In Supplementary Note 2, we calculate explicitly the Husmi-Q function and the NcFT coefficients for the transition Hamiltonian given by Eq. (16). We consider the slow enough time-varying function h(t) in Eq. (17) as follows

$$h(t)=\sin {\left[\frac{\pi }{2}{\sin }^{2}\left(\frac{\pi t}{2{t}_{f}}\right)\right]}^{2}.$$
(44)

This function meets the boundary conditions h(0) = 0 and h(tf) = 1, enabling a state transfer process in the adiabatic limit.

In Fig. 8, we first present numerical simulations of the transformation processes from the binomial code state \(| {\bar{0}}_{b}\left.\right\rangle\) (\(| {\bar{1}}_{b}\left.\right\rangle\)) to the cat code state \(| {\bar{0}}_{c}\left.\right\rangle\) (\(| {\bar{1}}_{c}\left.\right\rangle\)). The infidelity behaviors of these processes are respectively depicted by the black and red curves in Fig. 8a, while the time evolution of h(t) is shown in Fig. 8b. Note that the fidelities exhibit relatively high values even at the beginning of the transformation processes, which can be attributed to the significant overlap between the initial binomial state and the corresponding cat state at the first “sweet spot” (α2 ≈ 2.34) determined by Eq. (9). This overlap is also reflected in the similarity between the Wigner functions of binomial and cat code states shown in Fig. 2. Despite this large overlap, the fidelities continue to improve throughout the transformation process, reaching even higher values as they approach their respective target states.

Fig. 8: Code space transformation.
figure 8

a Infidelity time evolution of the prepared code states with respect to the cat code basis states \(| {\psi }_{T}\left.\right\rangle =| {\bar{0}}_{c}\left.\right\rangle\), \(| {\psi }_{T}\left.\right\rangle =| {\bar{1}}_{c}\left.\right\rangle\) and the time-dependent updated superposition state \(| {\tilde{\psi }}_{T}(t)\left.\right\rangle =(| {\bar{0}}_{c}\left.\right\rangle +{e}^{i\Delta \varphi (t)}| {\bar{1}}_{c}\left.\right\rangle )/\sqrt{2}\), cf. Eq. (45). b Envelope of the time-varying function h(t) introduced in the transition Hamiltonian, cf. Eq. (17). Parameters: λ = 0.25, tf = 5000 × 2π/ω0, Δ = 0.1ω0.

As discussed above, during the adiabatic transition processes, the two logical code words may accumulate different phases. As a result, an initial superposition binomial code state \(| {\psi }_{0}\left.\right\rangle ={c}_{0}| {\bar{0}}_{b}\left.\right\rangle +{c}_{1}| {\bar{1}}_{b}\left.\right\rangle\) is transferred to the following code state

$$| {\psi }_{t}\left.\right\rangle ={c}_{0}{e}^{i{\varphi }_{0}(t)}| {\bar{0}}_{t}\left.\right\rangle +{c}_{1}{e}^{i{\varphi }_{1}(t)}| {\bar{1}}_{t}\left.\right\rangle ,$$

where the transition code words \(| {\bar{0}}_{t}\left.\right\rangle\) and \(| {\bar{1}}_{t}\left.\right\rangle\) are given by Eq. (17). The time-dependent accumulated phases φ0(t) and φ1(t) include both dynamical and geometric contributions. Due to the designed form of the transition Hamiltonian given by Eq. (16), the transition code words \(| {\bar{0}}_{t}\left.\right\rangle\) and \(| {\bar{1}}_{t}\left.\right\rangle\) are degenerate. Thus, the phase difference Δφ(t) = φ1(t) − φ0(t) only contains geometric contribution and can be deterministically extracted from the numerical simulation via \(\Delta \varphi (t)=\arg (\langle {\bar{1}}_{t}| {\psi }_{t}\rangle \langle {\psi }_{t}| {\bar{0}}_{t}\rangle )\), which is independent of the superposition coefficients c0 and c1. By updating the target state with the extracted phase difference as a function of time (up to a global phase factor), i.e.,

$$| {\tilde{\psi }}_{T}(t)\left.\right\rangle \propto {c}_{0}| {\bar{0}}_{c}\left.\right\rangle +{c}_{1}{e}^{i\Delta \varphi (t)}| {\bar{1}}_{c}\left.\right\rangle ,$$
(45)

we display in Fig. 8a the infidelity of the transferred superposition code state (green curve) exhibiting high fidelities after the transition process.

Autonomous quantum error correction

QEC schemes with binomial codes typically require active error-syndrome measurements and adaptive recovery operations, which are hardware-intensive and prone to propagating errors of gate operations. In contrast, cat codes allow for autonomous quantum error correction schemes that can protect the encoded quantum information by autonomously correcting photon-loss errors based on parity measurements, without the need for feedback operations97. Here, we study in detail how to implement AQEC against single-photon loss errors using four-legged cat states based on our driving protocol.

To investigate the impact of a noisy environment and the error correction process, we extend the unitary Schrödinger equation, cf. Eq. (27), to the Lindblad master equation in the rotating frame as follows

$$\frac{\partial }{\partial t}{\hat{\rho }}_{c}(t)=\quad -\frac{i}{\lambda }\left[{\hat{H}}_{EC},{\hat{\rho }}_{c}(t)\right]+\kappa {{{\mathcal{L}}}}[\hat{a}]{\hat{\rho }}_{c}(t),$$
(46)

where \({\hat{H}}_{EC}\) is the Hamiltonian used for error correction and \({\hat{\rho }}_{c}(t)=| {\psi }_{c}(t)\left.\right\rangle \left\langle \right.{\psi }_{c}(t)|\) is the density operator of the code state. The second term on the right-hand side of Eq. (46) describes the decay of the cavity with a single-photon loss rate κ, where \({{{\mathcal{L}}}}[\hat{O}]\) is the Lindblad superoperator defined via \({{{\mathcal{L}}}}[\hat{O}]\hat{\rho }\equiv\hat{O}\hat{\rho }{\hat{O}}^{{{\dagger}} }-({\hat{O}}^{{{\dagger}} }\hat{O}\hat{\rho }+\hat{\rho }{\hat{O}}^{{{\dagger}} }\hat{O})/2\).

Following our method, to protect the quantum information encoded with cat states, we set the target Hamiltonian to be

$${\hat{H}}_{T}=-\Delta \left(| {\bar{0}}_{c}\left.\right\rangle \langle {\bar{0}}_{c}| +| {\bar{1}}_{c}\rangle \langle {\bar{1}}_{c}| +| {\bar{0}}_{e}\rangle \langle {\bar{0}}_{e}| +| {\bar{1}}_{e}\rangle \left\langle \right.{\bar{1}}_{e}| \right).$$
(47)

Here, \(| {\bar{0}}_{c}\left.\right\rangle\) and \(| {\bar{1}}_{c}\left.\right\rangle\) are the cat code words defined in Eq. (8), while \(| {\bar{0}}_{e}\left.\right\rangle \propto \hat{a}| {\bar{0}}_{c}\left.\right\rangle\) and \(| {\bar{1}}_{e}\left.\right\rangle \propto \hat{a}| {\bar{1}}_{c}\left.\right\rangle\) are the cat error states defined as

$$\left\{\begin{array}{l}| {\bar{0}}_{e}\left.\right\rangle \equiv \frac{1}{\sqrt{{{{{\mathcal{N}}}}}_{1}}}\left(| \alpha \left.\right\rangle -| -\alpha \left.\right\rangle -i| i\alpha \left.\right\rangle +i| -i\alpha \left.\right\rangle \right),\quad \\ | {\bar{1}}_{e}\left.\right\rangle \equiv \frac{1}{\sqrt{{{{{\mathcal{N}}}}}_{3}}}\left(| \alpha \left.\right\rangle -| -\alpha \left.\right\rangle +i| i\alpha \left.\right\rangle -i| -i\alpha \left.\right\rangle \right)\quad \end{array}\right.$$
(48)

with \({{{{\mathcal{N}}}}}_{m}=8{e}^{-{\alpha }^{2}}[\sinh {\alpha }^{2}+{(-1)}^{\frac{m-1}{2}}\sin {\alpha }^{2}]\) (m = 1, 3) the normalization factors, similar to \({{{{\mathcal{N}}}}}_{m}\) for m = 0, 2 in Eq. (8). The noncommutative Fourier coefficients for the cat error states are given in Supplementary Note 3.

We begin with an unknown state that is an arbitrary superposition of the cat code words

$$| {\psi }_{c}(0)\left.\right\rangle ={c}_{1}| {\bar{0}}_{c}\left.\right\rangle +{c}_{2}| {\bar{1}}_{c}\left.\right\rangle ,$$
(49)

which is the target state that we aim to protect with AQEC. The first step of our AQEC protocol is to periodically perform the projective parity measurement72,73,74

$${\hat{\Pi }}_{m}=\mathop{\sum }_{n=0}^{\infty }| 2n+m\left.\right\rangle \left\langle \right.2n+m|$$
(50)

with the parity index m = 0, 1. After each parity measurement, the code state collapses into a state with a certain parity. The probability of obtaining the m-parity outcome is given by

$${P}_{m}(t)=\,{{\mbox{Tr}}}\,[{\hat{\Pi }}_{m}{\hat{\rho }}_{c}(t){\hat{\Pi }}_{m}^{{{\dagger}} }]$$
(51)

with the constraint P0 + P1 = 1. Next, we divide the unit interval [0, 1] into two sections with lengths P0 and P1, respectively. Then one can generate a random number εr [0, 1] and identify its position within the unit interval. If the random number falls within the m-th section, the code state is projected into

$${\hat{\rho }}_{c}(t)\to \frac{{\hat{\Pi }}_{m}{\hat{\rho }}_{c}(t){\hat{\Pi }}_{m}}{{P}_{m}(t)},$$
(52)

and the fidelity is updated as

$${F}_{c}(t)\to \,{{\mbox{Tr}}}\,\left[\sqrt{{\hat{\rho }}_{T}^{1/2}{\hat{\rho }}_{c}(t){\hat{\rho }}_{T}^{1/2}}\right].$$
(53)

Here, \({\hat{\rho }}_{T}=| {\psi }_{T}\left.\right\rangle \left\langle \right.{\psi }_{T}|\) is the target density operator with \(| {\psi }_{T}\left.\right\rangle =| {\psi }_{c}(0)\left.\right\rangle\), cf. Eq. (49). The target state is updated in sequence as follows

(54)

whenever a parity change is detected (see Supplementary Note 4 for more details). By repeatedly performing parity measurements and updating the target state accordingly, the AQEC protocol is thus implemented using the cat states.

According to our proposal, the error correction Hamiltonian \({\hat{H}}_{EC}\) in Eq. (46) is set by the following time-dependent Hamiltonian, cf. Eq. (20),

$${\hat{H}}_{EC}=\hat{H}(t)=\beta V\left[\hat{x}\cos (\Omega t)+\hat{p}\sin (\Omega t),t\right],$$
(55)

where the driving field V(xt) is designed from Eqs. (25)–(26) with the target Hamiltonian given in Eq. (47). From Eqs. (34)–(35), the unitary time evolution over each Floquet period, i.e., t [sT, (s + 1)T] with \(s\in {{\mathbb{N}}}^{0}\) (nonnegative integers), can be decomposed into a sequence of quantum lattice gate operations as described by Eqs. (34) and (35).

In Fig. 9, we present the results of AQEC starting from a generic encoded state \(| {\psi }_{c}(0)\left.\right\rangle =\frac{1}{2}\left(\sqrt{\frac{5}{2}}| {\bar{0}}_{c}\left.\right\rangle +\sqrt{\frac{3}{2}}| {\bar{1}}_{c}\left.\right\rangle \right)\). The parity measurement is performed every five Floquet periods. In Fig. 9a, we plot the infidelity and parity probability P1 after each measurement for a single quantum trajectory with the Hamiltonian protection. It clearly shows that the infidelity decreases suddenly by nearly two orders of magnitude after each parity measurement, while it increases gradually between two successive measurements due to the photon loss, which causes the quantum state to evolve towards a mixed state. In Fig. 9b, we plot the infidelity of the protected state averaged over 103 trajectories. For comparison, we also show the result without the Hamiltonian protection (i.e., β = 0). It is clear that the fidelity of the protected state after each parity measurement is significantly enhanced by the Hamiltonian protection. We also plot the Wigner functions of the initial code state \(| {\psi }_{c}(0)\left.\right\rangle\) and the code states at the end of the time evolution when the updated target states after the last parity measurements recover the initial code state for different AQEC protocols, as shown by the insets of Fig. 9b.

Fig. 9: Autonomous quantum error correction.
figure 9

a Infidelity 1 − Fc(t) of the protected code state (blue solid) and the parity probability \({P}_{1}=\langle {\hat{\Pi }}_{1}\rangle\) (red solid) as functions of time for a single quantum trajectory with the driving Hamiltonian protection, i.e., \({\hat{H}}_{EC}=\hat{H}(t)\) in Eq. (55) with β = 0.02ω0. b Infidelity of the protected code state, averaged over 103 quantum trajectories, in the cases of no Hamiltonian protection (\({\hat{H}}_{EC}=0\); green solid), the ideal Hamiltonian protection [\({\hat{H}}_{EC}={\hat{H}}_{T}\) in Eq. (47) with Δ = 0.2ω0; black dashed], and the driving Hamiltonian protection [\({\hat{H}}_{EC}=\hat{H}(t)\) in Eq. (55) with β = 0.02ω0; blue solid]. In both panels, we start from the initial code state \(| {\psi }_{c}(0)\left.\right\rangle =\left(\sqrt{5}| {\bar{0}}_{c}\left.\right\rangle +\sqrt{3}| {\bar{1}}_{c}\left.\right\rangle \right)/2\sqrt{2}\) and perform parity measurements every five Floquet periods. The three insets in (b) show the Wigner functions of the code states at the beginning and end of the time evolution, illustrating how the updated target states after the last parity measurements recover the initial code state with different AQEC protocols. Other parameters: λ = 0.25, κ = 10−3ω0, wavenumber steps M = 20 and time steps N = 100, cf. Fig. 4a.

In our driving protocol, the Floquet Hamiltonian \({\hat{H}}_{F}\), corresponding to the time-periodic Hamiltonian in Eq. (55), approximates the target Hamiltonian \({\hat{H}}_{T}\) given by Eq. (47) under the RWA, cf. Eq. (22). To assess the errors induced by the RWA, we also show the results for the ideal AQEC protocol, where the error correction Hamiltonian in the quantum master equation (46) is set exactly as the target Hamiltonian (47), i.e., \({\hat{H}}_{EC}={\hat{H}}_{T}\), which establishes a lower boundary for the infidelity. Our results show that the AQEC performance with our driving Hamiltonian protection nearly overlaps that with the ideal Hamiltonian protection. The slowly increasing infidelity after each measurement arises from the multi-photon loss, which cannot be corrected by the four-fold symmetric cat code states42. Correcting multi-photon loss errors would require cat code states with higher-order discrete rotational symmetry in phase space38. Another error source arises from the fact that the cat error states \(| {\bar{0}}_{e}\left.\right\rangle\) and \(| {\bar{1}}_{e}\left.\right\rangle\) do not satisfy the Knill-Laflamme condition at the first sweet spot of the cat code states. This error can be mitigated by choosing cat codes at higher-order sweet spots (i.e., with larger α). However, since the decay rate is proportional to the average photon number, this comes at the cost of increasing the frequency of parity measurements.

Discussions

First, we discuss the experimental implementation of our proposal. In principle, the quantum lattice gate operations introduced in Eqs. (4) and (5) can be directly demonstrated with a cold atom in optical lattice potentials. Here, for practical scalable quantum computing, we propose to implement the quantum lattice gates in the superconducting circuit-QED architecture. One can implement the quantum lattice gate operations by introducing inductive couplings between the cavity and the JJ-based loops. As illustrated in Fig. 10, a superconducting cavity (left) couples to, via a mediated JJ-based circuit (middle) referred to as the “coupler”, the SQUID device (right) shunted by a large capacitance Cs with an external loop. The coupler induces a tunable effective mutual inductance M(Φx) between the cavity and the SQUID circuit98. The net fluxes Φx, Φext, and Φsq that are applied to the coupler, the external loop of the SQUID, and the SQUID itself, respectively, can be independently tuned via the coupler bias line, loop bias line, and SQUID bias line. In total, the effective dynamics of the superconducting cavity is described by the Hamiltonian \(\hat{{{{\mathcal{H}}}}}(t)=\hslash {\omega }_{0}{\hat{a}}^{{{\dagger}} }\hat{a}+{V}_{{{{\rm{SQUID}}}}}(\hat{\varphi })\) with \({\omega }_{0}=1/\sqrt{LC}\) the bare frequency of the cavity. Note that the charging energy of SQUID circuit can be suppressed by the large shunted capacitance Cs99, and the inductive energy of the SQUID circuit is given by

$${V}_{{{{\rm{SQUID}}}}}(\hat{\varphi })=-{E}_{J}\cos \left(\frac{{\Phi }_{{{{\rm{sq}}}}}}{{\Phi }_{0}}\pi \right)\cos \left[\frac{M({\Phi }_{x})}{L}\hat{\varphi }+\frac{{\Phi }_{{{{\rm{ext}}}}}}{{\Phi }_{0}}\right],$$
(56)

where EJ is the single JJ energy, Φ0 = 2π/2e is the magnetic flux quantum and \(\hat{\varphi }=\sqrt{\frac{\lambda }{2}}({\hat{a}}^{{{\dagger}} }+\hat{a})=2\pi \hat{\Phi }/{\Phi }_{0}\) with dimensionless Planck constant λ = 4e2ω0L/ is the cavity flux variable. The effective mutual inductance can be tuned, via the coupler bias, to be zero [M(Φx) = 0], arbitrarily large [M(Φx) → ], or even negative [M(Φx) < 0]98. By rapidly turning on/off the coupling between the cavity and the SQUID circuit via the mutual inductance or the Josephson energy of the SQUID, one can in principle implement the lattice gate \({{{\rm{XSL}}}}(\rho ,\gamma ,\delta )=\exp [i\gamma \cos (\rho \hat{\varphi }+\delta )],\) cf. Eq. (5), where the gate parameters γ, ρ and δ are tuned by the biased flux Φsq, Φx and Φext, respectively.

Fig. 10: Experimental implementation.
figure 10

A superconducting cavity (LC circuit, left) coupled to a SQUID loop (red circuit, right) shunted by a large capacitance Cs via a coupler (blue circuit, middle). M(Φx) is the effective mutual inductance between the cavity and SQUID circuit induced by the coupler. Φx, Φext and Φsq are the net fluxes applied to the coupler, SQUID external loop, and SQUID, respectively, and can be independently tuned via the coupler bias, loop bias, and SQUID bias lines.

However, compared to conventional SQUID circuits, the proposed architecture in Fig. 10 faces practical challenges in implementing our proposed lattice gates, particularly due to the need for precise control of the mutual inductance M(Φx) and the external flux Φext. Alternatively, one can implement the quantum lattice gates directly using a single SQUID circuit together with standard Gaussian operations. In this case, one can implement the quantum lattice gate \({{{\rm{XSL}}}}(1,\gamma ,\delta )={e}^{i\gamma \cos (\hat{\varphi }+\delta )}\) by turning the capacitive coupling between the cavity and the SQUID circuit. To implement a more general lattice gate \({{{\rm{XSL}}}}(\rho ,\gamma ,\delta )={e}^{i\gamma \cos (\rho \hat{\varphi }+\delta )}\), as defined in Eq. (5), one can apply additional squeezing operations, i.e.,

$${{{\rm{XSL}}}}(\rho ,\gamma ,\delta )={\hat{S}}^{{{\dagger}} }(\ln \rho ){{{\rm{XSL}}}}(1,\gamma ,\delta )\hat{S}(\ln \rho ).$$
(57)

Here, \(\hat{S}(z)\equiv {e}^{\frac{1}{2}({z}^{* }{\hat{a}}^{2}-z{\hat{a}}^{{{\dagger}} 2})}\) is the squeezing operator that satisfies \({\hat{S}}^{{{\dagger}} }(z)\hat{\varphi }\hat{S}(z)={e}^{-z}\hat{\varphi }\) for a real parameter \(z\in {\mathbb{R}}\). This squeezing operation is a standard Gaussian operation that has been demonstrated in superconducting circuits with a squeezed drive25,37,100.

Second, we discuss the time and scalability of the gate sequence. As mentioned above in the “Quantum lattice gate” subsection in Methods, compared to the popular SNAP gates and cubic phase gates, our proposed elementary PSL gate has fewer parameters and a shorter gate operation time. As shown in Fig. 1 and Eq. (5), one single PSL gate can be realized within one nanosecond for a characteristic cavity (oscillator) frequency ω0 = 2π × 5.0 GHz. However, due to the adiabatic ramping protocol adopted in this work, the total time of state preparation costs about 400 μs for the typical adiabatic preparation time tf = 5000 × 2π/ω0 and the discretization parameters NM = 20, cf. Fig. 4a. In comparison, state preparation protocols based on SNAP gates typically cost 3–5 μs per gate operation and require also a total execution time of about 300–500 μs26.

However, the slow execution time of our quantum lattice gates is not limited by the PSL gate operation itself, but rather by the adiabatic ramping protocol required to satisfy the adiabaticity condition in Floquet systems61, cf. Eq. (30). If we go beyond the adiabatic approximation, the total execution time can be significantly reduced. For example, by applying the efficient numerical optimization method developed previously57, which was originally developed to reduce the execution time of SNAP gates (from hundreds of microseconds to tens of microseconds), the code states discussed in this work can be prepared with less than one hundred PSL gates, reducing the total execution time from hundreds of microseconds to the order of tens of nanoseconds.

The gate sequence optimized numerically for the SNAP and cubic phase gates is tailored to a given set of parameters, while the analytical decomposition with PSL gates remains valid for a general class of parameter settings. Given that the adiabacity is satisfied, the scalability of the PSL gate sequence depends on the complexity of the target code state, and can be determined from the analytical expression of the driving field modulation over a single Floquet period, cf. Fig. 4a. In Fig. 11a, b, we plot and compare the charts of driving amplitude modulation A(kt) for the four-legged cat state \(| {\bar{0}}_{c}\left.\right\rangle\) at the first (α ≈ 1.538) and the second (α ≈ 2.345) sweet spots, cf. Eqs. (8) and (9). As the second-sweet-spot cat state is more complex than the first-sweet-spot one, cf. the two insets in Fig. 11, its driving chart exhibits a finer structure and spans a broader range of the wavenumber k. Therefore, more grids are needed for the discretized quantum lattice gates to prepare code states with the desired precision.

Fig. 11: Driving protocols for cat states at different sweet spots.
figure 11

Charts of the driving amplitude A(kt) for preparing the four-legged cat state \(| {\bar{0}}_{c}\left.\right\rangle\), cf. Eqs. (8) and (9), at a the first and b the second sweet spots. The insets in each chart show the Wigner functions of the corresponding target code states.

Third, we discuss the robustness of our proposal against various gate noises. Compared to other universal gate sets, the advantage of our proposed quantum lattice gates is the analytical ability to control the gate errors with deterministic parameters, albeit with a slow total execution time. For SNAP and cubic phase gates, the gate sequence can only be optimized by numerical methods such as gradient ascent pulse engineering (GRAPE)101. In contrast, the mathematical form of our proposed PSL gates enables the analytical decomposition of quantum circuits using the NcFT technique, and thus offers more flexibility in designing target Hamiltonians with desired properties. In fact, the gaps introduced in the target Hamiltonians, cf. Eqs. (13), (14), and (16) make the code state preparation robust against gate noises arising from the driving amplitude β(t) and frequency Ω(t) during the adiabatic ramp, cf. Eq. (29). To show this, we introduce the noisy driving amplitude and frequency by

$$\beta (t)\to \beta (t)+{\beta }_{{{{\rm{noise}}}}}\xi (t),\quad \Omega (t)\to \Omega (t)+{\Omega }_{{{{\rm{noise}}}}}\xi (t),$$

where ξ(t) represents the standard Gaussian noise with average 〈ξ(t)〉 = 0 and correlation \(\langle \xi (t)\xi ({t}^{{\prime} })\rangle =\delta (t-{t}^{{\prime} })\).

In Fig. 12a, we plot the infidelity 1 − F(tf) of the prepared binomial code state \(| {\bar{0}}_{b}\left.\right\rangle\) defined in Eq. (10) at the final adiabatic ramp time tf = 2000 × 2π/ω0 as a function of the driving noise strength βnoise/βf for different driving amplitudes βf = 0.025ω0 and βf = 0.05ω0. Thanks to the gap (βΔ) in the spectrum of the engineered Floquet Hamiltonian (\({\hat{H}}_{F}=\beta {\hat{H}}_{T}\)), the infidelity is robust against the increase of the driving amplitude noise. In Fig. 12b, we plot the infidelity 1 − F(tf) of the prepared binomial code state \(| {\bar{0}}_{b}\left.\right\rangle\) as a function of the frequency noise Ωnoise/ω0 with different detunings Ω(0) − ω0 at the beginning of the adiabatic ramp, cf. Eq. (29). Compared to the driving amplitude noise, the infidelity of the prepared state is more sensitive to the frequency noise. Considering the flux-induced energy fluctuations in superconducting circuits are mainly caused by the low-frequency 1/f noise102,103,104, the resulting infidelity remains below 10−2 for a frequency noise of Ωnoise ~ 10−6ω0105. The frequency noise can be further suppressed by increasing the initial detuning Ω(0) − ω0, as indicated by the red infidelity curve in Fig. 12b.

Fig. 12: Robustness against gate noises.
figure 12

a Infidelity 1 − F(tf) of the prepared binomial code state \(| {\bar{0}}_{b}\left.\right\rangle\), cf. Eq. (10), as a function of the amplitude noise strength βnoise/βf, for different final driving amplitudes: βf = 0.025ω0 (blue) and βf = 0.05ω0 (red). The initial driving frequency is detuned as Ω(0) = ω0/(1 + π × 10−3). b Infidelity 1 − F(tf) of the prepared binomial code state \(| {\bar{0}}_{b}\left.\right\rangle\) as a function of the frequency noise strength Ωnoise/ω0, for different initial driving frequencies: Ω(0) = ω0/(1 + π × 10−3) (blue) and Ω(0) = ω0/(1 + π × 10−2) (red). The final driving amplitude is set to be βf = 0.025ω0. c Infidelity 1 − F(tf) of the prepared binomial code state \(| {\bar{0}}_{b}\left.\right\rangle\) as a function of the phase-space lattice (PSL) gate phase noise strength ϵα/2π, with the driving frequency Ω(0) = ω0/(1 + π × 10−3) and the final driving amplitude βf = 0.025ω0. In all figures, the final time is fixed at tf = 2000 × 2π/ω0. Each data point represents an average over 40 noise samples, with the error bars representing one standard deviation.

In addition to quantum error correction, fault-tolerant quantum computing imposes an additional stringent requirement on gate operations: errors should not be amplified by the gate sequence in damaging ways38. In fact, the ladder operators of the cavity and the associated unitary error operators can be expanded in terms of small displacement operators18,38. When a displacement error operator \(\hat{D}(\alpha )={e}^{\alpha {\hat{a}}^{{{\dagger}} }-{\alpha }^{* }\hat{a}}\) passes through a PSL gate, the gate exhibits a small phase shift, i.e., \({{{\rm{PSL}}}}(\zeta ,\sigma ;\gamma ,\delta )\hat{D}(\alpha )=\hat{D}(\alpha ){{{\rm{PSL}}}}(\zeta ,\sigma ;\gamma ,\delta +\Delta \delta )\) with \(\Delta \delta =\sqrt{2\lambda }[\zeta {{{\rm{Re}}}}(\alpha )+\sigma {{{\rm{Im}}}}(\alpha )]\). We model such gate phase fluctuation Δδ with a standard Gaussian random number multiplied by a strength factor ϵα, and calculate the fidelity of the prepared state with such noise process. As shown in Fig. 12c, as long as the phase noise is smaller than a critical value (ϵαπ × 10−2), such errors do not accumulate in a way that severely impacts the gate sequence, even in deep quantum circuits based on quantum lattice gates.

Finally, we discuss the possibility of extending our proposal to multi-mode versions. As already mentioned in the “Universal gate sets for bosonic modes” subsection in Methods, a universal set of single-mode CV gates, together with a linear two-mode interaction such as the CSUM gate \({e}^{-i{\hat{x}}_{1}{\hat{p}}_{2}}\), suffice to enact universal quantum computing with multiple CV modes23,29. Here, we outline a general framework how to implement universal multi-mode CV quantum computing using our proposed PSL gates. As an example, we consider a two-mode unitary transformation generated by a target Hamiltonian \({\hat{H}}_{T}\equiv H({\hat{x}}_{1},{\hat{p}}_{1},{\hat{x}}_{2},{\hat{p}}_{2})\), and decompose it into a sequence of single-mode PSL gates and two-mode CSUM gates. To this end, we discretize the target Hamiltonian as follows:

$$\begin{array}{l}{\hat{H}}_{T}({\hat{x}}_{1},{\hat{p}}_{1},{\hat{x}}_{2},{\hat{p}}_{2})\approx {\sum }_{n=0}^{N}{\gamma }_{n}\cos \left[{\eta }_{1,n}{\hat{x}}_{1}+{\eta }_{2,n}{\hat{x}}_{2}+{\sigma }_{1,n}{\hat{p}}_{1}+{\sigma }_{2,n}{\hat{p}}_{2}+{\delta }_{n}\right].\end{array}$$
(58)

Here, the coefficients γnη1(2),nσ1(2),n, and δn can be obtained by extending the NcFT method for a single-mode plane-wave operator \({e}^{i({k}_{x}\hat{x}+{k}_{p}\hat{p})}\), cf. Eq. (23), to a two-mode correspondence \({e}^{i({k}_{1,x}{\hat{x}}_{1}+{k}_{1,p}{\hat{p}}_{1}+{k}_{2,x}{\hat{x}}_{2}+{k}_{2,p}{\hat{p}}_{2})}\). Note that, in the single-mode case, the success of our NcFT technique relies on a one-to-one correspondence between the periodic time and the phase degree of freedom of the bosonic mode. Extending this to the two-mode case introduces a complication: multiple phase degrees of freedom must be encoded using a single time parameter. This issue can be addressed by imposing proper constraints on the phases, which we will discuss in detail in a separate work.

Given the coefficients, the elementary Hamiltonian action from in Eq. (58), i.e.,

$${{{{\rm{PSL}}}}}_{1,2}\equiv {e}^{i{\gamma }_{n}\cos \left({\eta }_{1,n}{\hat{x}}_{1}+{\eta }_{2,n}{\hat{x}}_{2}+{\sigma }_{1,n}{\hat{p}}_{1}+{\sigma }_{2,n}{\hat{p}}_{2}+{\delta }_{n}\right)},$$

can be generated by a single-mode lattice gate operation \({{{{\rm{XSL}}}}}_{2}\equiv {e}^{i{\gamma }_{n}\cos {\hat{x}}_{2}}\) together with a two-mode CSUM gate. This can be achieved as follows. First, treating mode 1 as the control and mode 2 as the target, the CSUM gate enacts \({\hat{x}}_{2}\to {\hat{x}}_{2}+{\hat{x}}_{1}\), and thus transforms the gate XSL2 into \({{{{\rm{XSL}}}}}_{2,1}\equiv {e}^{i{\gamma }_{n}\cos ({\hat{x}}_{2}+{\hat{x}}_{1})}\). Then, by applying appropriate single-mode Gaussian gates (phase rotation, squeezing, and displacement gates) on each mode individually, one can generate the PSL1,2 gate defined above. Additionally, the two-mode CSUM gate can be replaced by a CZ gate \({e}^{i{\hat{x}}_{1}{\hat{x}}_{2}}\) if followed by a phase rotation gate that maps \({\hat{p}}_{2}\to {\hat{x}}_{2}\) on mode 2.

Conclusions

In this work, we have introduced a universal quantum lattice gate for bosonic quantum computation, which can be implemented in a JJ-based superconducting circuit architecture. Compared to other popular gates, such as the SNAP gate26 and the cubic phase gate23,28, which typically overlook higher-order nonlinearities in the JJ potential, our quantum lattice gates leverage the full nonlinearity of the JJ circuit as a powerful quantum resource. We have developed a comprehensive analytical framework to decompose a given bosonic code engineering process into sequences of quantum lattice gates through Floquet Hamiltonian engineering. Our proposal also offers a direct method for constructing target Hamiltonians aligned precisely with the desired target codes. We have demonstrated the versatility of our method across various applications, including code state preparation, code space embedding, and code space transformation, specifically with cat and binomial code states. We have also studied the autonomous quantum error correction of cat code states against photon loss errors with our driving protocol. The results not only advance the toolkit for bosonic quantum computation but also pave the way for more robust and scalable implementations in superconducting quantum circuits.

In general, our lattice gate decomposition for code state engineering is based on the adiabatic ramp protocol61, which limits the speedup of the total execution time. Nevertheless, the strength of our method lies in its analytical capability, which allows us to identify errors and suppress them. Specifically, the non-adiabatic errors and gate noises can be suppressed by enlarging the gaps in the designed target Hamiltonians; the non-RWA errors arising from higher-order terms in the Floquet-Magnus expansion can be mitigated via perturbative inverse Floquet Hamiltonian engineering96; the Trotterization errors from the gate discretization can be reduced by employing higher-order Suzuki-Trotter decompositions106.

Recently, the contribution of higher-order Josephson harmonics in tunnel junctions has been observed107. To account this, one can generalize the quantum lattice gate introduced in Eq. (5) as \({{{\rm{XSL}}}}(\rho ,\gamma ,\delta )\equiv \exp [i\gamma {\sum }_{m\ge 1}{E}_{Jm}\cos (m\rho \hat{x}+m\delta )]\), where ρ and δ are the squeezing and displacement parameters, respectively. The gate sequence for this modified quantum lattice gate can be identified via standard numerical optimization methods. In the future, we plan to explore methods to accelerate state engineering by incorporating counter-adiabatic terms, following Berry’s solution108. These terms can be engineered through additional driving fields using our technique of arbitrary Floquet Hamiltonian engineering63, or by directly optimizing the sequence and parameters of the quantum lattice gates with various numerical methods57,101,109. In short, it would be promising to enhance the efficiency and speed of our state engineering protocols, pushing the boundaries of what can be achieved in bosonic quantum computation.