Generating arbitrary superpositions of nonclassical quantum harmonic oscillator states
S. Saner O. Băzăvan D. J. Webb G. Araneda D. M. Lucas C. J. Ballance R. Srinivas








§ ABSTRACT 
Full coherent control and generation of superpositions of the quantum harmonic oscillator are not only of fundamental interest but are crucial for applications in quantum simulations, quantum-enhanced metrology and continuous-variable quantum computation. The extension of such superpositions to nonclassical states increases their power as a resource for such applications. Here, we create arbitrary superpositions of nonclassical and non-Gaussian states of a quantum harmonic oscillator using the motion of a trapped ion coupled to its internal spin states. We interleave spin-dependent nonlinear bosonic interactions and mid-circuit measurements of the spin that preserve the coherence of the oscillator. These techniques enable the creation of superpositions between squeezed, trisqueezed, and quadsqueezed states, which have never been demonstrated before, with independent control over the complex-valued squeezing parameter and the probability amplitude of each constituent, as well as their spatial separation. We directly observe the nonclassical nature of these states in the form of Wigner negativity following a full state reconstruction. Our methods apply to any system where a quantum harmonic oscillator is coupled to a spin.

Generating arbitrary superpositions of nonclassical quantum harmonic oscillator states
    S. Saner, O. Băzăvan, D. J. Webb, G. Araneda, D. M. Lucas, C. J. Ballance, R. Srinivas 

Department of Physics, University of Oxford, Clarendon Laboratory, 

Parks Road, Oxford OX1 3PU, United Kingdom 
Email: sebastian.saner@physics.ox.ac.uk
    September 9, 2024
=====================================================================================================================================================================================================================================================






One of the most fundamental and remarkable aspects of quantum mechanics is superposition, where quantum objects can exist in multiple states simultaneously. This phenomenon was famously described by Schrödinger using the example of a cat that is both dead and alive <cit.>. Superpositions are not just a curiosity of nature but are central to any application of quantum mechanics. In sensing <cit.>, it is the superposition of two states that respond differently to a given signal, such as a shift in frequency, that enables their use in atomic clocks <cit.>. 
In quantum computing, the ability to be in multiple states at once increases the amount of information that can be stored and processed simultaneously and is the basis of quantum advantage <cit.>. 

Superpositions of two-level systems can be controlled to almost arbitrary precision in many systems <cit.>. However, they ultimately only have two degrees of freedom: the relative amplitude and phase of their complex probabilities, which limits the amount of quantum information and how robustly it can be stored. Creating superpositions of a quantum harmonic oscillator, which is, in contrast, infinite-dimensional, has already led to a plethora of resource-efficient error correction protocols <cit.>, metrological applications <cit.>, and quantum simulations <cit.>. Nonetheless, constituents of those superpositions have been limited to two <cit.> or three <cit.> Fock states, classical states such as coherent states leading to cat states <cit.>, or Gaussian states such as displaced squeezed states leading to Gottesman-Kitaev-Preskill (GKP) states <cit.>. 



Instead, superpositions between arbitrary squeezed states <cit.>, or between non-Gaussian states such as trisqueezed states, exhibit increased Wigner negativities <cit.> and discrete rotational symmetries which can be a resource for new sensing <cit.> and error correction schemes <cit.>.
However, the lack of any strong, native interactions that create such superpositions has precluded any experimental realisation.



Here, we use a hybrid oscillator-spin system to create such superposition states <cit.> of the harmonic oscillator; we use our recently demonstrated strong spin-dependent nonlinear interactions of the oscillator <cit.> together with mid-circuit measurements of the spin. We create an arbitrary combination of generalised squeezed states <cit.> such as squeezed, trisqueezed, and quadsqueezed harmonic oscillator states.
These techniques not only enable arbitrary control of the complex probability amplitudes of the constituents forming the superposition, but also over the interactions that generate them. We can create constituents from interactions of the same type with different strengths and phases, or from entirely different nonlinear interactions. 


The hybrid system is central to our demonstration. 
Applying the spin-dependent nonlinear interaction to a superposition of spin states leads to
an entangled state of the spin and the nonclassical motional states.
A mid-circuit measurement of the spin subsequently disentangles the spin from the motion, leaving the harmonic oscillator in a superposition state. 
This method enables the repeated application of the nonlinear interaction within a single sequence, as the interaction is unitary and can act on any starting state without destroying it.

We create the spin-dependent nonlinear interaction

    Ĥ_k =ħΩ_k/2σ̂_z (â^k e^-iϕ + (â^†)^k e^iϕ),

by combining two noncommuting spin-dependent forces (SDF) <cit.> (see Supplement <cit.>). 
The spin Pauli operator is defined as σ̂_z = |1_s⟩⟨1_s| - |0_s⟩⟨0_s|, while â^† (â) describes the creation (annihilation) operator of a quantum harmonic oscillator. The order of the generalised squeezing interaction <cit.> is k; setting k=2, 3, 4 generates the squeezing, trisqueezing, and quadsqueezing interactions respectively. Squeezing interactions have been used extensively in metrology <cit.>, while the higher-order trisqueezing <cit.> and quadsqueezing <cit.> create non-Gaussian states of the harmonic oscillator <cit.>.
The coupling strength of the interaction is Ω_k, and ϕ is the phase relative to the oscillator, which also defines the generalised squeezing axis.
When applying this interaction to the initial state of the hybrid system |0_osc,0_s⟩, where |0_osc⟩ represents the harmonic oscillator in its ground state and |0_s⟩ denotes the spin configuration prepared in an eigenstate of σ̂_z, we obtain a generalised squeezed state |ζ_k, 0_s⟩:=exp(-iĤ_k t/ħ)|0_osc, 0_s⟩ with squeezing magnitude ζ_k = Ω_k e^iϕ t.

A superposition between two generalised squeezed states |ζ_k⟩ and |-ζ_k⟩, with orthogonal squeezing axes, can be created by applying the sequence shown in Fig. <ref>A. Initially, a R_y(π/2) rotation brings the spin part of the ion into a superposition state |+_s⟩= (|0_s⟩+|1_s⟩)/√(2) before applying the σ̂_z conditioned nonlinear coupling. After the nonlinear coupling and another R_y(π/2) rotation have been applied, the resulting entangled state of the hybrid system takes the form


    |ψ⟩ = N_k,+|ψ_k,+⟩|1_s⟩ + N_k,-|ψ_k,-⟩|0_s⟩,

where N_k,± denote the probability amplitudes of the even (|ψ_k,+⟩) and odd (|ψ_k,-⟩) superpositions of the oscillator, which arise directly from the normalisation coefficient of |ψ_k,±⟩ = 1/(2N_k,±)· (|ζ_k⟩±|-ζ_k⟩). The normalisation coefficients depend on the overlap between |ζ_k⟩ and |-ζ_k⟩ and, consequently, on the magnitude and order of the applied nonlinear interaction.

To project into the desired oscillator superposition state, we perform a mid-circuit measurement detecting |0_s⟩ or |1_s⟩, which occurs with a probability of |N_k,±|^2. This measurement collapses the state into the corresponding superposition of the oscillator, disentangled from the spin state. After the collapse, the created superposition can be employed in further experiments, or characterised through full tomography. 
Experimentally, we can select the outcomes with the desired superposition by postselecting based on the mid-circuit measurement result.


We experimentally demonstrate the creation of these superpositions on a trapped ^88Sr^+ ion in a 3D Paul trap <cit.>. The harmonic oscillator is formed by the axial motion of the confined ion with frequency 1.2. In practice, we cannot initialise a pure vacuum state of the oscillator |0_osc⟩ but rather a thermal state near its ground state with average occupation n̅=0.1. The coherence time of this oscillator is limited by its heating rate ṅ̅̇=300  quanta/s. 
The internal structure forms the spin system, with the states |0_s⟩≡|5S_1/2, m_j = -1/2⟩, and |1_s⟩≡|4D_5/2, m_j = -3/2⟩. The electronic state |2_s⟩≡|5S_1/2, m_j = 1/2⟩ is used as an auxiliary state.
 
The |0_s⟩↔|1_s⟩ quadrupole transition is driven by a 674 nm laser; the two noncommuting spin-dependent forces required to generate the nonlinear interaction are created by two 674 nm beams. The |2_s⟩ state is isolated from the interaction, and an RF antenna drives the |0_s⟩↔|2_s⟩ transition.

The projective measurement is achieved by driving the 5S_1/2↔ 5P_1/2 transition 
which distinguishes bright states (|0_s⟩ or |2_s⟩), which scatter photons, from the dark state (|1_s⟩), which does not. This measurement is only non-destructive to the motional state in the case where |1_s⟩ is detected, and hence no photons are actually scattered. Thus, we design our pulse sequences such that the desired oscillator state is heralded upon |1_s⟩, i.e., the desired harmonic superposition is correlated with |1_s⟩ just before the measurement. For example, for the state described in Eq. (<ref>), we would select the even superposition |ψ_k,+⟩.
The projective measurement is used mid-circuit for selecting the desired superposition as well as for the final detection step. In both cases, we collect photons for 200, yielding a combined dark state preparation and measurement fidelity of 0.993 <cit.>.


We create both even and odd superpositions between the state |ζ_k⟩ and |-ζ_k⟩ as illustrated in Eq. (<ref>) and Fig. <ref>. By measuring the qubit immediately after this creation step, we obtain the probability of collapsing in |0_s⟩ and |1_s⟩ respectively and, hence, the probability of generating the even and odd superposition. For example, for k=2 in Eq. (<ref>), (<ref>) we obtain the superposition between two oscillator states squeezed about orthogonal axes. The corresponding probabilities for the even (|ψ_2,+⟩) and odd (|ψ_2,-⟩) states are given by |N_2,±|^2 = (2± 2/√(cosh(2|ζ_2|)))/4, which depend on the magnitude of the squeezing parameter |ζ_2| (see Supplement  <cit.>).
In the limit of small squeezing parameter |ζ_2|, the generation of the odd superposition vanishes, which indicates that the overlap between |ζ_k⟩ and |-ζ_k⟩ tends to 1 as the oscillator state is unchanged from the initial state.
In the limit |ζ_2| →∞, the odd and even state have equal generation probability of 1/2 which indicates that |ζ_2⟩ and |-ζ_2⟩ are orthogonal and their overlap is 0. We confirm this behaviour experimentally; the data in Fig. <ref>A agree well with the theoretical prediction from an independent estimate of |ζ_2|. 

Extending this protocol to higher-order nonlinear interactions, we reconstruct the Wigner quasi-probability distribution of even and odd superpositions for squeezing (k=2), trisqueezing (k=3), quadsqueezing (k=4) (Fig. <ref> B-G). To do so, we measure the characteristic function χ(β) = ⟨D̂(β)|_⟩osc of the oscillator state and infer the Wigner distribution via a Fourier transform <cit.>. The characteristic function χ(β) is measured via the spin using a spin-dependent displacement D̂(σ̂_yβ/2) as shown in the tomography step in Fig. <ref>A. By varying the complex-valued displacement parameter β, we can sample different values of χ(β) (see Supplement <cit.>).   
 We find good agreement between the experimental results and numerical simulations using independently measured parameters (Fig. <ref> H-M). In the Supplement <cit.>, we show the experimental data and numeric predictions of |N_k,±|^2 as a function of |ζ_k| for k ∈{3,4}.

The superpositions created, especially the odd ones, have a large amount of Wigner negativity. In continuous variable quantum computation, large Wigner negativities are a requirement for quantum states that yield an advantage over classical computation <cit.> and hence highly sought after. 
We quantify the resourcefulness of these states by evaluating their Wigner logarithmic negativity (WLN) <cit.>, neglecting decoherence effects, and benchmark them against Fock and cat states. 
We find that for a fixed average Fock state occupation, the superpositions created in this work exhibit larger WLN than Fock or cat states <cit.> (see Supplement <cit.>).

Another important feature is that the non-vanishing Fock state occupations of the superpositions are spaced by 2k, where k is the order of the interaction (see Supplement <cit.>). As k increases, so does the spacing, and the states have increased rotational symmetries as seen in Fig. <ref>. These rotational symmetries can be used to encode more robust logical qubits <cit.>. For example, a single phonon loss in a (two-component) cat-qubit encoding results in an irreparable bit-flip error <cit.>. However, when the same loss channel acts on a qubit encoding formed by a superposition of squeezed states, it projects outside the code space, making the error correctable. 

Thus far, the superpositions have only involved constituents generated from the same nonlinear interactions. To create arbitrary superpositions with nonclassical constituents in the most general form 

    |ψ⟩∝ a |ζ_k⟩ + b |ζ_k'⟩,

we require individual control of a,b,ζ_k, ζ_k'∈ℂ. To create this state, we extend the Hilbert space of the spin to a qutrit (three-level system), {|0_s⟩,|1_s⟩,|2_s⟩} (see Fig. <ref>A, B).
We create the constituent states via the same interaction described in (<ref>), which only couples to |0_s⟩ when applied to the subspace {|0_s⟩, |2_s⟩}. Thus, we utilize this feature to temporarily hide one of the constituents |ζ_k,2_s⟩ while creating the other using a separate interaction <cit.>. 

The pulse sequence for creating the arbitrary superpositions is described in Fig. <ref>A; for simplicity, we first focus on the case a=b.
After the first nonlinear interaction is applied, the hybrid system is in a state |0_osc,2_s⟩ + |ζ_k, 0_s⟩. A subsequent R_y(π)_02 rotation, where the indices indicate the two spin sublevels |0_s⟩↔|2_s⟩ coupled, then swaps the spin states, hiding the state |ζ_k, 2_s⟩ from the effect of the second nonlinear interaction. The second nonlinear interaction then creates the state |ζ_k,2_s⟩ + |ζ_k'', 0_s⟩. One final R_y(π/2)_02 pulse followed by the shelving R_y(π)_01 pulse creates a state
described by Eq. (<ref>), with |ψ_k,±⟩∝|ζ_k⟩±|ζ_k''⟩ where ζ_k'' is no longer -ζ_k.

We then characterise the state using the same techniques as before. As mentioned previously, the unitarity of the nonlinear interaction is essential to prevent the destruction of the intermediate state by subsequent applications of the nonlinear interaction.
Splitting the sequence into two parts allows us to change the second nonlinear interaction relative to the first one <cit.>.

In Fig. <ref>C, we show that we can create a superposition of squeezed states with variable orientation of the respective squeezing axis, i.e. ζ_2'=e^i  2ϕζ_2, by adjusting the phase ϕ of the nonlinear interaction in the second part of the sequence relative to the first part. Hence, the created state takes the form |ζ_2⟩+|e^i  2ϕζ_2⟩. We measure the value of the characteristic function at a fixed magnitude |β| but vary its complex phase (β). At ϕ=π/2 we create a superposition of states 
squeezed about orthogonal axes (equal superposition of squeezed states shown in Fig. <ref>B, H), which gives rise to four equally spaced maxima in χ when scanning (β). By scanning the degree of freedom ϕ, the two orthogonal squeezing axes gradually merge into a single one at ϕ={0, π} (regular squeezed state), giving rise to two equally spaced maxima. This behaviour results in the sawtooth pattern displayed in Fig. <ref>C.

In Fig. <ref>D, we show that we can control the magnitude of squeezing of each constituent independently <cit.>.
We keep the squeezing parameter ζ_2 of the first constituent fixed while adjusting the magnitude of the second constituent by a factor c ∈ℝ which results in a state |ζ_2⟩+|c·ζ_2⟩ i.e. ζ_2'=cζ_2. For c=0, this corresponds to a superposition of a squeezed state with the initial ground state, whereas for c=-1, this corresponds to the equal superposition shown in Fig. <ref>B, H.
We measure the variance of the state along the squeezing axis of the first constituent and an axis orthogonal to it (i.e., two orthogonal quadratures). We compute the variances using the reconstructed Wigner functions for different settings of c. As shown in Fig. <ref>D, the variance closely follows the simulated prediction.

The variance of each constituent does not vary independently. In fact, altering the magnitude of the second constituent affects the variance of the resulting state in both quadratures: the quadrature dominated by the second constituent (blue) and the quadrature dominated by the first constituent (green), even though the first constituent's squeezing magnitude remains fixed. This behaviour is a direct consequence of the interference between both constituents.

To demonstrate that we can change the order k of the nonlinear interaction for each constituent, we combine a squeezed state (k=2) as the first constituent with a trisqueezed state (k=3) as the second constituent. Fig. <ref>E shows the resulting Wigner function of this superposition state.

While we do not include the data in the main text, we discuss in the Supplement <cit.> that the sequence can be amended to vary the relative amplitude and phase of the probability coefficients a and b. By replacing the first R_y(π/2) rotation by R_y(θ) to vary the probability amplitude ratio between the two constituents cos(θ/2)|ζ_k⟩ + sin(θ/2)|-ζ_k⟩. Further, it is also possible to replace the R_y(π) rotation with a rotation about an axis σ̂_γ = cos(γ)σ̂_x + sin(γ)σ̂_y to introduce a complex phase factor between the probability amplitude of the two constituents |ζ_k⟩ + |-ζ_k⟩ e^i2γ. Thus, we demonstrate complete control over a, b, ζ_k and ζ_k'' in Eq. (<ref>). 

Finally, we show that the sequence depicted in Fig. <ref>A can be extended further to construct spatially separated superposition states. Aside from a foundational quantum optics interest, such spatially-separated superpositions are essential in various applications, including continuous variable quantum error correction codes (e.g., cat <cit.> or GKP <cit.> states), and metrology <cit.>. We create these states using an additional spin-dependent displacement followed by a mid-circuit measurement as shown in Fig. <ref>. The second mid-circuit detection disentangles the spin from the harmonic oscillator, projecting the latter into a cat-like state, where the two spatially displaced states are themselves nonclassical and non-Gaussian. We use this protocol to create a spatially separated superposition of squeezed states 

    |ψ⟩∝ D(α)|ψ_2,+⟩+D(-α)|ψ_2,+⟩,

where |ψ_2,+⟩∝|ζ_2⟩+ |-ζ_2⟩.
Thus, we show that we can create an even broader class of oscillator superpositions by concatenating additional spin-dependent interactions and mid-circuit measurements.

Our work demonstrates the versatility of hybrid oscillator-spin systems in generating superpositions of oscillator states. By integrating nonlinear spin-dependent interactions with mid-circuit measurements, we have established arbitrary control of the phase, amplitude, and interaction type for the constituents of the superpositions. Using these techniques, we have experimentally demonstrated many nonclassical states that have until now only been explored theoretically <cit.>. Expanding complex harmonic oscillator states from cat <cit.> and GKP <cit.> states to superpositions of squeezed or superpositions of non-Gaussian states is not only of foundational interest but opens new avenues for error correction, sensing protocols, and continuous variable or hybrid oscillator-spin quantum computing <cit.>. For example, these superposition states have increased rotational symmetries, which enable robust error correction <cit.>.  Further, their large Wigner negativities indicate that they generally cannot be efficiently simulated classically <cit.>. Their utility applies not only to quantum computation but also to quantum-enhanced metrology. For example, the squeezed superposition state is suitable as a displacement sensor <cit.>, which, in the case of trapped ions, can be used for sensing small electric fields <cit.>.
The demonstrated techniques can be extended to multiple oscillators by coupling to additional motional modes of the ion, for example, to create two-mode squeezed superpositions <cit.> or in general to increase the size of the hybrid quantum system. Aside from trapped ions, the tools developed here can be applied to any physical system with a quantum harmonic oscillator coupled to a spin such as superconducting circuits <cit.>, nanoparticles <cit.>, or atoms coupled to a cavity <cit.> or in optical tweezers <cit.>. For systems with more massive oscillators, these superpositions could be used to test not only the boundaries of the classical and quantum world <cit.>, but how quantum physics interacts with gravity <cit.>. 





§ ACKNOWLEDGEMENTS

We would like to thank Alejandro Bermudez, Alexander Lvovsky, Andrew Daley, Joshua Combes and his research group, Mattia Walschaers, and Scott Parkins for very insightful discussions and comments on the manuscript.
This work was supported by the US Army Research Office (W911NF-20-1-0038) and the UK EPSRC Hub in Quantum Computing and Simulation (EP/T001062/1). GA acknowledges support from Wolfson College, Oxford. CJB acknowledges support from a UKRI FL Fellowship. RS acknowledges funding from the EPSRC Fellowship EP/W028026/1 and Balliol College, Oxford. 



§ AUTHOR CONTRIBUTIONS

SS and OB led the experiments and analysed the results
with assistance from DJW, GA, and RS; SS, OB, DJW, and GA updated and maintained the experimental apparatus; SS, OB, and RS conceived the experiments and wrote the manuscript with input from all authors; RS supervised the work with support from DML and CJB; DML and CJB secured funding for the work.



§ COMPETING INTERESTS

RS is partially employed by Oxford Ionics Ltd. CJB is a director of Oxford Ionics Ltd. All other authors declare no competing interests.



§ DATA AVAILABILITY

Source data for all plots and analysis code that support the plots are available from the corresponding authors upon reasonable request.

apsrev4-2






















Supplemental Material for:

Generating arbitrary superpositions of nonclassical quantum harmonic oscillator states







§ THEORY



 §.§ Generating the nonlinear spin-dependent interaction

While the method employed to generate the nonlinear interaction is described in detail in Refs. . We summarise the relevant details below.
We combine two spin-dependent forces (SDF) with non-commuting spin-conditionings [σ̂_α, σ̂_α'] ≠ 0 and ϕ the oscillator phase of the second SDF which results in a Hamiltonian

    Ĥ   = ħΩ_α/2σ̂_α (â  e^-iΔ t+ â^† e^iΔ t)
       +ħΩ_α'/2σ̂_α' (â  e^-i (mΔ  t+ϕ) + â^† e^i(mΔ  t+ ϕ)).

The two forces are detuned from resonance with the oscillator by Δ and mΔ, respectively. By choosing k = 1-m, we can resonantly drive effective generalised squeezing interactions

    Ĥ_k =ħΩ_k/2σ̂_β (â^k e^-iϕ + (â^†)^k e^iϕ),

which for σ̂_β = σ̂_z corresponds to Eq. (<ref>) in the main text, where

    σ̂_β   ∝
        [σ̂_α, σ̂_α']   if  k  mod 2=0
    σ̂_α'   otherwise.

While we can choose our spin bases such that σ̂_β = σ̂_z, and do so for the majority of the data,
it was convenient for Fig. <ref>DE to choose σ̂_β = σ̂_x which maximises the strength and avoids two single qubit rotations. We have indicated the native spin-conditioning for every sequence in Sec. <ref>. 
The coupling strength Ω_k 

    Ω_2,3,4 = {Ω_α'Ω_α/Δ, Ω_α'Ω_α^2/2 Δ^2,Ω_α'Ω_α^3/8 Δ^3}·sin(θ_α,α'),

which is proportional to sin(θ_α, α') where θ_α,α' denotes the angle between the spin bases α and α' and is maximised if the two spin bases are orthogonal. We utilise this angle to adjust the squeezing magnitude for the constituent |c·ζ_2⟩ in the superposition of Fig. <ref>D, while keeping the optical power and the duration of the interaction fixed.
For the quantum circuits presented in the main text, we adjust Δ, Ω_α, Ω_α', and σ̂_β. The settings are tabulated in Sec. <ref>.  



 §.§ Normalisation coefficient

We require the squeezed state superposition to be normalised i.e. ⟨ψ_± |ψ_±|=⟩ 1, and find a normalisation coefficient

    |N_2,±|^2    = 1/4(⟨ζ_2|±⟨ζ_2|)(|ζ_2⟩±|ζ_2⟩)
       =1/4⟨ζ_2|ζ_2|+⟩⟨-ζ_2|-ζ_2|±⟩2 Re(⟨-ζ_2|ζ_2|)⟩
       = (2 ±2/√(cosh(2|ζ_2|)))/4,

where we assume that each squeezed state is normalised. The inner product ⟨-ζ_2|ζ_2|$⟩ can be determined from the Fock state expression of a vacuum squeezed state <cit.>:

    |ζ_2⟩   = 1/√(cosh(|ζ_2|))∑_n=0^∞ (- e^i (ζ_2)tanh(|ζ_2|))^n √((2n)!)/2^n n!|2n⟩,
       ⟨-ζ_2|ζ_2|=⟩1/cosh(|ζ_2|)∑_n=0^∞ e^i π ntanh(|ζ_2|)^2n(2n)!/2^2n (n!)^2
       ⟨-ζ_2|ζ_2|⟩= 1/cosh(|ζ_2|)1/√(tanh(|ζ_2|)^2 +1 ) = 1/√(cosh(2|ζ_2|)).


For the superpositions involving non-Gaussian (i.e. trisqueezed or quadsqueezed) constituents, it is unclear if closed-form solutions exist <cit.>.



 §.§ Native gate sequences


In this section, we give a more detailed description of the pulse sequence used to generate the superposition states. For simplicity, we omitted the spin-echoR(π)pulses described below from the main text and Fig. <ref>.

All nonlinear interactionsÛ_NLexcept for trisqueezing in Fig. <ref> are conditioned onσ̂_z= |1_s⟩⟨1_s| - |0_s⟩⟨0_s|.
For squeezing and quadsqueezing, we setσ̂_α = σ_xandσ_α'= σ̂_ywhich results inσ̂_β= σ̂_z. 
For the trisqueezing superposition shown in Fig. <ref>, we keep this setting which results in the spin basis of the nonlinear interactionσ̂_β= σ̂_x.
For trisqueezing shown in Fig. <ref>, we change the trisqueezing spin basis compared to Fig. <ref>. We set the spin conditioning of the SDF toσ̂_α'=σ̂_zsuch that the resulting nonlinear interaction is conditioned onσ̂_β=σ̂_z.


    |0_osc,0_s⟩   R_y(π/2)→|0_osc, 0_s⟩ +|0_osc, 1_s⟩/√(2)Û_ NL→|ζ/2, 0_s⟩ + |-ζ/2, 1_s⟩/√(2)
       R_y(π)→-|ζ/2, 1_s⟩ + |-ζ/2, 0_s⟩/√(2)Û_ NL|_π→-|ζ, 1_s⟩ + |-ζ, 0_s⟩/√(2)
       R_y(π/2)→ - N_+ (|ζ⟩ +|-ζ⟩) |0_s⟩ + N_- (|ζ⟩ -|-ζ⟩) |1_s⟩.

We can adjust the phases of the qubit rotations to adjust which harmonic oscillator superposition is heralded by|1_s⟩. For example, we can change the axis of theR_y(π)pulse toR_x(π). The sequence continued from Eq. (<ref>) then looks like

    R_x(π)→-i|ζ/2, 1_s⟩ -i|-ζ/2, 0_s⟩/√(2)Û_ NL|_π→-i|ζ, 1_s⟩ - i|-ζ, 0_s⟩/√(2)
       R_y(π/2)→ i N_- (|ζ⟩ -|-ζ⟩) |0_s⟩ - i N_+ (|ζ⟩ +|-ζ⟩) |1_s⟩.


For the trisqueezing superposition shown in Fig. <ref>, we omit theR(π/2)andR(π)pulses.
To herald the even (odd) superposition upon detecting|1_s⟩, we initialise the starting state in|1_s⟩(|0_s⟩)  before applying the interaction.


For the sequence in Fig. <ref>, we require all of the applied nonlinear interaction to be conditioned onσ̂_zsuch that we can treat its impact on the subsystem{|0_s⟩, |2_s⟩}as a closed system. 



 §.§ Adjusting the probability amplitude coefficient of the superpositions


We can arbitrarily set the probability amplitudes of each constituent instead of only balanced superpositions. 
To adjust the relative amplitudes, we can replace the firstR_y(π/2)rotation byR_y(θ)rotation such that the initial spin superposition is the form 

    |ψ⟩ = cos(θ/2)|0_osc, 0_s⟩ + sin(θ/2)|0_osc, 1_s⟩.

The rest of the sequence remains unchanged. This results in an amplitude ratiocos(θ/2)/sin(θ/2)between the constituents of the superposition. For example, forθ=π/2, the amplitudes are equal and the sequence is the same as in Eq. <ref>. 
As a proof of principle, we produce a squeezed state superposition withθ=π/4i.e.a |ζ_2⟩ + b|-ζ_2⟩with an amplitude ratioa/b = √(2)+1(see Fig. <ref>).


To adjust the relative phase of the complex probability amplitudes, we vary the axisσ̂_γ= cos(γ)σ̂_x + sin(γ)σ̂_yof theR_γ(π)rotation. After the full sequence, we obtain

    |ψ⟩   = -N_-^(γ)(|ζ⟩i e^iγ -|-ζ⟩ie^-iγ)|0_s⟩
       +N_+^(γ)(|ζ⟩i e^iγ +|-ζ⟩ie^-iγ)|1_s⟩.


 After rearranging, we obtain

    |ψ⟩ = -ie^iγ N_-^(γ)(|ζ⟩ -|-ζ⟩e^-i2γ)|0_s⟩
    
        + ie^iγ N_+^(γ)(|ζ⟩ +|-ζ⟩e^-i2γ)|1_s⟩.


In Fig. <ref>, we show the herald probability|N_+^(γ)|^2and the resulting Wigner functions when changing the relative phase of the probability amplitudeγ.


 §.§ Measuring the variance of squeezed superpositions

To estimate the variance in Fig. <ref>D based on the experimentally reconstructed Wigner functionW, we estimate the probability distributionP_xandP_pwhich can be found by integrating the Wigner function along the orthogonal quadraturespandxrespectively:

    P_x(x)    = ∫  W(x,p) dp,
    
        P_p(p)    = ∫  W(x,p) dx.

The variances are then

    Var(x̂)    = ⟨x̂^2|-⟩⟨x̂|^⟩2,
     Var(p̂)    = ⟨p̂^2|-⟩⟨p̂|^⟩2,

where the moments ofxandpare given by

    ⟨x̂^k|   ⟩= ∫ x^k P_x(x) dx,
    ⟨p̂^k|   ⟩= ∫ p^k P_p(p) dp.

For the experimental data, we rotate the matrix representing the discretised Wigner function such that the principle axes are aligned withxandp. For the simulated data, we extract a discretised Wigner function before performing the same analysis. We calculate the probability distributions and the moments by approximating the integration via the trapezoidal rule. 


§ EXTENDED EXPERIMENTAL DATA AND PARAMETERS


 §.§ Experimental parameters for the results in the main text

For Figs. <ref>, <ref>, and <ref> we used the experimental parameters tabulated in Tab. <ref>.


For Fig. <ref> we use the same parameter configuration as for Fig. <ref>B and adjust the number of shots such that we have more than300successful heralds per point.
For Fig. <ref> we use the same parameter configuration as for Fig. <ref>C and adjust the number of shots such that we have more than300successful heralds per point. 



 §.§ Mid-circuit detection error

The fidelity of the mid-circuit detection directly determines the degree of certainty to which we know the heralded state. We consider the imperfect detection in the case of heralding the odd superposition|ψ_k,-⟩. Just before detection, the wave function is 

    |ψ⟩ = N_k,-|ψ_k,-⟩|1_s⟩ + N_k,+|ψ_k,+⟩|0_s⟩.

With perfect detection, we could obtain the true description of|ψ_k,-⟩. With detection errors, if we herald on|1_s⟩(dark) the state reads

    ρ_ detect   ∝ P( dark| dark) |N_k, -|^2 |ψ_k,-⟩⟨ψ_k,-|
       + P( dark| bright) |N_k, +|^2 |ψ_k,+⟩⟨ψ_k,+|,

whereP(dark|dark)is the detection probability of dark given dark, which is the complementary event ofP(bright|dark)(false bright) andP(dark|bright)the probability of dark given bright (false dark).
Hence, the ratio of the mixturemof making the desired state, to the false dark state is given by

    m = P( dark| dark) |N_k, -|^2/P( dark| bright)|N_k, +|^2 ≈|N_k, -|^2/P( dark| bright)|N_k, +|^2 .

In Fig. <ref>, we measure and fit the histograms of our detection sequence. From the fits (see Fig. <ref>), we estimate the false dark probability to beP(dark|bright) = 4.9e-5. 
For the interactions considered the true herald ratio|N_k,-|^2/|N_k,+|^2 ≥0.1. Hence, given our readout errors, the resultingm > 10^3. If the true herald ratio is even lower, the experimental threshold can be adjusted at the expense of a small loss in overall detection events.
We can also measure a cumulative state-preparation and readout (SPAM) error by calculating the dark and bright probability based on the actual histogram data. We findP(dark) = 0.993. Hence, the probability of preparing and measuring dark is dominated by the state preparation process.

Another error contribution in our experiment was leakage of the 1033 deshelving laser. 

This results in a decay of the excited|1_s⟩state that is faster than the lifetime of the transition. This leakage is observed in the baseline measurement of Fig. <ref>A (green crosses). At 880, we observe a contrast loss of∼5%. This could be resolved by improving the extinction of the 1033 when it is switched off.



 §.§ Effect of mid-circuit detection on ion motion

We verify that performing a mid-circuit detection on a dark state (|1_s⟩) leaves the motional state unperturbed. We first prepare a regular squeezed state|ζ_2, 1_s⟩, followed by the mid-circuit readout before finally measuring the characteristic function along the principle axes of the squeezed state.
We investigate three different versions of the mid-circuit measurement: regular measurement (i.e. collect fluorescence for 200), adding 200 delay but not applying the fluorescence laser, and omitting the mid-circuit measurement altogether. The measurement of the characteristic function for the three different settings is shown in Fig. <ref>.  We conclude that the mid-circuit measurement does not perturb the motional state beyond increasing the sequence length by the mid-circuit measurement duration.


 §.§ Detection probabilities for equal superpositions shown in Fig. <ref>.

We repeat the detection probability measurements shown in Fig. <ref>A for the trisqueezed and quadsqueezed superpositions (see Fig. <ref>). A closed-form analytic solution does not exist for these superpositions <cit.>; we replace the theory line with a numerically simulated curve. 
We compare the detection probabilities of these superpositions to the probabilities obtained when no interaction is applied (i.e. system left idle).
Especially for the quadsqueezed state, the change in probability from leaving the system idle becomes comparable to the detection probability. This explains why the quadsqueezed superposition substantially deviates from the numerical prediction.


§ CHARACTERISING THE SUPERPOSITIONS


 §.§ Wigner negativity of different ideal superpositions

There are several measures to quantify the utility, or more precisely, the resourcefulness of the oscillator state superpositions. Here, we choose the Wigner logarithmic negativity (WLN) <cit.>
    𝐖(ρ) = log(∫ dx dp|W_ρ(x,p)|),

whereWis the Wigner-function of stateρ. 
WLN is a resource monotone, but it is still an active area of research under which circumstances it adequately quantifies the difficulty of simulating the state classically. 
We calculate the expected WLN𝐖of the ideal version of the equal superposition states studied in this work (coherent, squeezed and trisqueezed superpositions) and pure Fock states using a discrete version of Eq.(<ref>). We measure the negativity as a function of the mean phonon numbern̅ = ⟨N̂|_⟩oscfor the specific oscillator state.


 §.§ Comparing the realistic and ideal states

While we have good agreement between experiment and numerical simulation with realistic system parameters i.e. finite initial oscillator occupation ofn̅=0.1and non-negligible heating rate ofṅ̅̇= 300 quanta/s, there is a difference compared to an ideal system with perfect initial groundstate and no heating rate.
The macroscopic properties, such as the overall shape of the Wigner function, stay the same, but certain details or features become more or less accentuated. We qualitatively compare the realistic and ideal states in Fig. <ref>. The Wigner functions with finite temperature and heating effects are identical to the plots shown for Fig. <ref>HI.


We extract the WLN described in the previous section for the numerically simulated and experimentally reconstructed states.
The WLN for the experimental states is close to the value of the realistic simulation but reduced compared to the ideal states due to experimental imperfections. We present the extracted values in Table <ref>. However, we emphasise that these values should give a qualitative estimation but do not represent a quantitative characterisation. 
We would need to improve the Wigner function reconstruction protocol to deal with spurious ripples in the Wigner function for largerx,  pvalues, which arise due to shot noise and other experimental imperfections during the tomography. We have windowed the Wigner function around the origin to avoid biasing the results. The window size is set such that(∫_window W(x,p) dx dp)/(∫W(x,p) dx dp) =0.95. For reference, we also extract the minimum value of the Wigner functionmin(W), which is occasionally stated in the literature. However, we stress that this measure is bounded tomin(W)≥-1/πand only confirms that the state has Wigner negativity but does not measure resourcefulness. For example, the ideal Fock state|1⟩has a minimum value ofmin(W) = -1/π.


 §.§ Reconstruction of the Wigner function

We obtain the Wigner function via Fourier transform of the characteristic functionχ(β) ∈ℂ, following Refs. .
As an example, we provide the characteristic function of the even (Fig. <ref>B) and odd squeezed superposition states (Fig. <ref>C).
In Fig. <ref>, we show the experimentally measured characteristic function.
In the experimentally measured characteristic function, and consequently, in the reconstructed Wigner function, we find a tilt relative to thex,paxes of∼0.44, with daily variations on the order of0.06. This tilt is the result of a constant motional phase offset between the SDF used for reconstruction and the axis of the squeezing interaction. To make the comparison to the numerically simulated data easier, we correct the shift by digitally rotating the reconstructed Wigner function by this amount.



 §.§ Density matrices of equal superpositions in Fig. <ref>

To provide more intuition about the odd and even superposition states, we show the density matrices in the Fock basis (see Fig. <ref>). As expected only the Fock states|k·2n⟩for the even superpositions and|k·(2n + 1)⟩for the odd superpositions, wheren∈ℕ_0andkthe order of the nonlinear interaction, are populated.