Quantum dynamics of a few-photon parametric oscillator

Modulating the frequency of a harmonic oscillator at nearly twice its natural frequency leads to amplification and self-oscillation. Above the oscillation threshold, the field settles into a coherent oscillating state with a well-defined phase of either $0$ or $\pi$. We demonstrate a quantum parametric oscillator operating at microwave frequencies and drive it into oscillating states containing only a few photons. The small number of photons present in the system and the coherent nature of the nonlinearity prevents the environment from learning the randomly chosen phase of the oscillator. This allows the system to oscillate briefly in a quantum superposition of both phases at once - effectively generating a nonclassical Schr\"{o}dinger's cat state. We characterize the dynamics and states of the system by analyzing the output field emitted by the oscillator and implementing quantum state tomography suited for nonlinear resonators. By demonstrating a quantum parametric oscillator and the requisite techniques for characterizing its quantum state, we set the groundwork for new schemes of quantum and classical information processing and extend the reach of these ubiquitous devices deep into the quantum regime.

Parametric amplifiers and oscillators are quintessential devices used to amplify small electromagnetic signals [1,3,5], convert radiation from one frequency to another [4,5], create squeezed light and entangled photons [6,7], and realize new information processing architectures [8][9][10][11][12]. They operate by modulating a parameter, the natural frequency of the resonant circuit ω c , at approximately twice its frequency ω p ≈ 2ω c . The modulation amplifies one of the field quadratures at the half-harmonic frequency ω p /2. A sufficiently large amplification overtakes the detuning and decay present in the system and leads to an exponential increase in the half-harmonic cavity field amplitude. Nonlinearities clamp this exponential growth and cause the system to enter an oscillating steady-state. These nonlinearities can be dissipative or dispersive. An example of the latter is the Kerr nonlinearity that induces a change in the cavity frequency proportional to the intracavity field intensity or photon number. Furthermore, the self-oscillation amplitude scales inversely with the magnitude of the nonlinearity. These nonlinearities have been exceedingly small in parametric oscillators to date, leading to large oscillation amplitudes that result in rapid decay of quantum coherence and the appearance of classical dynamics [13]. The quantum regime of nonlinear parametric oscillators has been extensively studied in theory [6,14,15] but received only limited attention in experiments [17].
We experimentally realize a quantum Kerr parametric oscillator (KPO) by implementing an on-chip superconducting nonlinear resonator and investigate its quantum dynamics under parametric driving. In contrast to previously demonstrated optical and microwave parametric oscillators, our device operates in the quantum regime Quantum or Classical FIG. 1. Schematic of a parametric oscillator above the instability threshold. An LC circuit with a harmonically modulated inductance is an example of a parametric oscillator. Above threshold, the system is oscillating at half the driving frequency ωp/2, which can be described as a parametric down-conversion process. The phase of the oscillator with respect to the external driving can be 0 or π. Classically, this symmetry is broken and one of the two cases is realized at random. Quantum mechanically, the oscillator can be in a superposition of both states.
with a self-oscillating state containing only a few photons. The resonator is implemented as an LC circuit with an array of Josephson junctions in place of the inductor [18]. for different driving amplitudes at ∆/2π = 11.2 MHz. The spectrum goes from a double-peaked shape at small β to a single narrow peak in the self-oscillation regime, passing through a multi-peaked spectrum at intermediate values of β. The white dot indicates the approximate on-set of self-oscillation. The spacing between the two peaks in the weak drive limit is 2∆, as explained by the two-photon emission process shown in (E).
down from 8 GHz to below 4 GHz by an on-chip flux line and most of the measurements are done with the resonator in the 6 − 8 GHz frequency range. This tunability also enables parametric driving of the formβ(t)(â +â † ) 2 whereβ(t) is proportional to the voltage V (t) applied to the flux line. With parametric drivingβ(t) = 2β(t) cos ω p t at frequency ω p = 2(ω c − ∆), which is slightly detuned from the parametric resonance 2ω c , the dynamics of the resonator in a rotating frame at half the driving frequency is well described by the Hamiltonian where β(t) is the slowly varying amplitude of the parametric driving. We first characterize the energy level structure of the Kerr nonlinear resonator without parametric driving (β = 0). To do this, we prepare the resonator in an excited stateρ 0 and let it relax back to the vacuum state while we collect the emitted microwave signal. The power spectral density (PSD) of this signal, which we call "transient PSD" to distinguish it from the PSD measured at steady state, contains multiple peaks evenly spaced by χ (Fig. 2B) due to the nonlinear energy level structure. The nth peak captures photons at frequency ω c −(n−1)χ, emitted during relaxation from |n to |n − 1 . We measure the transient PSDs for coherent statesρ 0 = |α α|, which we prepare by driving the system resonantly at ω c with 1 ns pulses of different amplitudes. Since the pulses are much shorter than 1/χ, they simply displace the state of the resonator from vacuum into a coherent state |α , with α proportional to the pulse amplitude. In agreement with theory, the number of peaks we observe in the spectrum grows with α as higher Fock states are populated (Fig. 2B).
The transient PSD measurement provides a way to infer the Fock occupations p n = n|ρ 0 |n of the initial statê ρ 0 . The probability that a transition from |n to |n − 1 occurs during free relaxation to |0 is exactly equal to the total population of all Fock states |n and higher. When the peaks are resolved (χ κ), the total power in the nth peak is proportional to ∞ k=n p k (see SI). The measured signal is related to the field emitted from the resonator by a frequency-dependent gain factor. The calibrated gain allows us to calculate the Fock occupations {p n } in the initial state from the transient PSD.
The quantum parametric oscillator does not have a sharp self-oscillation threshold. To study the transition to self-oscillation, we modulate the flux line with a continuous wave at a fixed frequency ω p = 2(ω c − ∆) and measure the mean photon number (Fig. 2C) as well as the steady state PSD (Fig. 2D) for increasing parametric drive amplitudes β. The mean photon number n ( Fig. 2C) is obtained by turning the drive off and measuring the transient PSDs as the state relaxes. At large β, both the classical and the quantum model predict photon numbers very close to the measured results. For smaller β, the measurements show a gradual transition into self-oscillation, smoothed by quantum fluctuations. This behavior deviates from the classical prediction but is reproduced well by the quantum model.
The steady state PSD (Fig. 2D) spectrum reduces to two peaks spaced by 2∆ in the weak drive limit and can be understood as the result of a two-photon emission process, with one photon emitted at the cavity frequency ω c = ω p /2 + ∆ and the other, by energy conservation, at ω p /2 − ∆ (Fig. 2E). For large drive amplitude β the resonator enters the self-oscillation regime and is phaselocked to the external parametric drive leading to the narrow peak seen in Fig. 2D. The narrowing of this peak results from the reduction of the switching rate between the two possible self-oscillating states (0 and π phase) as we increase β. Characterizing the full quantum state of the parametric oscillator is challenging and calls for a distinct approach to quantum tomography. There are several methods based on measuring Fock state populations or parities of displaced states. These approaches were developed for qubit-resonator systems [19][20][21][22] or qubits with tunable coupling to the environment [23] and require more complicated devices with auxiliary cavities and control parameters. Other methods have been developed in quantum optics based on studying the statistics of the output fieldâ out of linear resonators [24][25][26], whereâ out is linearly related to the resonator fieldâ(t = 0) at some initial time. These methods are not suitable for our system since χ κ leads to a highly nonlinear relationship between the output field and the intracavity mode.
We develop a state tomography method suited for the quantum parametric oscillator that does not need any auxiliary systems. Our method is based on measuring the transient PSD using a purely linear detection of the output field. The transient PSD measurement provides us with the diagonal elements ofρ. By displacing the state in phase space and detecting the transient PSD, we find the diagonal elements of a displaced density matrix DρD † that contain information about the off-diagonal elements ofρ. Repeating this for several different displacements, we obtain enough information to estimate the full density matrix by a maximum likelihood method. This aspect of the tomography method is conceptually similar to an existing technique using generalized Q functions [19,27]. We arrive at the estimate ofρ by minimizing a loss function L(ρ est ), which quantifies the difference between the measured PSD and the PSD simulated for the stateρ est (see SI). This convex optimization problem can be solved efficiently by semi-definite programming. The broadening of the higher Fock state peaks due to their shorter lifetime limits the maximum size of an unknown state that we can reconstruct. Additionally, we have observed systematic errors which become more severe for larger measured states. We suspect this to be due to excitation of unwanted off-resonant transitions by strong drive pulses and other nonlinear effects.
We use this tomography procedure to study the free dynamics of the system. After displacing the resonator from vacuum to a coherent state, we let the system evolve freely for a time τ before performing the tomography (Fig. 3A). Here as well as in the rest of the paper, the displacements used in the tomography procedure are arranged in concentric rings as illustrated in Fig. 3B. Fig. 3C shows an example of a raw tomography dataset, consisting of the transient PSDs for all displacements. We reconstruct the state at different time slices τ and observe the evolution of the Wigner distributions due to the Kerr nonlinearity (Fig. 3D). In agreement with simulations, this evolution involves a "shearing" distortion, which has a classical counterpart explained by the amplitude-dependence of the resonator frequency, and also exhibits negative Wigner function values, a signature of quantum mechanical behavior [19].
The system can behave nonclassically under a parametric drive, but only before photons leaking out of the oscillator cause the loss of quantum coherence. In our system, quantum dynamics persist in this transient regime for a sufficiently long time to allow detailed observation. This nonclassical evolution is important to understand in the context of emerging applications for parametric oscillators in quantum information processing [10,11]. To investigate the transient dynamics, we turn on the parametric drive suddenly and measure the time evolution of the mean photon number for three different drive amplitudes, β/2π = 3.5, 5.8 and 11.5 MHz ( Fig. 4A). The observed time dependence of the mean photon number is in good agreement with a theory fit to all three data sets simultaneously, using only the detuning ∆, the loss rate κ and a drive conversion factor β/V as free fit parameters. Using the previously described tomography procedure, we also reconstruct the quantum state for β/2π = 5.8 MHz at t = 20 ns (transient state) and t = 2000 ns (steady state) (Fig. 4A). Comparison between the theoretically predicted and experimentally obtained Wigner functions shows relatively good agreement with fidelities of 0.93 and 0.94 for the transient and steady states respectively. We attribute the discrepancies to systematic errors in the tomography process which we believe could be mitigated by improving its calibration. The only fit parameters in the theoretical prediction are the overall rotation angle and a short delay (t d = 2.5 ns) between the end of the parametric drive pulse and the start of the tomography (see SI).
These results demonstrate that the state of the oscillator can be engineered by designing the parametric drive β(t). For example, by adiabatically changing β(t), we can prepare even-parity energy eigenstates of the Hamiltonian for different values of β as long as losses are negligible [6]. Intriguingly, for ∆ < 0, as β approaches and exceeds χ, the energy eigenstate adiabatically connected to the vacuum state begins to closely approximate the even-parity Schrödinger's cat state. We set the pump detuning in the experiment to ∆/2π = −6.7 MHz, and begin with the resonator in the vacuum state with β = 0.
We increase β slowly so the system follows the eigenstate of the Hamiltonian. Numerical simulations suggest that given our system's parameters, with t max = 22 ns (Fig. 4B) is a reasonable parametric driving profile for preparing a cat state. The length of this signal is much shorter than the cavity decay time 1/κ ≈ 150 ns but long enough to ensure approximately adiabatic evolution. We perform the experiment ramping to different values of β max and verify the result with the state tomography procedure. The results are compared to the Wigner functions found theoretically (Fig. 4D). In the simulations, the rotation angle and the drive conversion factor β/V are the only fit parameters and common to all four data sets. As described above, we again assume a short (t d = 2.5 ns) period of free evolution between the end of state preparation and the start of our tomography, which causes a small distortion (due to χ) of the reconstructed state with respect to the eigenstate of the driven system. In our experiment, we ramp the drive up to β max ∼ χ, which is necessary to see the emergence of the cat state. To a good approximation, the generated states after the short free evolution, described byÛ 0 (t d ), are with α = 0.64, 0.88, 1.08, and 1.2 for the data shown in Fig. 4D. The largest of these corresponds to a 4|α| 2 = 5.8 photon Schrödinger's cat state [22,28]. In the β max χ regime, the relevant eigenstate exponentially approaches a cat state of size 8β/χ due to the double-well shape of the system's effective potential (Fig. 4C) (see SI for more details). For very large cat states, imperfections in the tomography process that grow with the number of photons in the analyzed state prevent its faithful reconstruction. In comparison to other schemes for cat state generation [22,[28][29][30][31], our method is significantly more hardware efficient as it requires only a resonator with one input and one output line for both state generation and read-out.
The parametric oscillator is one of the paradigmatic systems in quantum optics and has found an enormous range of applications over the years. We have experimentally demonstrated the few-photon quantum dynamics of a parametric oscillator by introducing a large Kerr nonlinearity in the microwave frequency regime. We show that the system can adiabatically generate cat states of five to six photons and have developed a tomography method suited for the characterization of its state and dynamics. Our work demonstrates that nontrivial quantum states can be engineered and characterized with nearly minimal hardware complexity. The quantum coherence and hardware efficiency of the system bode well for the prospects of scaling these devices to larger networks in emerging applications of parametric oscillators for quantum information processing and optimization [8][9][10][11][12]. The device, shown in Fig. S1A, was fabricated using a 5 mask lithography process on a 500-µm high-resistivity (> 10 kΩ · cm) Si substrate. First, the aluminum ground planes and feed lines are defined in photolithography using a liftoff process. Next, palladium marks are added in preparation for aligning subsequent electron-beam (ebeam) lithography masks. The SQUID array shown in Fig. S1B is fabricated using a Dolan-bridge double-angle technique for growing Al/AlO x /Al junctions via in situ oxidation [S1, S2]. After junction growth, the resonator capacitor is defined using e-beam lithography; narrow wires and an unconventional capacitor design were chosen to accommodate an array of nanomechanical resonators introduced in later devices [S3], and are not essential to the design. Finally, a superconducting connection between the SQUID array and capacitor leads is formed using a bandage process [S4]. An equivalent circuit diagram of the device is shown in Fig. S1C. Table I gives the device parameters for the Kerr parametric oscillator. Here, the maximum resonator frequency ω c,max is determined from a fit to the flux-bias tuning curve ω c (Φ e ) = ω c,max | cos(πΦ e /Φ 0 )|, where Φ e is the externally applied magnetic flux and Φ 0 is the flux quantum. The maximum Josephson energy of each of the N = 10 SQUIDs, denoted by E J,max is determined from normal-state resistance measurements. Together, E J,max and ω c,max are used to extract the resonator charging energy E C , which closely matches predictions from finiteelement capacitance simulations. The resonator intrinsic and extrinsic decay rates are denoted as κ i and κ e , respectively. χ is the resonator frequency shift per photon, determined from the peak-to-peak splitting in transient PSD measurements of a coherent state. The maximum parametric drive amplitude β/2π used in the experiment is about 20 MHz. Larger β are possible but lead to large states that cannot be faithfully reconstructed by the employed tomography procedure due to the nonlinearity of the system.

Up-conversion board
We use two separate channels of a Tektronix series 5200 arbitrary waveform generator (AWG) to synthesize the temporal profile of the pulses at an intermediate frequency around 4 GHz and then further up-convert them to the desired frequencies close to the first (for the displacement) and second (for the parametric drive) harmonic of the resonator frequency using single sideband mixers (Fig. S2). Both pulses are then combined together and sent to the sample through the flux line. Thanks to a weak but nonzero direct coupling of the flux line to the resonator, it can effectively double as a weakly coupled charge line which we can use instead of sending the displacement pulses through the resonator input/output line. This way, we avoid saturation of the measurement setup by the reflection of the strong resonant pulse. The input/output line is not used for driving the system, except in initial characterization measurements of the resonator frequency using a vector network analyzer.

Phase locking
For a state prepared by parametric driving, the orientation of the reconstructed quasiprobability distribution in phase space is determined by the difference between the phase φ d of the tomography pulses and the phase φ p /2 of the ω p /2 subharmonic of the parametric driving. These in turn depend on the absolute phases of the up-conversion local oscillators. Slow changes of this relative phase due to phase drifts of the signal generators would lead to gradually accumulating errors over long measurements. To mitigate this, we monitor the phase over the course of the measurement using a second downconversion board. Here we somewhat unconventionally feed both the parametric driving pulse at ω p and the displacement pulse at ω d into the RF port of the mixer while the LO port is 50 Ohm terminated (Fig. S2), relying on intermodulation to produce a signal at ω p −2ω d . By measuring this signal, we can determine the phase difference φ = φ p − 2φ d which needs to be constant to ensure correct performance of the tomography measurements even with very long acquisition times. We observe a slow drift of φ of about 1 radian per hour which we then eliminate by measuring φ in approximately 1 minute intervals and appropriately adjusting the phase of the up-conversion LO.

Down-conversion and digitization of the signal
The signal emitted by the resonator at about 7 GHz is first amplified by a traveling wave parametric amplifier (TWPA) [S5] and a low-noise high electron mobility transistor (HEMT) amplifier inside the cryostat. At room temperature, it is further amplified and then converted to an intermediate frequency of 125 MHz by a singlesideband mixer (Fig. S2). This is then recorded by a digitizer card (AlazarTech ATS9350) with a 12-bit resolution and a sampling rate of 500 MS/s. The acquired data is first saved in the on-board memory buffer and then transferred to GPU for real-time data processing.

A. Hamiltonian of a parametric oscillator
The Hamiltonian of the SQUID array resonator including the parametric driving iŝ wheren is the number of Cooper pairs andφ the overall phase across the junction array. E C is the resonator's charging energy, N is the number of SQUIDs in the array and E J is the Josephson energy for a single SQUID in the array which depends on the external flux Φ(t). The flux is harmonically modulated around its mean value with a small amplitude such that E J (Φ(t)) can be approximated as E J + δE J cos ω p t. After Taylor-expanding cos(φ/N ) to fourth order, we get The quadratic time-independent part of the Hamiltonian can be diagonalized by defininĝ We also drop c-valued terms in the expression above and get where ω c (0) = 8E C E J /N . We then transform the Hamiltonian into a rotating frame at the frequency ω p /2, perform a rotating wave approximation assuming |ω c (0) − ω p /2| ω c (0) and normal-order the resulting expression, which giveŝ where we have defined the Kerr nonlinearity χ = E C /N 2 , the parametric drive strength β = ω c (0) δE J /8E J , the dressed resonator frequency ω c = ω c (0) − χ and the detuning ∆ = ω c − ω p /2.

Effective potential
An intuitive way to understand some aspects of the dynamics of this system is to define the following effective potential [S6] in phase space where α = |α|e iθ . In the regime 2β > |∆|, the effective potential has two symmetric local maxima (corresponding to classical stationary points of the system) at |α| > 0 and θ ∈ {0, π}. Therefore in this case the effective potential has the shape of an inverted double-well. (Fig. S3) When its maxima ±α are sufficiently well separated, i.e., β is large, the potential can be approximated by a quadratic function near ±α. Therefore the eigenstates of the system are close to Fock states displaced by ±α, as described in more detail below.

Eigenstates of the Hamiltonian
Understanding the structure of the eigenvalues and eigenstates of the Hamiltonian at different β helps us see how adiabatic ramping up of the drive amplitude can lead to Schrödinger's cat state generation [S6]. At β = 0, the eigenstates of the Hamiltonian are simply all Fock states {|n }. As we increase β, two adjacent energy levels get closer and eventually merge together at very large β (Fig. S4). In the limit where β/χ 1, as suggested by the effective potential argument above, the eigenstates where there is only one global maximum while (B) represents the regime β χ (far above threshold) where the potential has two symmetrically placed local maxima.
of the Hamiltonian form many two-dimensional nearly degenerate subspaces spanned by {D(±α) |n } for each n, where ±α = ± 2β/χ are the locations of the local maxima. To show this, we first note that the coupling between states of the formD(+α) |n andD(−α) |m un- S4. Eigenenergy of the Kerr parametric oscillator with different driving strength. The plot is made in the classical limit χ ∆ where the oscillation threshold happens at β = |∆|/2. The nearby energy levels merges together and becomes nearly degenerate at large β where the degenerate subspaces are spanned by displaced Fock states. When we adiabatically increase β from 0 to above threshold, the vacuum state will evolve into an even cat state since parametric driving preserves parity.
Intuitively, this follows from the large separation of the two potential wells. The Hilbert space therefore effectively decomposes into two nearly decoupled subspaces consisting of states localized around +α and −α. Next, we observe that the Hamiltonian within each of these subspaces is close to a harmonic oscillator. That is, D † (α)ĤD(α) ≈D † (−α)ĤD(−α) ≈ −2χα 2â †â in the limit of α 1. Here we have used the relation α ≈ 2β/χ and only kept the highest order terms in α. This confirms that displaced Fock states indeed approximate the pairwise degenerate eigenstates ofĤ in the large α limit.
It follows from the adiabatic theorem that if we prepare the system in the vacuum state and adiabatically increase β, the system will follow the eigenstate and end up within the corresponding degenerate subspace, as determined by the order of eigenenergies of the Fock states at β = 0, which in turn depends on ∆ and χ. For ∆ < 0, |0 has the highest eigenenergy at β = 0 and evolves into {D(±α) |0 } at large β. Since parametric driving preserves the parity of the state, the final state is an even cat state |α +|−α . The case ∆ > 0 is more complicated since |0 may not have the highest eigenenergy at β = 0 and can thus evolve into some [D(α) +D(−α)] |n where n = 0. Therefore for cat state generation, having a negative detuning is helpful since that gives a lager energy gap between |0 and all other higher Fock states and therefore allows a faster adiabatic tuning of β. Driving with an appropriately chosen positive detuning or non-adiabatic drive variations could on the other hand be useful for generating displaced Fock states or their superpositions [S6].

B. Displacement pulse calibration
The displacement pulses used in the experiments are 1 ns short pulses created by an arbitrary waveform generator. To calibrate them, we apply pulses with different amplitudes {V i } to the vacuum state of the resonator and measure the transient PSDs for the generated states. In order to simplify the analysis of the PSD measurements, we preprocess the raw data {S(ω; V i )} by integrating over frequency bins centered around each of the individual transition peaks, thus effectively reducing the dimensionality of the analyzed data to n × m, where n is the number of bins and m the number of different pulse amplitudes.S The number of bins which can be usefully analyzed is limited by the increasing overlaps between the peaks corresponding to higher transitions with larger linewidths. When multiple transitions contribute to the same bin, the assumptions we use in our model to arrive at Eq. (S16) fail and a more complex model needs to be used. In most of our measurements, we have used n = 4 to 5 bins.
The calibration is done under the following assumptions: • The state generated by a single pulse with voltage V i is a coherent state |α i .
• α i depends linearly on the voltage V i , i.e., α i = kV i , where k is a single fit parameter common to all pulses.
• The bin powers {S j (V i )} calculated from the measured PSDs are related to the theoretical predictions {S j (|α i )} by a gain factor c j which may in principle be different for each bin j.
The calibration parameters k and c = (c 1 , . . . , c n ) are obtained by minimizing the loss function (Fig. S5A) For a given k, finding optimal c reduces to a simple linear fitting problem and the value of k is then calculated by minimizing L(k) = min c L(k, c).
To evaluate the loss function above, we need to calculate the theoretically expected total power S j (ρ 0 ) in each transient PSD peak for a given stateρ 0 . The result represented by Eq. (S16) was outlined in the main text and its full derivation follows below.

Analytical formula for the transient PSD
The transient PSD for an initial stateρ 0 is given formally by From the quantum regression theorem, we can rewrite the two-time correlation function as â † (t + τ )â(t) = Tr â † eL τâ eL tρ 0 (S12) whereL is the Liouvillian for the nonlinear resonator with the HamiltonianĤ = − χ 2â †â †ââ (S13) and energy decay rate κ.
By solving the master equation and Fourier transforming the two-time correlation function, the analytical result for the transient PSD is where u = iχ/κ and p n = n|ρ 0 |n . Notably, S(ω;ρ 0 ) depends only on the diagonal part ofρ 0 and this dependence is linear.
Our system satisfies |u| = χ/κ 1. In this regime, the transient PSD takes the form of a sum of relatively well separated peaks. The last equation can be approximated as S(ω;ρ 0 ) ≈ ∞ n=0 p n n−1 j=0 (2j + 1)κ (ω + jχ) 2 + [(2j + 1)κ/2] 2 (S15) where the transient PSD is decomposed into a sum of Lorentzians with different linewidths and different cen-ter frequencies corresponding to the different peaks we measured. Therefore the total power in the jth peak is expected to be (S16)

C. State tomography
To reconstruct an unknown quantum stateρ, we displace it in phase space (Fig. 3A) using short calibrated pulses with complex voltages {V i } (including phase), then measure the transient PSDs and calculate the corresponding integrated powersS j (V i ) for each of the bins j = 1, . . . , n after each displacement i = 1, . . . , m. The density matrixρ can be estimated by minimizing the difference between the predicted transient PSDs and the measured ones, which is expressed by the loss function ( Fig. S5B) (S17) under the linear constraint Tr(ρ) = 1 and under the condition thatρ is positive semidefinite. Notice that both the displacement and the map from a density matrix to a transient PSD are linear transformations and therefore this minimization problem is convex and can be efficiently solved by the Matlab package CVX [S7, S8].

Parity constraint
The form of the master equation of the parametrically driven system implies that all states which can be generated from a vacuum state are mixtures of states with even and odd parity. To see this, we only need to observe that the Hamiltonian conserves parity and the collapse operator √ κâ flips the parity of a state. When reconstructing states prepared by parametric driving, we use this condition as an additional constraint on the unknown density matrixρ, requiring that whereP = e iπâ †â is the parity operator. We justify this assumption by verifying the corresponding symmetry of the transient PSDs measured in the tomography process under a rotation of the applied displacement by π. For a stateρ satisfying Eq. (S18), we expect that S(ω; +α) =S(ω; −α) and we check this by plotting the differenceS(ω; +α) −S(ω; −α) and observing that it is negligible when compared with the PSDsS(ω; α) themselves (Fig. S6). but for the tomography dataset in (Fig. 4B right).

D. Data fitting
Some parameters of the system such as its nonlinearity χ and linewidth κ do not change significantly among different experiments and we therefore assume constant values for them which are obtained from initial characterization measurements. Other parameters like the detuning ∆ and the parametric driving amplitude β vary between experiments and therefore their values are determined separately in each instance, either directly from the settings of the experiment or by fitting.
The nonlinearity χ/2π = 17.3 MHz was calculated from transient PSD measurements as the mean spacing between adjacent peaks. To fit the time evolution of the photon number (Fig. 4A), we introduce the detuning ∆, the linewidth κ and the drive conversion factor β/V that relates the voltage amplitude V applied to the flux line to the parametric driving amplitude β as fit parameters and get ∆/2π = 24.6 MHz and κ/2π = 1.1 MHz by minimizing the L 2 distance between the simulation results and the measured data. The value of ∆/2π found by fitting is very close to the value 24 MHz set in the experiment and the small difference is likely due to slow magnetic flux noise which causes variations in the resonator frequency. The linewidth is also consistent with direct VNA measurements, which give values around 1 MHz, slightly depending on the resonator frequency. Since none of the results in this work are very sensitive to the exact value of the linewidth, we fix κ/2π to be 1.1 MHz in all subsequent theory fits.
The process of fitting the steady-state mean photon number n (Fig. 2C) is similar to the case of the timedependent n measurement described above (Fig. 4A) except that κ is fixed and the detuning resulting in the best fit is ∆/2π = 25.3 MHz. For fitting the PSDs at steady state (Fig. 2D), we fix ∆/2π to the value 11.2 MHz set in experiment. The only fit parameter is the drive conversion factor β/V . The tomography measurement of the freely evolving state (Fig. 3C) has two fit parameters: the size and phase of the coherent state at τ = 0. In principle, the observed phase θ of the state should be easily predictable since it only depends on the relative phase between the preparation pulse and the tomography pulse, both of which are generated by the AWG and undergo the same path through the up-conversion chain. A calculation based on the used experimental settings gives a predicted phase of θ ≈ 1.19π. The size of the state α can also be estimated from the pulse calibration parameter k, which gives α ≈ 1.5. Treating both θ and α as unknown fit parameters, we get θ ≈ 1.24π and α ≈ 1.0. In this fitting, to achieve simultaneous match to the different states at each of the different evolution times τ , we choose the objective function to be the geometric mean of the fidelities between each measured state and the corresponding theoretical prediction.
For states prepared by parametric driving, their phase θ depends on the absolute phase of the signal generator. Through the feedback loop described in Section I C 2, we stabilize the phase over the measurement time at a fixed value. This value in principle depends on the frequency of the signal in a complex way due to the variation of the system's S-parameters with frequency. Since the different tomography measurements are mostly performed at different resonator frequencies, we treat the phase of the state θ at each of these frequencies as a fit parameter. Another fit parameter we introduce for these measurements is the possible small time interval t d of free evolution between the preparation pulse and the tomography pulses. This is to account for a potential delay between the two pulses which are generated by different channels of the AWG and processed by separate up-conversion boards.
For the tomography measurements of the steady state under parametric driving (Fig. 4B right), both the detuning ∆ and the parametric driving amplitude β have been fixed through fitting to the time evolution of the photon number (Fig. 4A). Therefore the only remaining fit parameters are the phase of the state and the delay time t d . By maximizing the fidelity between the reconstructed density matrix and the simulation result, we get t d = 2.5 ns, which is then kept fixed for all other tomography measurements with parametric driving. Consequently, for the transient state at 20 ns (Fig. 4B left), the only fit parameter is its phase θ.
For the tomography measurements of the adiabatically prepared cat states, the phase θ is again unknown but should be the same for all four states. Therefore the phase θ and the drive conversion factor β/V are the only fit parameters for all four data sets. We again assume that the final state is distorted by a short period of free evolution whose length we fix at t d = 2.5 ns based on previous measurements.