Engineering Arbitrary Hamiltonians in Phase Space

We introduce a general method to engineer arbitrary Hamiltonians in the Floquet phase space of a periodically driven oscillator, based on the non-commutative Fourier transformation (NcFT) technique. We establish the relationship between an arbitrary target Floquet Hamiltonian in phase space and the periodic driving potential in real space. We obtain analytical expressions for the driving potentials in real space that can generate novel Hamiltonians in phase space, e.g., rotational lattices and sharp-boundary well. Our protocol can be realised in a range of experimental platforms for nonclassical states generation and bosonic quantum computation.

Another area of interest for Hamiltonian engineering is topology.Due to the non-commutative nature of phase space, a quantum particle moving on a closed phasespace loop acquires a geometric phase analogous to the Aharonov-Bohm phase for a particles in a magnetic field.As a consequence, a gapped lattice Hamiltonian in phase space can support non-trivial Chern numbers [16,[32][33][34][35][36][37][38][39][40].This is an appealing feature because in a system with a physical boundary, it would lead to topologically robust edge transport.While it has been shown how to generate arbitrary lattice potentials in phase-space [41], so far it was unclear how to combine such a potential with a sharp phase-space confinement.
It is well known that the stroboscopic dynamics of any periodically driven system can be described in terms of a time-independent Floquet Hamiltonian ĤF defined via exp 1 iλ ĤF T ≡ Û (T, 0) = T exp 1 iλ Here, Û (T, 0) is the time-evolution operator with T the time-period of the system's time-dependent Hamiltonian Ĥ(t).In adddition, λ is an effective dimensionless Planck constant, and T is the time-ordering operator.Except for very few models, it is impossible to obtain a closed form of the Floquet Hamiltonian ĤF from the time-dependent Hamiltonian Ĥ(t).Instead, one often evaluates the Floquet Hamiltonian relying on a high-frequency expansion [42][43][44], e.g. the Magnus expansion theory [45,46], the van Vleck degenerate perturbation theory [47] and the Brillouin-Wigner perturbation theory [48].In this work, we focus on the inverse problem, that is, to find the timedependent Hamiltonian Ĥ(t) that synthetizes a target Floquet Hamiltonian Ĥ(T ) F .This is the realm of Floquet engineering which is a very developed and active field [39,[49][50][51][52].Most of the work so far has focused on implementing specific Floquet Hamiltonians of interest.However, a systematic constructive method to solve the inverse Floquet problem for a single quantum particle is still missing.In this work, we provide such a method.

II. MODEL AND GOAL
As a starting point, we consider a periodically driven oscillator with lab-frame Hamiltonian Here, ω 0 is the oscillator natural frequency, β is the amplitude of the nonlinear driving potential V (x, t) which has time-period T d and might contain also static terms.
In order to introduce an effective dimensionless Planck constant λ [53][54][55], the position x, the momentum p and Ĥ(t) have been rescaled such that [x, p] = iλ and at the same time the Schrödinger equation reads iλ ψ = Ĥ(t)ψ.Parameter λ measures the quantumness of our system and λ → 0 corresponds to the classical limit.
The Floquet Hamiltonian to be designed has timeperiod T = 2π/ω 0 and is defined via Ô(t) = exp[iâ † âω 0 t], where â = (x + ip)/ √ 2λ is the annihilation operator.In other words, Ĥ(t) in Eq. ( 1) is the rotating-frame Hamiltonian given by Ĥ(t) = Ô(t) Ĥ(t) Ô † (t) + iλ Ȯ(t) Ô † (t), arXiv:2302.04257v3[quant-ph] 19 Nov 2023 which in our case reads Ĥ(t) = βV [x cos(ω 0 t) + p sin(ω 0 t), t]. (3) We enforce the time-periodicity, Ĥ(t) = Ĥ(t + T ), by setting T = qT d with q ∈ N ≥ 1, corresponding to a qphoton resonance.Any detuning from the multiphoton resonance is formally incorporated in the driving potential V (x, t).For weak nonlinearity, β ≪ ω 0 , the evolution in the rotating frame is slow.Thus, we are in the realm of application of the Floquet high-frequency expansions, here, with the small parameter β/ω 0 .This allows us to approximate the Floquet Hamiltonian with the leading order of the Floquet-Magnus expansion corresponding to the rotating wave approximation (RWA), Our goal is to engineer an arbitrary target Floquet Hamiltonian Ĥ(T ) F in phase space by properly designing the driving potential V (x, t) in real space.Up to leading order (RWA) we, thus, require that the right-hand side of Eq. ( 4) coincides with the target Hamiltonian Ĥ(T ) F .The ensuing solution becomes exact in the highfrequency limit ω 0 /β → ∞.

III. NcFT TECHNIQUE
As a preliminary step towards deriving a suitable driving potential V (x, t), we introduce a useful decomposition of the target Hamiltonian Ĥ(T ) F in the form of a noncommutative Fourier transformation (NcFT).This can be viewed as a variant of quantum distribution theory [56].We wish to decompose the target Hamiltonian Ĥ(T ) F as a sum of plane-wave operators Ĥ(T ) It can be shown that the Fourier coefficients f T (k x , k p ) are given by the inverse transformation [see App.A] where the phase-space function H (T ) Q (x, p) is the equivalent of the Husimi Q-function, here, for a Hamiltonian instead of the density operator.We remind that the Q-function of an operator evaluated at a phase space point (x, p) is simply its expectation value in the corresponding coherent state, H (T ) Q (x, p) = ⟨α| Ĥ(T ) F |α⟩ with â|α⟩ = α|α⟩ and α = (x+ip)/ √ 2λ.The latter mean value can be calculated by normal ordering the target Hamiltonian Ĥ(T ) F (â † , â).We point out three important features of the Hamiltonian Q-function: (i) For fixed λ, the  9); (Right) the engineered real-space potential V (x, t) for parameters q = 6 and λ = 0.01.The white contours indicate the minima of the instantaneous real-space potential.
mapping between Floquet Hamiltonians and Q-functions is one-to-one; (ii) The Hamiltonian Q-function has the same phase-space symmetries as the corresponding Floquet Hamiltonian, and (iii) a well-defined classical limit

IV. DESIGNING DRIVING POTENTIAL
The driving potential V (x, t) that generates the target Floquet Hamiltonian H (T ) F (x, p) can be readily obtained from its Fourier coefficient f T (k x , k p ).We can formally write the solution as a superposition of sinusoidal potentials (7) with time-varying amplitudes A(k, τ ) and phases ϕ(k, τ ) determined from the Fourier coefficients in polar coordi- This solution can be readily verified by plugging it into Eqs.( 3) and (4), and changing the integration variables back to cartesian coordinates to arrive at Eq. ( 5) [see App.C].In the remainder of this paper, we demonstrate the flexibility of our method by calculating the potential V (x, t) for a range of interesting Floquet Hamiltonians.In passing, we will also highlight more general features of our solution and comment on certain associated subtleties.

V. EXAMPLES A. Rotational lattice
We now apply our method to engineer a particularly interesting Floquet Hamiltonian with q-fold symmetry in phase space Ĥ(T ) The discrete rotational symmetry can be described by R( 2π q ) Ĥ(T ) F R † ( 2π q ) = Ĥ(T ) F , where R(θ) = exp(iâ † âθ) is a phase-space rotation by an angle θ [20,37].This Hamiltonian supports q global minima, cf. the Q-function in Fig. 1 (left) for q = 6.Here, we have rescaled the phasespace coordinates such that the q global minima fulfill |x + ip| = 1 corresponding to different classical solutions.Remarkably, quantum fluctuations do not introduce any tunneling between these solutions as the corresponding coherent states |α m ⟩ = |e im 2π q / √ 2λ⟩ with m = 0, 1, • • • , q − 1 are exact zero-energy eigenstates.In other words, the groundstate manifold is q-dimensional space spanned by q q-legged cat states.
Note that since the Hamiltonian Q-function is a polynomial, its Fourier transform Eq. ( 6) is divergent.To solve this problem, we renormalize the divergence introducing the bounded Hamiltonian Ĥ(T ) We can calculate analytically f T (k x , k p ) and V (x, t) for Ĥ(T ) F γ for any arbitrary positive integer q and γ > 0. This allows us to arrive at a closed expression for the driving potential in the limit γ → 0 (see App. D) ] .We note that for q = 2 we recover a well-known result: Eq. ( 10) corresponds to a parametrically driven Duffing oscillator [53,54,57,58].We further note that the driving period is T d = T /q which directly follows from the q-fold rotational symmetry of the Floquet Hamiltonian.
Realizing Hamiltonian ( 9) is appealing in view of quantum computation because weak photon decay, with rate κ ≪ 2βq 2 , steers the oscillator towards its groundstate manifold containing code and error spaces of cat code [8], see Ref. [28] for q = 2 and Appendix for the general case.In the App.E, we numerically verify the quality of the weak dissipation and rotating wave approximation for realistic parameters.

B. Sharp-boundary well
Next, we demonstrate that our method allows us to engineer wells with a sharp boundary in phase space.For concreteness we choose an elliptical shape, i.e.H (T ) Q (x, p) = −β inside the white dashed line in Fig. 2(a) and H (T ) Q = 0 otherwise.In the classical limit λ → 0, our method allows us to find a closed-form solution for V (x, t) (see App. G).However, our solution is divergent at two time-dependent positions.In addition, it does not directly apply to the quantum regime, λ ̸ = 0, because the dependence of V (x, t) on λ is not analytical.This is due to the exponential factor in Eq. ( 6) leading to divergent NcFT coefficients f T (k x , k p ) in the limit of large wavevectors, k 2 x + k 2 p → ∞, for any λ ̸ = 0. We remove these unphysical features by smoothing out the target Floquet Hamiltonian by applying a convolution with a Gaussian kernel with standard deviation σ, cf.Fig. 2(a).For σ above a threshold, σ > λ/2, the NcFT spectrum f T (k x , k p ) becomes integrable and, thus, leads to a smooth solution for V (x, t), cf.Fig. 2(b) and the closed expression in the App.G.This implies that we can implement a potential step that is arbitrarily sharp compared to the typical dimensions of phase-space well, but should remain smooth on the scale of the oscillator quantum fluctuations.Note that Floquet Hamiltonians with sharper boundaries (σ < λ/2) are well-defined but cannot be realized using our method (see App. B and App.G).The spectrum and first few eigenmodes are also shown in Fig. 2(c).The latter are squeezed non-gaussian states.

C. Moiré superlattice
In Ref [41], we have shown how to synthesize arbitrary lattices in phase space.We can use our method to combine a lattice potential with a sharp confinement realizing a finite-size lattice.For concreteness we focus on a Moiré superlattice, cf.Fig. 3(a).This is the phase-space equivalent of the 2D potential for electrons in twisted graphene [59][60][61][62].The Moiré superlattice is formed by overlay- ing two honeycomb lattices with a relative twist angle θ 0 in a finite region of radius R. Outside of this region H (T ) F (x, p) = 0.The resulting Hamiltonian Q-function for the twist angle θ 0 = 10 • is shown in Fig. 3(a).As discussed above, overall the Floquet Hamiltonian should be smooth on the scale of the oscillator quantum fluctuations.As for the phase-space well example above, this can be implemented by applying a convolution with a Gaussian kernel to the initially discontinuous Floquet Hamiltonian.The ensuing transition between the Floquet lattice potential and the phase-space region with H (T ) F (x, p) = 0 can be arbitrarily sharp compared to R or to the honeycomb lattice constant.A closed formula for the Floquet Hamiltonian is given in the App.H.
Applying our method, we calculate the NcFT spectrum f T (k x , k p ) shown in Fig. 3(b).It is formed by three groups of twelve peaks.Each group of peaks is obtained from a single peak by applying one of the six-fold phasespace rotations and/or the rotation by the twist angle θ 0 , cf Fig. 3(a).The width of all the peaks is ∝ R −1 .All these features as well as the exact locations of the peaks can be read out from a closed-form solution for f T (k x , k p ) given in the App.H.In Fig. 3(c), we plot the ensuing driving potential V (x, t) for 0 ≤ t < T d .[In this case, the driving period T d is one-sixth of the natural period, T d = T /6, reflecting the 6-fold rotational-symmetry of our target Floquet Hamiltonian.]In Fig. 3(d), we also plot the instant driving potential at t = 0 and t = T d /2 (or t = T /12).We note that the real-space driving potential is a sequence of discrete lattice potentials localized in a finite region of real space that are switched on for a short time interval.We note further that in the limit R → ∞, the peaks in (k x , k p )-space become δ-functions, and the driving potential reduces to a discrete sequence of stroboscopic lattices with specific amplitudes, wavelengths, and phases [39,41,63,64].Considering that the contact interaction of cold atoms turns into a longdistance Coulomb-like interaction in the rotating frame [39,41,[64][65][66][67][68][69][70], many atoms in the phase space Moiré superlattice would mimic the behavior of electrons in twisted bilayer graphene [59][60][61][62].

D.
Artificial atomic spectrum Our method can be leveraged to implement a target spectrum {E n } as well as desired target eigenstates {|ψ n ⟩}.As mentioned above, this could be useful for quantum simulations with interacting atoms.In this scenario, our method could be straightforwardly applied to the target Floquet Hamiltonian Ĥ(T ) For concreteness, we consider |ψ n ⟩ = |n⟩ where {|n⟩} is the harmonic oscillator (Fock states) eigenbasis.In this example, the Hamiltonian Q-function and the NcFT spectrum can be easily expressed as a sum over the excitation number n, and respectively.Here, 1 F 1 (a; b; z) is the Kummer confluent hypergeometric function.The driving potential V (x) can be straightforwardly calculated by plugging Eq. ( 12) into Eqs.(7) and (8).Note that since the NcFT spectrum f T (k cos τ, k sin τ ) is independent of the angular coordinate τ , the driving potential V (x) is static.This, in turn, follows from our choice of eigenbasis leading to a target Floquet Hamiltonian invariant under arbitrary phase-space rotations, cf.Eq. (11).Note further that the asymptotic behavior for k → ∞ ensures that the integral in Eq. ( 7) is well defined.In Fig. 4 we display the potential V (x) for two interesting choices of the spectrum {E n }.In panel (a), we fix {E n } to be the spectrum of the hydrogen atom ) while all other levels are zero, E n>1 = 0. Thus, at λ = 3/4, the energies E 0 and E 1 of the second spectrum display an exact crossing.

VI. STATE PREPARATION
Our method combined with an adiabatic ramp protocol following Ref.[31] can be exploited to prepare a desired quantum state.As example, we demonstrate the preparation of a cat state in the groundstate manifold of Hamiltonian Eq. ( 9), including also the effects of dissipation, see App.F.

VII. EXPERIMENTAL IMPLEMENTATIONS
In order to design arbitrary Hamiltonians in phase space, one needs the ability to engineer the driving realspace potential V (x, t) in experiments.This might be difficult in practice.An alternative route is to directly use Eq.(7).In the App.F 2, we show that the target Floquet Hamiltonian can be well-approximated by replacing the integral with sum of a finite number of cosine lattice potentials.For example, we demonstrate the preparation of a three-legged cat state with 99% fidelity using only 5 such potentials.In cold atom experiments, the building block cosine lattice is formed by laser beams intersecting at an angle [41,71,72].In ex-periments with superconducting circuits [73][74][75], a microwave cavity in series with a Josephson junction (JJ) biased by a dc voltage (V ) is described by the Hamiltonian Ĥ , where E J is the JJ energy, ω J = 2eV /ℏ is the Josephson frequency and ∆ = 2e 2 /(ℏω 0 C) with C the cavity capacitance [76][77][78][79][80][81][82][83][84][85][86][87].

VIII. SUMMARY AND OUTLOOK
In this work, we have introduced a general constructive method to derive the driving potential, up to leading order in the Floquet-Magnus expansion, generating any arbitrary Floquet Hamiltonian of a single Bosonic mode.We have also shown that, in App.E and App.F, it can be transferred to state-of-the-art experimental platforms to efficiently prepare quantum states as part of a longlived quantum memory.A natural extension of our work would be to include higher-order perturbative corrections as the inverse problem of the Floquet-Magnus theory.Another exciting prospect is to extend our method to a many-body scenario by upgrading the single-particle plane-wave operator exp[i(k x x + k p p)] used in Eq. ( 5) to a many-body equivalent exp[ j i(k j x xj + k j p pj )].In experiments with superconducting circuits, this could be implemented coupling a dc-voltage biased JJ to multiple superconducting cavities [76,81,82,84,85,87].

Appendix A: Noncommutative Fourier transformation
In this section, we provide detailed calculation of the noncomutative Fourier transformation (NcFT) coefficient for a given target Floquet Hamiltonian operator Ĥ(T ) F = Ĥ(T ) F (x, p).We start from writing the target Hamiltonian as a sum of plane-wave operators, cf.Eq. ( 5) in the main text, In order to calculate the Fourier coefficient f T (k x , k p ), we first express the target Hamiltonian with reordered ladder operators Note that the ordering here keeps all the terms from commutators.By defining the coherent state |α⟩ as the eigenstate of lowering operator â|α⟩ = α|α⟩, we calculate the operator in the diagonal coherent representation Q (α, α * ) can also be written as In order to calculation the NcFT coefficient f T (k x , k p ) in Eq. (A1), we need to calculate the coherent diagonal element of the plane-wave operator ⟨α|e i(kx x+kp p) |α⟩.For this purpose, we introduce the displacement operator Dα ≡ e αâ † −α * â with the following relationship [69] Dα Dβ = e iIm(αβ * ) Dα+β , Dα |β⟩ = e iIm(αβ * ) |α + β⟩.
From the hermicity of Hamiltonian operator, we have the following important relationship Note out that here we present a general way to calculation the NcFT coefficient.In practice, for some specific target Hamiltonians, there may exist a simpler and more direct way to obtain the result as for the rotational lattice shown below.
Appendix B: One-to-one correspondence between Hamiltonian operator and Q-function In the above derivation of the NcFT coefficient for a given target Hamiltonian operator, we perform Fourier transformation of the Hamiltonian Q-function that only takes the diagonal elements of Hamiltonian operator in the coherent state representation, cf.Eq. (A2).One may wonder if some information is lost by neglecting the off-diagonal elements.In this section, we will prove the Hamiltonian operator Ĥ(T ) F (x, p), given by Eqs.(A1) and (A7), is fully determined by its Hamiltonian Q-function where f nm (k x , k p ) is the NcFT coefficient of the operator |n⟩⟨m| given by with the Q-function of operator |n⟩⟨m| given by Because the target Hamiltonian is the liner superposition of |n⟩⟨m| with n, m = 0, 1, • • • , we just need to prove By introducing x = r cos θ, p = r sin θ and k x = k cos τ, k p = k sin τ , we have the Fourier component where J n−m (z) is the Bessel function with order of n − m, and 1 F 1 (a; b; z) is the Kummer confluent hypergeometric function.We introduce the marix element [88] ⟨n where (z) is the generalized Laguerre polynomial.Then, we have the matrix element of H nm (x, p) in Fock representation This is the identity (B4) we aim to prove.As a result, for a fixed λ, the mapping between Floquet Hamiltonians and Q-functions is one to one.Note that the exponentially suppression factor e − λ Appendix C: Designing driving potential In this section, we show how to construct the driving potential from the NcFT coefficient such that its Floquet Hamiltonian equals to the target Hamiltonian in the leading order (RWA).We introduce the polar coordinate system in (k x , k p ) space via (k x = k cos τ, k p = k sin τ ), and write the Fourier expansion Eq. (A1) as Here, we have defined the Fourier component in the polar coordinate system via f T (k, τ ) ≡ f T (k x , k p ), and allow for negative k via the relation f T (k, τ ) ≡ f * T (−k, τ ), cf.Eq.(A8).From the Fourier component of the target Hamiltonian f T (k, τ ), we design the real space driving potential as follows where the phase variable ω 0 t in time domain plays the role of angle τ in (k x , k p ) space.By setting driving period T d = 2π/(ω 0 q), the Hamiltonian in the rotating frame, cf.Eq. ( 3) in the main text, becomes By time averaging the above Hamiltonian, cf.Eq. ( 4) in the main text, and comparing the averaged result to Eq. (C1), one can directly find that the lowest-order Floquet-Magnus expansion (RWA) gives the target Hamiltonian Ĥ(T ) F (x, p).We can also write the engineered driving potential in real space as where we have introduced phase ϕ(k, t) = Arg[f T (k, ω 0 t)] and used the property f (−k, ω 0 t) = f * (k, ω 0 t).Thus, the driving potential can be engineered by superposing a series of cosine lattice potentials with tunable amplitudes |kf T (k, ω 0 t)| and phases ϕ(k, t).

Appendix D: Rotational lattice Hamiltonian
We apply our method to engineer the target Floquet Hamiltonian with q-fold discrete rotational lattice symmetry in phase space where the factor e −γâ † â with γ > 0 is introduced to suppress the divergence of Hamiltonian in phase space.The above Hamiltonian is a generalised version of the rotational lattice Hamiltonian discussed in the main text, and it goes back to Eq. ( 9) by setting α 0 = 1/ √ 2λ.Using the identity we have the Hamiltonian Q-function as follows Here, we have defined the parameter In order to obtain the analytical expression for the NcFT coefficient of Hamiltonian Q-function, we transform into the polar coordinate system by introducing (x = r cos ϕ, p = r sin ϕ) and (k x = k cos τ, k p = k sin τ ).Plugging Eq. (D2) into Eq.(A7), we have Here, J q (•) is the Bessel function of q-th order and 1 F 1 (a; b; •) is the Kummer confluent hypergeometric function.
In order to obtain an analytical expression for the driving potential V (x, t), we introduce the following identities: Identity II : Identity III : Here, Γ(•) is the Gamma function and D(z) = e −z 2 z 0 e t 2 dt is the Dawson function.Identity Eq. (D5) is the special case of identity Eq. (D4) by setting n = 0.In identity Eq. (D6), 2 F 1 (a, b; c; z) is the hypergeometric function given by the path integral in the complex ζ-plane [89] 2 The above integral is valid for For |z| < 1, the hypergeometric function can be written by the power series with the Pochhammer symbol (x) k = Γ(x + k)/Γ(x).The analytical continuation of 2 F 1 (a, b; c; z) into the domain |z| > 1 can be realised from the following relationship [89] 2 The above relationship is valid for | arg(−z)| ≤ π − ϵ (0 < ϵ < π) and a − b / ∈ Z. From Eqs. (C2) and (D3)-(D6), we obtain the analytical expression for the designed driving potential for finite values of γ > 0 and arbitrary complex number α 0 as follows Next, we discuss how to calculate driving potential V γ (x, t) in the limit of γ → 0 (σ γ = 1/ √ 1 − e −2γ → +∞).Using Eq. (D8), (D9) and the Euler's reflection formula Γ(z)Γ(1 − z) = π/ sin(πz), we obtain the following series expansion of the hypergeometric function for |z| > 1 and q Note that the parameter m actually can take the whole real values, i.e., m ∈ R.Although the above expansion is not defined for m at negative integers m ∈ Z − , but the limit values do exist and can be defined as the values of the expansion for m ∈ Z − .By plugging the above series expansion Eq. (D11) and also the confluent hypergeometric ), we obtain the driving potential in the limit of γ → 0 (σ γ → +∞) In line (D12), only terms that satisfy k +m−q = 0 give nonzero contribution otherwise 1 Γ(m+k−q+1) 1 2σ 2 γ k+m−q = 0 in the limit of σ γ = +∞ (note that Gamma function Γ(z) = ∞ for nonpositive integer argument z ∈ Z − 0 ).In line (D14), only terms with even integer q and m ≤ q/2 give nonzero contribution.Furthermore, we emphasise that the driving potential V (x, t) is obtained from the RWA.In the rotating frame, the oscillating terms from m < q/2 cannot cancel the time-dependent parts given by terms that contain e ±iqτ in line (D13).Therefore, the only nontrivial contribution comes from the term with m = q/2 in line (D14).For the same reason, only nontrivial contribution comes from the term with m = (q − 1)/2 in line (D15).By neglecting these terms, we have the designed driving potential By taking the value of α 0 = 1/ √ 2λ, we obtain the driving potential given by Eq. (10) shown in the main text, i.e., B q,m λ q−m x 2m − C q cos(qω 0 t)x q + 1. (D18) Here, we have defined the coefficients .
Appendix E: Dissipative dynamics for the Bosonic code Floquet Hamiltonian In this section, we discuss the dissipative dynamics for the Floquet Hamiltonian Eq. ( 9) of the main text, which holds q-fold rotational symmetry in phase space.We prove that our arbitrary phase-space Hamiltonian engineering method can indeed allow us to design driving potential such that photon loss (dissipation) naturally leads the system state into the code subspace.Here, we focus on the experimentally relevant scenario occurring at a photon loss rate κ that is much smaller than the typical frequency of the oscillations in the rotating frame (set by driving strength β).We mention in passing that the special case of q = 2 (parametric oscillator) has been already discussed in Ref. [28].
Note that the Floquet Hamiltonian Ĥ(T )

F
given by Eq. ( 9) in the main text is the special case of the Floquet Hamiltonian Ĥ(T ) F γ given by Eq. (D1) in this Supplementary Material by taking parameters γ = 0 and α 0 = 1/ √ 2λ.In fact, the Floquet Hamiltonian Ĥ(T ) F is somewhat "unphysical" because it is divergent when the phase-space coordinates approach infinity.In contrast, the Floquet Hamiltonian Ĥ(T ) F γ is more relevant to the real experiments and more convenient to be studied by numerical simulations since the divergence is suppressed by the exponential factor e −γâ † â.For this reason, below, we discuss the dissipative dynamics for the more general Floquet Hamiltonian Ĥ(T ) F γ .However, we anticipate that the same basic physics will apply to any sufficiently small γ and, thus, also to the limiting case γ = 0, corresponding to the Floquet Hamiltonian Ĥ(T ) F of the main text.

1.
Groundstate manifold and engineered driving potential Before delving into the dissipative dynamics of the oscillator, we introduce a basis of rotationally invariant states for the "ground-state" manifold of Hamiltonian Ĥ(T ) F γ and give more details about the implementation of Ĥ(T ) F γ .We remind that the q-fold symmetry of the Hamiltonian Ĥ(T ) F γ in phase space is described by Rq Ĥ(T ) where Rq ≡ e iâ † â 2π q is the discrete rotational operator [20,37].According to the Bloch theorem extended in phase space [20,37], the eigentsates of q-fold rotational Hamiltonian can be written in form of where index l labels the Bloch bands, m is called quasinumber representing the parity of state, |ϕ l ⟩ is the cell state of l-th Bloch band and N l,m is the normalized factor.For the Hamiltonian Ĥ(T ) F γ given by Eq. (D1), the q standard coherent states |α 0 e γ+i 2πm q ⟩ with m = 0, 1, • • • , q − 1 are the degenerate exact zero-energy eigenstates.We thus choose the coherent state |α 0 e γ ⟩ as the cell state for the lowest band, |ϕ 0 ⟩ = α 0 e γ ⟩, and construct the Bloch states of lowest band Here, we have omitted the band index l = 0 for simplicity.The above q q-legged cat states |ψ m ⟩ construct the code subspace.For example for q = 3, the code space is panned by the following three code states Next, we give more details about the implementation of Ĥ(T ) F γ .Also in this case we take as an example for concrete numerical results the case q = 3, but our more qualitative discussion will apply equally well to other integer values of q.In Fig. 5(a), we calculate and plot the corresponding Hamiltonian Q-function H (T ) Qγ (x, p)/β, cf.Eq. (D2), for given parameters.The analytical expression for the driving potential V γ (x, ω 0 t) that generates the target Hamiltonian Ĥ(T ) F γ is given by Eq. (D10).In Fig. 5(b), we plot V γ (x, ω 0 t) in one time period T = 2π/ω 0 .Different from the power-law driving potential given by Eq. ( 10) in the main text, the driving potential here for generating Ĥ(T ) F γ is confined in a finite region in real space due to the suppression factor γ. As discussed in the main text, the driving potential V γ (x, ω 0 t) can be formally written as a superposition of cosine potentials Here, the time-varying amplitudes A(k, ω 0 t) and phases ϕ(k, ω 0 t) are determined from the non-commutative Fourier coefficients f T (k, ω 0 t) given by Eq. (D3) in polar coordinates (k x = k cos ω 0 t, k p = k sin ω 0 t) In Fig. 5(c) and (d), we plot the time-varying amplitude A(k, ω 0 t) and the phase ϕ(k, ω 0 t) respectively as functions of time 0 ≤ t ≤ 2π/ω 0 and wavenumber 0 ≤ k ≤ 10k c , where we have defined characteristic wavenumber 2.
Dissipative dynamics leading to the groundstate manifold Here, we go back to the main goal of this section, i.e. to study the dissipative dynamics of the oscillator when the Floquet Hamiltonian Ĥ(T ) F γ is prepared using our method.In the presence of weak photon loss and pure dephasing, the dissipative dynamics of the density matrix ρ(t) is described by the Lindblad master equation [90] where { Â, B} ≡ Â B + B Â is the anticommutator, κ is the single-photon loss rate and η is the dephasing rate.Here, Ĥ(t) is the full Hamiltonian with the potential V γ (x, t) as derived using our method, cf Eq. (E7).We note in passing that the master equation (E10) is valid for any weakly non-linear high-quality-factor oscillator [90].In our work, these two conditions translate into β ≪ ω 0 , and κ, η ≪ ω 0 , respectively.Below we show that in the realistic parameter regime the oscillator tends to relax into the groundstate manifold of the Bosonic code Floquet Hamiltonian.We will corroborate our analytical derivation with numerical results obtained by directly solving Eq. (E11), cf Fig. 2. For our initial analytical treatment, we consider the dissipative dynamics when photon decay is the only decay channel (setting the dephasing rate η = 0) and switch to the rotating frame arriving at the Lindblad master equation [8,18,22,28,53,[90][91][92][93] Going from the lab-frame master equation (E10) to the rotating-frame master equation (E13), we have further simplified the description by applying the RWA to the Hamiltonian in the rotating frame.This approximation [45,46] is standard and is consistent with the assumption β ≪ ω 0 leading to Eq. ( 47).We recall that the Floquet Hamiltonian Ĥ(T ) F γ comprises q separate wells.In the limit λ ≪ 1, the transition between different wells occurs on a very large time scale (exponentially large in λ −1 ).This means the total process towards equilibrium via photon loss can be divided into two distinguished stages: a fast process of time scale κ −1 relaxing to local equilibrium point (x e , p e ) in each local well and a slow transition (tunneling) process between local wells.The total probability in the code space is governed by the first quench process because the second slow transition process only adjusts the distribution over the states inside the code space.This makes it convenient to first analyze the linearized dynamics within one well.Below we follow the general treatment of Ref. [94].
For concreteness, we consider the well about α e = (x e + ip e )/ √ 2λ, cf. the minima of the Hamiltonian in Fig. 5(a).Up to leading order in δâ ≡ â − α e , we can approximate the local rotating wave Hamiltonian as a harmonic oscillator Thus, the oscillations about each well have frequency 2βq 2 , giving rise to the quantized spectrum E n = λ2βq 2 (n+1/2).We note that the broadening of the Floquet spectrum due to photon decay is of the order λκ.much smaller than the typical level spacing, κ ≪ 2q 2 β, the dissipative dynamics is well approximated by a rate equation for the Floquet states [53,94] Here, P l = ⟨ϕ l |ρ|ϕ l ⟩ is the occupation of the local Floquet states |ϕ l ⟩, i.e., the cell state of l-th Bloch band defined in Eq. (E2).The transition rates W ll ′ according to Fermi golden rule are given by Because the local ground state for the quasienergy well is simply the standard coherent state |ϕ 0 ⟩ = |α e ⟩ according to our target setting, the transition rates from a local ground to the excited states are exactly zero, i.e., W l0 = 0 for l > 0. Thus, photon decay induces transitions only towards the bottom of the well.Obviously, the same analysis applies to all q quasienergy wells.Thus, any superposition of states each localized about different quasi-energy wells will relax towards the groundstate manifold spanned by the coherent states |α m ⟩ = |e im 2π q α e ⟩ with m = 0, . . .q − 1.We note in passing that Eq. (E15) can be easily modified to account for finite photon dephasing with rate η corresponding to the RWA Lindblad master equation [8,18,22,28,53,92,93] In the limit 2q 2 β ≫ κ, η/(2λ), we arrive at Eq. (E15) with modified transition rates By making harmonic approximation near the bottom of quasienergy well, the local Floquet levels are simple displaced Fock states |ϕ l ⟩ = δâ †l |α e ⟩/ √ l!.The modified transition rates are approximately leading to a steady state Boltzmann distribution over the quasienergy states with the effective thermal occupation number from rate equation (E15) This corresponds to the steady state groundstate manifold occupation probability If the dephasing rate is weak enough compared to the photon loss rate satisfying η|α e | 2 /κ ≪ 1, the leakage probability is also small (1−P 0 ≈ η|α e | 2 /κ).However, if the dephasing rate is strong so that η|α e | 2 /κ ≫ 1, the leakage probability is large (P 0 ≈ 0) indicating that the state preparation via photon loss does not work anymore.Summarizing our analysis so far, we can conclude that the oscillator tends to relax into the groundstate manifold of the Bosonic code Floquet Hamiltonian in the parameter regime identified by the set of conditions shown in Eq. (E12).The first condition ensures that photon decay is the dominant decay channel, the second that the broadening of the Floquet levels is small, and the third that the RWA is valid.This is a realistic parameter regime in experiments with superconducting circuits.
This conclusion is also corroborated by numerical results, shown in Fig. 2, obtained by simulating the full dissipative dynamics as defined by the master equation (E10), which includes the full time-dependence of the coherent Hamiltonian Eq. (E11).In these simulations, for the case q = 3, we start from the zero-excitation Fock state and choose the driving strength β = 1/40.In Fig. 6(left), we show the time evolution of total probability of density matrix over the three-fold rotational code states given by Eq. (E4), i.e., for different photon loss rates at zero dephasing rate η/ω 0 = 0.The probability P code (t) increases monotonously from an initial small value P code (0) ≈ 0.1 with faster speed for stronger photon loss.P code (t) is larger as the photon loss rate becomes stronger till an optimal value κ/ω 0 = 2.0 × 10 −3 (up to P code > 0.95), then drops again as photon loss rate continues increasing.This is because when the photon loss rate becomes large enough so that the condition κ ≪ 2βq 2 ω 0 ≈ 0.4ω 0 is not satisfied, the off-diagonal matrix elements of density matrix play the role and the validity of rate equation (E15) breaks down.In Fig. 6(right), we show the time evolution of P code (t) for different photon loss rates with finite dephasing rate η/ω 0 = 10 −5 .The results are qualitatively the same as that in Fig. 6(left) but with a slightly optimal point of photon loss rate and a lower probability in code space due to the effective temperature from pure dephasing, cf.Eq. (F8).

Dissipative dynamics within the code and error spaces
Next, we briefly discuss the dynamics within the groundstate manifold of the Floquet Hamiltonian (D1).As we explain below this manifold can be chosen to host the code and error spaces for a cat code.Here, we assume that the conditions in Eq. (E12) are fulfilled and go back to the simplified description without counter-rotating terms and dephasing.It is convenient to rewrite the Lindblad master equation projected onto the groundstate manifold into a basis of q-legged cat states, cf.Eq. (E3), Here, N m are normalization constants and the quantum number m can be interpreted as the quasi-angularmomentum or, equivalently, the photon number modolus q.We note that in the semi-classical limit, λ → 0, all normalization converge exponentially fast to the same value, lim λ→0 N m = q.We note further that the annihilation operator increases the value of the quasi-angularmomentum, Thus, the annihilation operator â projected onto the groundstate manifold reads where Π 0 is the projector on the groundstate manifold.If λ is small enough such that the last approximation is accurate, the dissipative dinamycs of density matrix ρ code within the groundstate manifold is described by the simple master Lindblad Master equation with the simple jump operator L = q−1 m=0 |ψ m+1 ⟩⟨ψ m | and damping rate κ ′ = κ|α 0 e γ | 2 .The groundstate manifold can be chosen to host the code and error spaces for a bosonic code if q = 2(N + 1).In this case, the code space is spanned by the cat states with quasi-angular momentun m = 0 and m = N + 1.The remaining cat states span the error space.Such a code allows to correct the simultaneous decay of up to N photons.This can be straightforwardly verified by applying the corresponding Knill-Laflamme conditions [95].
The standard cat code (correcting single photon decay) corresponds to the case q = 4 (or N = 1) [22].In this case, the syndrome is a simple parity measurement (code states have even parity) [96].When every error is detected by applying repeated parity measurements, they can also be corrected by updating the definition of the code and error states using the mapping m → m + 1, e.g. after the first error is detected the odd (even) states become the code (error) states.

Adiabatic ramp
To prepare the target bosonic code states, we adopt the recently proposed "adiabatic ramp" method introduced in Ref. [31].Following this method, we set the system Hamiltonian with periodic driving potential as where V γ takes the form given by Eq. (E7) with time-dependent driving amplitude β(t) and frequency ω(t).The main idea is that the driving potential in Eq. (F1) is turned on adiabatically from β(0) = 0, ω(0) ̸ = ω 0 within a preparation time t pre to the values of β(t pre ) = β 0 , ω(t pre ) = ω 0 .In such a way, the initial cavity state is adiabatically ramped to the target code states.Following the ramp receipt present in Ref. [31], we modulate the driving amplitude β(t) and frequency ω(t) in form of the sigmoidal function with the center of ramp t c , the slope of ramp s.The parameters C and D are determined by the boundary conditions As in Ref. [31], we set g(0) = β(0) = 0, g(t pre ) = β(t pre ) = 0.025ω 0 , t c = t pre /6 and s = 60/t pre for getting the profile of driving amplitude β(t).Note that the Floquet adiabaticity may be lost when two Floquet states have equal quasienergies modulo ω 0 = 2π/T [31].This problem can be circumvented by making the driving frequency ω(t) incommensurate with but close to the harmonic frequency [31].Here, we set g(0) = ω(0) = ω 0 /(1 + π × 10 −3 ), g(t pre ) = ω 0 , t c = 2t pre /3 and s = 30/t pre to fix the profile of driving frequency ω(t).
In the lower panel of Fig. 7(a), we plot the time-varying driving amplitude β(t) and frequency ω(t) for a preparation time t pre = 5000T .The preparation process starts from an initial cavity state |ψ pre (0)⟩, and the prepared state |ψ pre (t)⟩ is obtained by Schroedinger equation Note that the three code states given by Eqs.(E4)-(E6) are three different eigenstates of phase-space rotational operator Rq ≡ e iâ † â 2π q with different parity m, which is kept unchanged during the adiabatic preparation process by our designed driving potential.According to Eq. (E2), the code states |ψ 0 ⟩, |ψ 1 ⟩ and |ψ 2 ⟩ can be adiabatically achieved from Fock states |0⟩, |2⟩ and |1⟩ respectively.As the cavity vacuum state is typically easier to start with than other Fock states, we focus on preparing the code state |ψ 0 ⟩ given by Eq. (E4) below.
To show the quality of prepared state, we define the fidelity F [ρ 0 , ρ pre (t)] of prepared state ρ pre (t) with respect to the target code state ρ 0 = |ψ 0 ⟩⟨ψ 0 | by [97] In the middle panel of Fig. 7(a), we plot the time evolution of fidelity calculated at the stroboscopic time moments (ω 0 t/2π ∈ Z + ).It clearly shows that the prepared state starts with a low fidelity (F < 0.4) and approaches the target state with high fidelity (F > 0.99).In the upper panel of Fig. 7(a), we plot several snapshots of Husmi Q-functions of prepared states.From an initial vacuum state of cavity, the prepared state begins transiting to the target state when the driving amplitude β(t) starts to ramp up around t = 700T .Thereafter, the prepared state actually already achieves the target state with high fidelity but keeps oscillating due to the finite detuning (ω(t) ̸ = ω 0 ).When the driving amplitude β(t) is close to the final value 0.025ω 0 , the driving frequency ω(t) starts to ramp up and the fidelity of prepared state becomes lower (see the snapshot at t = 3400T ).Finally, when the driving frequency ω(t) approaches the value of ω 0 , the prepared state is stabilized to the target state with high fidelity again.
According to the Floquet adiabatic condition given in Ref. [31], the quality of prepared state becomes better when the preparation time is longer.In Fig. 7(b), we plot the fidelity of final prepared state F [ρ 0 , ρ pre (t pre )] as a function of preparation time t pre .The plot verifies that the quality of prepared state is already very good when t pre > 3000T .For the typical cavity frequency ω 0 /2π = 5GHz in the superconducting circuits, the preparation time is less than 1µs which is much faster than the preparation time using other protocols in circuit QED [13,98,99].

Discretization of wavenumbers
We now discuss another possible errors from implementation of our driving protocol.As given by Eq. ( E7), the engineered driving potential is written in a superposition of cosine potentials with modulated amplitudes and phases in time.In the real experiments, such driving potential can be created by laser beams for cold atoms [41] or by superconducting circuits with Josephson junctions [31].In both experiments, the potential is approximated by finite number of cosine functions with discretized wavenumbers (F5) In our numerical simulation, we choose the cutoff of wave number k cut = 10k c , where the characteristic wavenumber is defined by k c ≡ (1 − e −2γ )/4λ, cf.Eq. (E9).
In Fig. 8, we plot the time evolution of fidelity for different choice of discretized wavenumber step ∆k = k c /10, ∆k = k c and ∆k = 2k c corresponding to 100, 10 and 5 tunable Josephson junctions (optical lattices) in the superconducting circuits experiment (cold atoms experiment) respectively.This means the number of Josephson junctions (laser beams) in the circuit-QED (cold atom) experiments and thus the complexity of total operations can be significantly reduced.
Our numerical results show that the discretization of the wavenumbers causes some discrepancies during the initial phase of the preparation and also small oscillations in the long-time behavior.Nevertheless, the averaged final fidelity of the prepared state can reach up to 99.64%, 99.63% and 99.06% for ∆k = 0.1k c , ∆k = k c and∆k = 2k c respectively.The fidelity of the state prepared by using our scheme is higher than the fidelity obtained for cat states of similar amplitude using a sequence of interleaved Selective Number-dependent Arbitrary Phase (SNAP) gates and displacement gates [27].We note in passing that the infidelity obtained using our method could be further increased without adding any additional experimental complexity.In fact, we expect the residual infidelity in our numerical results to be mainly due to high-order Floquet-Mafnus corrections (beyond the RWA) creating a deviation between the implemented and the target Floquet Hamiltonian.This deviation could be reduced by fine-tuning the driving potential to account for higher-order terms.Such an extension of our work is in preparation.
Coherently controlling multiple tunable Josephson Junctions (JJs) for designing functional quantum devices and quantum computation/simulation is a well-established technology in circuit-QED architectures.Examples include the Josephson ring modulator architecture [100,101] with 4 JJ (one for each transmon qubit), the quantum-statepreservation superconducting circuit [102] with 9 transmons, the Google programmable superconducting processor Sycamore [103] with 54 transmon qubits and the recent IBM quantum processor Eagle [104] with 127 transmons qubit.We also note in passing that, in the spirit of Trotter discretization [105], our driving scheme could even be realized even with a single transmon by decomposing the multiple JJs unitary operation into a sequence of discrete gate operations.We leave a detailed study of this scenario to a future work [106].

Photon loss
We first set pure dephasing rate η = 0 and discuss the effects of finite photon loss rate κ > 0. To analyse the errors from photon loss during state preparation, we unravel the master equation Eq. (E10) with wave-function Monte Carlo method [107] in the framework of quantum trajectory theory [108].For an initial pure quantum state, each stochastic quantum trajectory of prepared state |ψ together with a normalization procedure of |ψ   pre (t)⟩ and dW n is Wiener noise increment satisfying (dW n ) 2 = dt independent of each quantum trajectory.The density matrix of the prepared state at time t is approximated by N quantum trajectories The first term on the right hand side (RHS) of Eq. (F6) represents the unitary time evolution by adiabatic ramp.The second and third terms on the RHS of Eq. (F6) describe dissipation and fluctuation (quantum jump) respectively.
The quantum trajectory Eq. ( F6) is used to analyze the effects of noises during preparation.To verify the above analysis, we obtain numerical results directly from master equation (E10).We define the probability of the prepared state over the target state and compare it to the total probability in the code space P code (t) defined by Eq. (E22).In Fig. 9(a), we plot P 0 (t pre ) and P code (t pre ) of prepared state as function of preparation time t pre with finite photon loss rate κ/ω 0 = 10 −6 .As expected, although the probability P 0 (t pre ) of prepared state over the target state at final time t pre = 5000 × 2π/ω 0 is indeed lowered a bit, the total probability in the code space P code (t pre ) for long reparation time t pre = 5000 × 2π/ω 0 is close to one.In Fig. 9(d), we plot P 0 (t pre ) and P code (t pre ) of prepared state at final time t pre = 5000 × 2π/ω 0 as function of photon loss rate κ/ω 0 .Our result shows that the quantum information is well protected inside the code space, i.e., P code (t pre ) keeps close to one, as long as κ/ω 0 < 10 −5 which is satidfied in the real circuit-QED experiments [98,99].For strong photon loss rate κ/ω 0 > 10 −5 , the quantum information will start to leak outside the code space due to the nonnegligible excitation of prepared state into high-band states (P code < 1).The main channel is â † â|ψ l,m ⟩ in the dissipation term on the RHS of Eq. (F6), i.e., â † â|ψ l,m ⟩ ∝ q−1 p=0 e imp 2π q â † â Rp q |ϕ l ⟩ ∝ q−1 p=0 e imp 2π q Rp q â † â|ϕ l ⟩. (F9) In the second step, we have used the property R † q â † â Rq = â † â.Because â † â|α⟩ = αâ † |α⟩ is no longer a standard coherent state, the resultant state â † â|ψ l,m ⟩ can be a superposition of the states with the same parity m in all Bloch bands.Thus, the operator â † â can introduce inter-band transition making the quantum information leaking outside of the code space.Other channels introducing such leakage errors may include the non-RWA effects, the adiabatic approximation and the deformation of the cell state |ϕ l ⟩ from the standard coherent state due to the finite detuning (ω(t) ̸ = ω 0 ), cf. the snapshot in the upper panel of Fig. 7(a) at earlier time of preparation process t = 700 × 2π/ω 0 .
The photon loss leads the system state to the code subspaces without distinguishing the states in the code space.To prepare a specific target state in code space, one can detect the photon number loss and track the change of parity during preparation.As single-photon loss only alternates the Bloch index m → m + 1 during the whole preparation time according to Eq. (F7), we can correct such flipping errors inside the code space by updating the knowledge of quantum state when one single photon is detected without backaction on the encoded system.This error correction scheme can be extended to mult-photon loss because of ân |ψ l,m ⟩ ∝ |ψ l,m+n ⟩.All in all, we only need to counter how many photons are detected in total, and update our knowledge accordingly in the end end of preparation process.As long as the efficiency of photon number detection is perfect, the errors from photon loss can be 100% tracked and corrected.Alternatively, one can directly perform parity measurement and post select the code state [22].

5.
Pure dephasing Now, we set photon loss rate κ = 0 and discuss the effects of finite pure dephasing rate η > 0. In this case, the stochastic quantum trajectory of prepared state |ψ )+tanh(α 2 Driving    pre (t)⟩.Accordding to Eq. (F9), starting from a rotational state |ψ l,m ⟩, the noisy terms â † â|ψ l,m ⟩ and (â † â) 2 |ψ l,m ⟩ on the RHS of Eq. (F10) keep the parity m of prepared state exactly unchanged during the whole preparation process.As a result, the probability of the prepared state over the target state P 0 and the total probability in the code space P code must be exactly the same (no state flipping in the code space).
In Fig. 9(b), we plot P 0 (t pre ) and P code (t pre ) of prepared state as a function of preparation time t pre with finite pure dephasing rate η/ω 0 = 10 −7 and zero photon loss κ/ω 0 = 0.It clearly shows that the probability of prepared state over the target state P 0 (t pre ) coincides with total probability in the code space P code (t pre ) for any preparation time.In Fig. 9(e), we further plot P 0 (t pre ) and P code (t pre ) of prepared state at final time t pre = 5000 × 2π/ω 0 as function of pure dephasing rate η/ω 0 ∈ (10 −8 , 10 −4 ).The probabilities P 0 (t pre ) and P code (t pre ) still coincide with each other showing that the pure dephasing noise does not introduce intra-band transition.
Our code preparation protocol is robust as long as η/ω 0 < 10 −7 corresponding to pure dephasing time T ϕ < 0.32ms for a microwave superconducting cavity with frequency ω 0 /2π = 5GHz.The significant leakage of quantum information outside the code space for η/ω 0 > 10 −7 comes from the consequence of â † â|ψ l,m ⟩ and (â † â) 2 |ψ l,m ⟩ in the dissipation term on the RHS of Eq. (F10) as discussed below Eq. (F9).Finally, in Fig. 9(c), we plot P 0 (t pre ) and P code (t pre ) of prepared state as a function of preparation time t pre with finite photon loss rate κ/ω 0 = 10 −6 and finite pure dephasing rate η/ω 0 = 10 −7 showing that our code preparation protocol is indeed robust against decoherence in the noisy environment.the driving potential V (x, t) is the convolution of −x −2 and u(x, t), i.e., V (x, t) = −x −2 * v(x, t) Here, we should first replace the convolution function 1/z 2 by (z 2 − ϵ 2 )/(z 2 + ϵ 2 ) 2 to get converged integral, and then take the limit of ϵ → 0 to obtain V (x, t).
In the (x, p) phase space, we set a new coordinate (x ′ , p ′ ) system which is rotated by an angle τ between x and x ′ axes given by the following orthogonal transformation x = x ′ cos τ − p ′ sin τ ; p = x ′ sin τ + p ′ cos τ. (G4) We can express the target Hamiltonian Q-function with rotated coordinates by H ′ (T ) Q (x, p).Then, we project the Hamiltonian on the x ′ axis by the so-called Radon transformation The 1D Fourier transformation of the above projected function is given by The 2D Fourier transformation of H ′ (T ) Q (x ′ , p ′ )e −i(k x ′ x ′ +k p ′ p ′ ) dx ′ dp ′ .(G7) By plugging Eq. (G5) into Eq.(G6) and comparing with Eq. (G7), we have This is the so-called projection-slice theorem [110,111].Comparing to Eq. (A7) and using the orthogonal transformation (G4), we have
B and App.I].