# Engineering Arbitrary Hamiltonians in Phase Space

Lingzhen Guo  Max Planck Institute for the Science of Light, Staudtstrasse 2,
91058 Erlangen, Germany  Vittorio Peano  Max Planck Institute for the Science
of Light, Staudtstrasse 2, 91058 Erlangen, Germany

###### Abstract

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.

##  I Introduction

Generation of nonclassical bosonic states Gerry and Knight (2004); Strekalov
and Leuchs (2019); Kubala et al. (2015a), e.g., squeezed lights, Fock states
and Schrödinger’s cat states, is important not only for fundamental studies of
quantum mechanics but also for applications in quantum technologies Braunstein
and van Loock (2005); Pan et al. (2012); Strekalov and Leuchs (2019); Yan et
al. (2021). For example, bosonic states with discrete translational or
rotational symmetries in phase space Leghtas et al. (2013); Mirrahimi et al.
(2014); Heeres et al. (2017); Rosenblum et al. (2018); Flühmann et al. (2019);
Hu et al. (2019); Campagne-Ibarcq et al. (2020); Gertler et al. (2021) have
been proposed to encode quantum information Cochrane et al. (1999); Gottesman
et al. (2001); Travaglione and Milburn (2002); Michael et al. (2016); Albert
et al. (2018); Grimsmo et al. (2020), paving the way for hardware efficient
quantum error correction Tzitrin et al. (2020); Terhal et al. (2020); Joshi et
al. (2021); Cai et al. (2021). Bosonic code states can be prepared and
stabilized against dissipation via a sequence of universal gates, e.g.
interleaved selective number-dependent arbitrary phase (SNAP) gates and
displacement gates Krastanov et al. (2015); Fösel et al. (2020); Kudra et al.
(2022). A fundamental limitation of this approach is that the state can leak
out of the code and error subspaces Gottesman et al. (2001); Grimsmo et al.
(2020). Hamiltonian engineering can solve this problem, allowing to design
systems in which photon decay naturally leads into these subspaces Gottesman
et al. (2001); Puri et al. (2019); Rymarz et al. (2021); Conrad (2021). A
further advantage of this approach is that the dephasing error for the ensuing
logical qubit can be exponentially suppressed, allowing for fault tolerant
syndrome detection Puri et al. (2019).

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
phase-space loop acquires a geometric phase analogous to the Ahronov-Bohn
phase for a particles in a magnetic field. As a consequence, a gapped lattice
Hamiltonian in phase space can support non-trivial Chern numbers Zaslavskii et
al. (1986); Berman et al. (1991); Leboeuf et al. (1990, 1992); Gottesman et
al. (2001); Billam and Gardiner (2009); Guo et al. (2013); Zhang et al.
(2017); Liang et al. (2018); Lörch et al. (2019). 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 Guo et al. (2022), 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 Hamiltonian, known as
_Floquet Hamiltonian_ $\hat{H}_{F}$ which is defined via

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

Here, $\hat{U}(T,0)$ is the time-evolution operator from $t=0$ to $t=T$ where
$T$ is the time-period of the system’s time-dependent Hamiltonian
$\hat{H}(t)$. In adddition, $\lambda$ is an effective dimensionless Planck
constant, and $\mathcal{T}$ is the time-ordering operator. For a single
Bosonic mode or, equivalently, quantum mechanical particle, the Floquet
Hamiltonian $H_{F}(\hat{x},\hat{p})$ can be any arbitrary function in phase
space. Except for very few models, it is impossible to obtain a closed form of
the Floquet Hamiltonian $H_{F}(\hat{x},\hat{p})$ from the time-dependent
Hamiltonian $\hat{H}(t)$. Instead, one often evaluates the Floquet Hamiltonian
relying on a high-frequency expansion Rahav et al. (2003); Goldman and
Dalibard (2014); Bukov et al. (2015a), e.g. the Magnus expansion theory Casas
et al. (2001); Blanes et al. (2009), the van Vleck degenerate perturbation
theory Eckardt and Anisimovas (2015) and the Brillouin-Wigner perturbation
theory Mikami et al. (2016). In this work, we focus on the inverse problem,
that is, to find the time-dependent Hamiltonian $\hat{H}(t)$ that synthetizes
a target Floquet Hamiltonian $H^{(T)}_{F}(\hat{x},\hat{p})$. This is the realm
of Floquet engineering which is a very developed and active field Bukov et al.
(2015b); Liang et al. (2018); Rudner and Lindner (2020). 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

$\displaystyle\hat{\mathcal{H}}(t)=\frac{\omega_{0}}{2}(\hat{p}^{2}+\hat{x}^{2})+\beta
V(\hat{x},t).$ (2)

Here, $\omega_{0}$ is the oscillator natural frequency, $\beta$ is the
amplitude of the nonlinear driving potential $V(\hat{x},t)$ which has time-
period $T_{d}$ and might contain also static terms. In order to introduce an
effective dimensionless Planck constant $\lambda$, the position $\hat{x}$, the
momentum $\hat{p}$ and $\hat{H}_{\rm LF}(t)$ have been rescaled such that
$[\hat{x},\hat{p}]=i\lambda$ and at the same time the Schrödinger equation
reads $i\lambda\dot{\psi}=\hat{\mathcal{H}}(t)\psi$. The Floquet Hamiltonian
we want to design will be defined in a frame rotating with the oscillator
natural period $T=2\pi/\omega_{0}$ via the transformation
$\hat{O}(t)=\exp(i\hat{a}^{\dagger}\hat{a}\omega_{0}t)$, where $\hat{a}$ is
the oscillator annihilation operator,
$\hat{a}=(\hat{x}+i\hat{p})/\sqrt{2\lambda}$. We assume to be in the regime of
weak non-linearity, $\beta\ll\omega_{0}$, and that the natural period $T$ is
an integer multiple of the driving period $T_{d}$, $T=qT_{d}$. Then, the
Hamiltonian in the rotating frame
$\hat{H}(t)=\hat{O}(t)\hat{\mathcal{H}}(t)\hat{O}^{\dagger}(t)-i\lambda
O(t)\dot{O}^{\dagger}(t)$ reads

$\displaystyle\hat{H}(t)=\beta
V[\hat{x}\cos(\omega_{0}t)+\hat{p}\sin(\omega_{0}t),t].$ (3)

This Hamiltonian has time-period $T$ allowing to use Eq. (1) to define the
Floquet Hamiltonian. We note that a small detuning from the multiphoton
resonance can be incorporated in the driving potential $V(x,t)$. In addition,
the system in this frame evolves on the slow time-scale $\beta^{-1}$. Thus, we
are in the realm of application of the Floquet high-frequency expansions,
here, with the small parameter $\beta/\omega_{0}$. This allows to approximate
the Floquet Hamiltonian with the leading-order of the Floquet-Magnus expansion
corresponding to the rotating wave approximation (RWA),

$\displaystyle\lim_{\omega_{0}\to\infty}H_{F}(\hat{x},\hat{p})$
$\displaystyle=$ $\displaystyle\frac{1}{T}\int_{0}^{T}dt\hat{H}(t).$ (4)

Our goal is to engineer an arbitrary target Floquet Hamiltonian
$H^{(T)}_{F}(\hat{x},\hat{p})$ in phase space by properly designing the
driving potential $V(x,t)$ in real space, see Fig. 1. Up to leading order (or
within the RWA) we, thus, require that the right-hand side of Eq. (4)
coincides with the target Hamiltonian $H^{(T)}_{F}(\hat{x},\hat{p})$. The
ensuing solution becomes exact in the high-frequency limit
$\beta/\omega_{0}\to 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
$H^{(T)}_{F}(\hat{x},\hat{p})$ in the form of a _noncommutative Fourier
transformation_ (NcFT). This can be viewed as a variant of quantum
distribution theory Scully and Zubairy (1997). We wish to decompose the target
Hamiltonian $H^{(T)}_{F}(\hat{x},\hat{p})$ as a sum of plane-wave operators

$\displaystyle H^{(T)}_{F}(\hat{x},\hat{p})=\frac{\beta}{2\pi}\iint
dk_{x}dk_{p}f_{T}(k_{x},k_{p})e^{i(k_{x}\hat{x}+k_{p}\hat{p})}.$ (5)

It can be shown that the Fourier coefficients $f_{T}(k_{x},k_{p})$ are given
by the inverse transformation [see the Appendix A]

$\displaystyle
f_{T}(k_{x},k_{p})=\frac{e^{\frac{\lambda}{4}(k^{2}_{x}+k^{2}_{p})}}{2\pi\beta}\iint
dxdpH^{(T)}_{Q}(x,p)e^{-i(k_{x}x+k_{p}p)},\ $ (6)

where the phase-space function $H^{(T)}_{Q}(x,p)$ is the equivalent of the
Husimi Q-function, here, for an 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
$H^{(T)}_{Q}(x,p)=\langle\alpha|\hat{H}^{(T)}_{F}|\alpha\rangle$ for the
corresponding coherent state
$\hat{a}|\alpha\rangle=(x+ip)|\alpha\rangle/\sqrt{2\lambda}$. The latter mean
value can be calculated by normal ordering the target Hamiltonian
$\hat{H}^{(T)}_{F}(\hat{a}^{\dagger},\hat{a})$. We point out three important
features of the Hamiltonian $Q$-function: (i) For fixed $\lambda$, the 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) $\lim_{\lambda\to
0}H^{(T)}_{Q}(x,p)=H^{(T)}_{F}(x,p)$. For more details see the Appendix B and
Appendix H.

Figure 1:  Rotational lattice Hamiltonian in phase space. (Left) Q-function
$H^{(T)}_{Q}(x,p)$ for the target Floquet Hamiltonian Eq. (9) and (Right) the
engineered real-space potential $V(x,t)$ for parameters $q=6$, $\lambda=0.01$.
The white curves in the right figure indicate the locations of minimum
(stable) points of $V(x,t)$ at a fixed time moment.

##  IV Designing driving potential

The driving potential $V(x,t)$ that generates the target Floquet Hamiltonian
$H^{(T)}_{F}(\hat{x},\hat{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

$\displaystyle V(x,t)$ $\displaystyle=$
$\displaystyle\int_{0}^{+\infty}A(k,\omega_{0}t)\cos[kx+\phi(k,\omega_{0}t)]dk$
(7)

with time-varying amplitudes $A(k,\tau)$ and phases $\phi(k,\tau)$ determined
from the Fourier coefficients in polar coordinates ($k_{x}=k\cos\tau$,
$k_{p}=k\sin\tau$)

$\displaystyle A=k|f_{T}(k\cos\tau,k\sin\tau)|,\
\phi=\mathrm{Arg}f_{T}(k\cos\tau,k\sin\tau).\ $ (8)

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). For more details see the Appendix C. In the reminder of this
paper, we demonstrate the flexibility of our method 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
subtleties associated with it.

##  V Examples

###  V.1 Rotational lattice

We now apply our method to synthetize a particularly interesting Floquet
Hamiltonian with $q$-fold rotational lattice symmetry in phase space

$\displaystyle\hat{H}^{(T)}_{F}=\beta\left[(\hat{x}-i\hat{p})^{q}-1\right]\left[(\hat{x}+i\hat{p})^{q}-1\right].$
(9)

The discrete rotational symmetry can be described by
$\hat{R}(\frac{2\pi}{q})\hat{H}^{(T)}_{F}\hat{R}(\frac{2\pi}{q})=\hat{H}^{(T)}_{F}$,
where $\hat{R}(\theta)=\exp(i\hat{a}^{\dagger}\hat{a}\theta)$ is a phase-space
rotation by an angle $\theta$ Guo et al. (2013); Grimsmo et al. (2020). This
Hamiltonian supports $q$ global minima, cf. the $Q$-function in Fig. 1 (left)
for $q=6$. Here, we have rescaled the phase-space coordinates such that the
global minima fulfill $|x+ip|=1$. Each minimum corresponds to a different
classical solution. Remarkably, quantum fluctuations do not introduce any
tunneling between these solutions as the corresponding coherent states
$|\alpha\rangle=|e^{im\frac{2\pi}{q}}/\sqrt{2\lambda}\rangle$ with
$m=0,1,\cdots,q-1$ are exact zero-energy eigenstates. In other words, the
groundstate manifold is $q$-dimensional and is spanned by $q$ $q$-legged cat
states. The oscillator is steered into this manifold by photon decay (see the
Appendix E). The case $q=2$ has been proposed as a platform for fault tolerant
syndrome detection Puri et al. (2019).

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
$\hat{H}^{(T)}_{F\gamma}=U_{\gamma}\hat{H}^{(T)}_{F}U_{\gamma}^{\dagger}$ with
$U_{\gamma}=\exp[-\gamma\hat{a}^{\dagger}\hat{a}]$. Obviously,
$\lim_{\gamma\to 0}\hat{H}^{(T)}_{F\gamma}=\hat{H}^{(T)}_{F}$. We can
calculate analytically $f_{T}(k_{x},k_{p})$ and $V(x,t)$ for
$\hat{H}^{(T)}_{F\gamma}$ for any arbitrary positive integer $q$ and
$\gamma>0$. In the limit $\gamma\to 0$, the Fourier transform can be viewed as
a generalized function. This allows to perform the integral in Eq. (7) and
arrive at a closed expression for the driving potential (for more details see
the Appendix D)

$\displaystyle V(x,t)$ $\displaystyle=$
$\displaystyle\sum_{m=1}^{q}B_{q,m}\lambda^{q-m}x^{2m}-C_{q}\cos(q\omega_{0}t)x^{q},\
\ \ \ \ \ $ (10)

with $B_{q,m}=\frac{(2^{m}q!)^{2}(-1)^{q+m}}{(2m)!(q-m)!}$ and
$C_{q}=\frac{2\sqrt{\pi}q!}{\Gamma[(2q+3-(-1)^{q})/4]}$. We note that for
$q=2$ we recover a well-known result: Eq. (10) corresponds to a parametrically
driven Duffing oscillator Marthaler and Dykman (2006); Peano et al. (2012);
Dykman (2012); Bachtold et al. (2022). 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.

Figure 2:  Elliptical well in phase space. (a) Hamiltonian Q-function
$H^{(T)}_{Q}(x,p)/\beta$ with long axis $a=2$, short axis $b=1$, $\lambda=0.1$
and convolution factor $\sigma=\sqrt{\lambda}$. (b) Designed driving potential
$V(x,t)$ in one Floquet period (lower) and at instants $t=0$, $t=\pi/2$
(upper). (c)  Energy spectrum and Husimi Q-functions of ground state ($i=1$),
first excited state ($i=2$), fourth ($i=5$) excited eigenstates. The dashed
circles in (a) and (c) indicate the boundary of ellipitical well.

###  V.2 Sharp-boundary well

Next, we demonstrate that our method allows to engineer wells with a sharp
boundary in phase space. For concreteness we choose an elliptical shape, i.e.
$H^{(T)}_{Q}(x,p)=-\beta$ inside the white dashed line in Fig. 2(a) and
$H^{(T)}_{Q}=0$ otherwise. In the classical limit $\lambda\to 0$, our method
allows to find a closed form solution for $V(x,t)$ (see Appendix F). However,
our solution is divergent at two time-dependent positions. In addition, it
does not directly apply to the quantum regime, $\lambda\neq 0$, because the
dependence of $V(x,t)$ on $\lambda$ 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_{p}^{2}\to\infty$, for any $\lambda\neq 0$. We eliminate these
unphysical features by smoothing out the target Floquet Hamiltonian applying a
convolution with a gaussian kernel with standard deviation $\sigma$, cf. Fig.
2(a). For $\sigma$ above a threshold, $\sigma>\sqrt{\lambda/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 Appendix
F. This implies that we can implement a potential step that is arbitrarily
sharp when compared to the typical dimensions of the phase-space well, but
should remain smooth on the scale of the oscillator quantum fluctuations. We
point out that a sharp boundary narrower than quantum fluctuations
($\sigma<\sqrt{\lambda/2}$) does exist (Hamiltonian matrix elements in Fock
presentation are well defined, see Appendix B and Appendix F) but cannot be
realised by our present method. The spectrum and first few eigenmodes are also
shown in Fig. 2(c). The latter are squeezed non-gaussian states.

Figure 3:  Moiré superlattice in phase space. (a) Hamiltonian Q-function of
Moiré superlattice with twisted angle $\theta=10^{\circ}$ and confined in a
region with radius $R=20\pi$. (b) Density plot of NcFT coefficient $\beta
f_{T}(k_{x},k_{p})$. (c) Designed driving potential $V(x,t)$ for
$t\in[0,T/6)$. (d) $V(x,t)$ at fixed time instants $t=0$ (upper) and $t=T/12$
(lower).

###  V.3 Moiré superlattice

In Ref Guo et al. (2022), we have shown how to synthetize 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 Bistritzer and MacDonald
(2011); Cao et al. (2018a, b); Thomson et al. (2018). The Moiré superlattice
is formed by overlaying two honeycomb lattices with a relative twist angle
$\theta_{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
$\theta_{0}=10^{\circ}$ 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 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 (in our
dimensionless phase-space coordinates the lattice constant is $\pi$). A closed
formula for the Floquet Hamiltonian is given in the Appendix G.

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 phase-
space rotations and/or the rotation by the twist angle $\theta_{0}$, cf Fig.
3(a). The width of all the peaks is $\propto 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 Appendix G. In Fig. 3(c), we
plot the ensuing driving potential $V(x,t)$ for $0\leq 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\to\infty$, the peaks in $(k_{x},k_{p})$-space
become $\delta$-functions and the driving potential reduces to a discrete
sequence of stroboscopic lattices with specific amplitudes, wavelengths and
phases Guo et al. (2016); Liang et al. (2018); Guo et al. (2022); Guo (2021).
Considering that the contact interaction of cold atoms turns into a long-
distance Coulomb-like interaction in the rotating frame Sacha (2015);
Mierzejewski et al. (2017); Giergiel et al. (2018); Liang et al. (2018); Sacha
(2020); Guo and Liang (2020); Guo et al. (2022); Guo (2021); Hannaford and
Sacha (2022), many atoms in the phase space Moiré superlattice would mimic the
behaviour of electrons in twisted bilayer graphene Bistritzer and MacDonald
(2011); Cao et al. (2018a, b); Thomson et al. (2018).

###  V.4  Artificial atomic spectrum

As a final application, we show that our method allows to implement a target
spectrum $\\{E_{n}\\}$ as well as desired target eigenstates
$\\{|\psi_{n}\rangle\\}$. 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
$\hat{H}^{(T)}_{F}=\sum_{n}E_{n}|\psi_{n}\rangle\langle\psi_{n}|.$ For
concreteness, we consider $|\psi_{n}\rangle=|n\rangle$ where $\\{|n\rangle\\}$
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$,

$\displaystyle
H^{(T)}_{Q}(x,p)=e^{-\frac{x^{2}+p^{2}}{2\lambda}}\sum_{n=0}^{\infty}\frac{E_{n}}{n!}\Big{(}\frac{x^{2}+p^{2}}{2\lambda}\Big{)}^{n},$
(11)

and

$\displaystyle
f_{T}(k_{x},k_{p})=\sum_{n=0}^{\infty}\lambda\frac{E_{n}}{\beta}e^{\lambda\frac{k^{2}}{4}}{{}_{1}F_{1}}(1+n;1;-\lambda\frac{k^{2}}{2}),$
(12)

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\tau,k\sin\tau)$ is independent of the angular
coordinate $\tau$, 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
${{}_{1}F_{1}}(1+n;1;-\frac{k^{2}}{2})\sim(-\frac{\lambda
k^{2}}{2n!})^{n}e^{-\frac{\lambda k^{2}}{2}}$ for $k\to\infty$ 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
$E_{n}=-{\beta\lambda}/{(n+1)^{2}}$. In panel (b) we choose
$E_{1}=-\beta\lambda$ and $E_{0}=E_{1}-\beta\lambda(\lambda-\frac{3}{4})$
while all other levels are zero, $E_{n>1}=0$. Thus, at $\lambda=3/4$, the
energies $E_{0}$ and $E_{1}$ of the second spectrum display an exact crossing.

Figure 4:  Artificial spectrum: (a) designed hydrogen atomic levels with
parameter $\lambda=\frac{1}{4}$; (b) levels $|0\rangle$ and $|1\rangle$ gapped
from other degenerate levels with $\lambda=1$. In both figures, the
eigenstates $|n\rangle$ are the harmonic Fock states.

##  VI Experimental implementations

In order to design arbitrary Hamiltonians in phase space, one needs the
ability to engineer the driving real-space potential $V(x,t)$ in experiments.
As discussed above, cf. Eq. (7), the driving potential can be engineered by
superposing a series of cosine lattice potentials with tunable amplitudes,
wave vectors and phases. In cold atom experiments, this can be done with
optical lattices that are formed by laser beams intersecting at an angle
Moritz et al. (2003); Hadzibabic et al. (2004); Guo et al. (2022). In
experiments with sperconducting circuits Chen et al. (2014); Hofheinz et al.
(2011); Chen et al. (2011), a microwave cavity in series with a Josephson
junction (JJ) biased by a dc voltage ($V$) is described by the Hamiltonian
$\hat{\mathcal{H}}(t)=\hbar\omega_{0}\hat{a}^{\dagger}\hat{a}-E_{J}\cos[\omega_{J}t+\Delta(\hat{a}^{\dagger}+\hat{a})]$,
where $E_{J}$ is the JJ energy, $\omega_{J}=2eV/\hbar$ is the Josephson
frequency and $\Delta=\sqrt{2e^{2}/(\hbar\omega_{0}C)}$ with $C$ the cavity
capacitance Armour et al. (2013); Gramich et al. (2013); Leppäkangas et al.
(2013, 2015, 2016); Armour et al. (2015); Trif and Simon (2015); Kubala et al.
(2015b); Hofer et al. (2016); Dambach et al. (2017); Lang and Armour (2021);
Lang et al. (2022). In principle, more complicated driving potential can be
realized using more JJs and electric elements.

##  VII Summary and Outlook

In this work, we have introduced a powerful tool for Floquet engineering: A
general constructive method to derive the driving potential $V(x,t)$
generating any arbitrary target Floquet Hamiltonian $H^{(T)}_{F}(x,p)$ of a
single Bosonic mode. Here, we provide an approximation for $V(x,t)$ up to
leading order in the high-frequency Floquet-Magnus expansion. A natural
extension of our work would be to include higher-order perturbative
corrections, which can be viewed as the inverse problem of the Floquet-Magnus
theory. Another exciting prospect is to extend our method to a many-body
scenario. A building block for this extension is to upgrade the single-
particle plane-wave operator $\exp[{i(k_{x}\hat{x}+k_{p}\hat{p})}]$ used in
the NcFT Eq. (5) to a many-body equivalent,
$\exp[{\sum_{j}i(k^{j}_{x}\hat{x}_{j}+k^{j}_{p}\hat{p}_{j}})]$ where
($\hat{x}_{j},\hat{p}_{j}$) are the phase-space coordinates of the $j$-th
particle. In experiments with superconduncting circuits, this could be
implemented coupling a dc-voltage biased JJ to multiple superconducting
cavities Armour et al. (2013, 2015); Trif and Simon (2015); Hofer et al.
(2016); Dambach et al. (2017); Lang et al. (2022).

Acknowledgements

We acknowledge helpful discussions with Florian Marquardt and Muxin Han.

##  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 $\hat{H}^{(T)}_{F}=H^{(T)}_{F}(\hat{x},\hat{p})$. We start from
writing the target Hamiltonian as a sum of plane-wave operators, cf. Eq. (5)
in the main text,

$\displaystyle H^{(T)}_{F}(\hat{x},\hat{p})=\frac{\beta}{2\pi}\int\int
dk_{x}dk_{p}f_{T}(k_{x},k_{p})e^{i(k_{x}\hat{x}+k_{p}\hat{p})}.$ (13)

In order to calculate the Fourier coefficient $f_{T}(k_{x},k_{p})$, we first
express the target Hamiltonian with reordered ladder operators

$H^{(T)}_{F}(\hat{a}^{\dagger},\hat{a})\equiv\sum_{n,m}\chi_{nm}(\hat{a}^{\dagger})^{n}\hat{a}^{m}.$

Note that the ordering here keeps all the terms from commutators. By defining
the coherent state $|\alpha\rangle$ as the eigenstate of lowering operator
$\hat{a}|\alpha\rangle=\alpha|\alpha\rangle$, we calculate the operator in the
diagonal coherent representation

$\displaystyle
H^{(T)}_{Q}(\alpha,\alpha^{*})\equiv\langle\alpha|\hat{H}^{(T)}_{F}|\alpha\rangle=\sum_{n,m}\chi_{nm}(\alpha^{*})^{n}\alpha^{m}.$
(14)

Function $H^{(T)}_{Q}(\alpha,\alpha^{*})$ can also be written as
$H^{(T)}_{Q}(x,p)$ by identifying $\alpha=(x+ip)/\sqrt{2\lambda}$ where

$\displaystyle\left\\{\begin{array}[]{lll}x\equiv\langle\alpha|\hat{x}|\alpha\rangle=\sqrt{\frac{\lambda}{2}}(\alpha^{*}+\alpha)\\\
p\equiv\langle\alpha|\hat{p}|\alpha\rangle=i\sqrt{\frac{\lambda}{2}}(\alpha^{*}-\alpha).\\\
\end{array}\right.$ (17)

In order to calculation the NcFT coefficient $f_{T}(k_{x},k_{p})$ in Eq. (13),
we need to calculate the coherent diagonal element of the plane-wave operator
$\langle\alpha|e^{i(k_{x}\hat{x}+k_{p}\hat{p})}|\alpha\rangle$. For this
purpose, we introduce the displacement operator $\hat{D}_{\alpha}\equiv
e^{\alpha\hat{a}^{\dagger}-\alpha^{*}\hat{a}}$ with the following relationship
Guo and Liang (2020)

$\displaystyle
D_{\alpha}D_{\beta}=e^{i\mathrm{Im}(\alpha\beta^{*})}D_{\alpha+\beta},\ \
D_{\alpha}|\beta\rangle=e^{i\mathrm{Im}(\alpha\beta^{*})}|\alpha+\beta\rangle.$
(18)

We then write the plane-wave operator as
$e^{i(k_{x}\hat{x}+k_{p}\hat{p})}=\hat{D}_{-\sqrt{\frac{\lambda}{2}}(k_{p}-ik_{x})}$.
Using the relationship (18), we have the matrix element of plane-wave operator
$e^{i(k_{x}\hat{x}+k_{p}\hat{p})}$ in coherent state representation

$\displaystyle\langle\alpha|e^{i(k_{x}\hat{x}+k_{p}\hat{p})}|\beta\rangle$
$\displaystyle=$
$\displaystyle\langle\alpha|\hat{D}_{-\sqrt{\frac{\lambda}{2}}(k_{p}-ik_{x})}|\beta\rangle$
(19) $\displaystyle=$
$\displaystyle\Big{\langle}\alpha\Big{|}\beta-\sqrt{\frac{\lambda}{2}}(k_{p}-ik_{x})\Big{\rangle}e^{-i\sqrt{\frac{\lambda}{2}}\mathrm{Im}[(k_{p}-ik_{x})\beta^{*}]}$
$\displaystyle=$ $\displaystyle
e^{-\frac{1}{2}|\alpha|^{2}-\frac{1}{2}|\beta-\sqrt{\frac{\lambda}{2}}(k_{p}-ik_{x})|^{2}+\alpha^{*}\beta-\alpha^{*}\sqrt{\frac{\lambda}{2}}(k_{p}-ik_{x})-i\sqrt{\frac{\lambda}{2}}\mathrm{Im}[(k_{p}-ik_{x})\beta^{*}]}$

In the last step, we have used and the identity
$\langle\alpha|\beta\rangle=e^{-|\alpha|^{2}/2-|\beta|^{2}/2+\alpha^{*}\beta}$.
Thus, we have the diagonal elements of plane-wave operator from Eq. (19)

$\displaystyle\langle\alpha|e^{i(k_{x}\hat{x}+k_{p}\hat{p})}|\alpha\rangle=\exp\big{(}-\frac{\lambda}{4}|k_{p}-ik_{x}|^{2}\big{)}e^{i(k_{x}x+k_{p}p)}=e^{-\frac{\lambda}{4}(k_{x}^{2}+k_{p}^{2})}e^{i(k_{x}x+k_{p}p)}.$
(20)

Using Eqs. (14) and (20), we have the Fourier coefficient from Eq. (13) as
follows

$\displaystyle
f_{T}(k_{x},k_{p})=\frac{e^{\frac{\lambda}{4}(k^{2}_{x}+k^{2}_{p})}}{2\pi\beta}\int\int
dxdpH^{(T)}_{Q}(x,p)e^{-i(k_{x}x+k_{p}p)}.\ $ (21)

Eqs. (13) and (21) construct the noncommutative Fourier transformation (NcFT)
technique introduced in this paper. From the hermicity of Hamiltonian
operator, we have the following important relationship

$\displaystyle f_{T}(k_{x},k_{p})=f_{T}^{*}(-k_{x},-k_{p}).$ (22)

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. (14). One may wonder if some information is lost by
neglecting the off-diagonal elements. In this section, we will prove the
Hamiltonian operator $H^{(T)}_{F}(\hat{x},\hat{p})$, given by Eqs. (13) and
(21), is fully determined by its Hamiltonian Q-function
$H^{(T)}_{Q}(x,p)\equiv\langle\alpha|H^{(T)}_{F}(\hat{x},\hat{p})|\alpha\rangle$
together with commutator $[\hat{x},\hat{p}]=i\lambda$.

We write the Hamiltonian in the Fock representation
$\hat{H}^{(T)}_{F}=\sum_{n,m}\xi_{nm}|n\rangle\langle m|$ with
$n,m=0,1,\cdots$ and define the following operator

$\displaystyle H_{nm}(\hat{x},\hat{p})$ $\displaystyle\equiv$
$\displaystyle\frac{1}{2\pi}\int\int
dk_{x}dk_{p}f_{nm}(k_{x},k_{p})e^{i(k_{x}\hat{x}+k_{p}\hat{p})}.$ (23)

where $f_{nm}(k_{x},k_{p})$ is the NcFT coefficient of the operator
$|n\rangle\langle m|$ given by

$\displaystyle f_{nm}(k_{x},k_{p})$ $\displaystyle=$
$\displaystyle\frac{e^{\frac{\lambda}{4}(k^{2}_{x}+k^{2}_{p})}}{2\pi}\int\int
dxdpH_{nm}(x,p)e^{-i(k_{x}x+k_{p}p)}$ (24)

with the Q-function of operator $|n\rangle\langle m|$ given by

$\displaystyle H_{nm}(x,p)=\langle\alpha|n\rangle\langle m|\alpha\rangle.$
(25)

Because the target Hamiltonian is the liner superposition of $|n\rangle\langle
m|$ with $n,m=0,1,\cdots$, we just need to prove

$\displaystyle\langle
n^{\prime}|H_{nm}(\hat{x},\hat{p})|m^{\prime}\rangle=\delta_{n,n^{\prime}}\delta_{m,m^{\prime}}.$
(26)

Using coherent state in the basis of Fock states
$|\alpha\rangle=e^{-\frac{|\alpha|^{2}}{2}}\sum_{n}\frac{\alpha^{n}}{\sqrt{n!}}|n\rangle$,
we calculate the Q-function of $|n\rangle\langle m|$

$\displaystyle H_{nm}(x,p)$ $\displaystyle=$
$\displaystyle\langle\alpha|n\rangle\langle
m|\alpha\rangle=\frac{1}{\sqrt{n!m!}}e^{-\frac{x^{2}+p^{2}}{2\lambda}}\Big{(}\frac{x-ip}{\sqrt{2\lambda}}\Big{)}^{n}\Big{(}\frac{x+ip}{\sqrt{2\lambda}}\Big{)}^{m}.$
(27)

By introducing $x=r\cos\theta,\ p=r\sin\theta$ and $x=k\cos\tau,\
p=k\sin\tau$, we have the Fourier component

$\displaystyle f_{nm}(k,\tau)$ $\displaystyle=$
$\displaystyle\frac{e^{\frac{\lambda}{4}k^{2}}}{2\pi}\int\int rdrd\theta
H_{nm}(r\cos\theta,r\sin\theta)e^{-ikr\cos(\theta-\tau)}$ (28)
$\displaystyle=$
$\displaystyle\frac{e^{i(m-n)\tau}}{\sqrt{n!m!}}\Big{(}\frac{1}{\sqrt{2\lambda}}\Big{)}^{m+n}\frac{e^{\frac{\lambda}{4}k^{2}}}{2\pi}\int_{0}^{\infty}r^{m+n+1}e^{-\frac{r^{2}}{2\lambda}}dr\int_{0}^{2\pi}d\theta
e^{-ikr\cos(\theta-\tau)+i(m-n)(\theta-\tau)}$ $\displaystyle=$
$\displaystyle\frac{e^{i(m-n)\tau}}{\sqrt{n!m!}}\Big{(}\frac{1}{\sqrt{2\lambda}}\Big{)}^{m+n}e^{\frac{\lambda}{4}k^{2}}i^{n-m}\int_{0}^{\infty}r^{m+n+1}e^{-\frac{r^{2}}{2\lambda}}J_{n-m}(-kr)dr$
$\displaystyle=$ $\displaystyle
e^{\frac{\lambda}{4}k^{2}}\sqrt{\frac{n!}{m!}}\Big{(}ie^{i\tau}\frac{1}{k}\sqrt{\frac{2}{\lambda}}\Big{)}^{m-n}\frac{\lambda}{\Gamma(1-m+n)}{{}_{1}F_{1}}(1+n;1-m+n;-\frac{\lambda}{2}k^{2})$

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 Guo and Marthaler (2016)

$\displaystyle\langle
n^{\prime}|e^{i(k_{x}\hat{x}+k_{p}\hat{p})}|m^{\prime}\rangle=e^{-\frac{\lambda}{4}(k^{2}_{x}+k^{2}_{p})}\Big{[}\sqrt{\frac{\lambda}{2}}(k_{p}+ik_{x})\Big{]}^{m^{\prime}-n^{\prime}}\sqrt{\frac{n^{\prime}!}{m^{\prime}!}}L_{n^{\prime}}^{m^{\prime}-n^{\prime}}\big{[}\frac{\lambda}{2}(k^{2}_{x}+k^{2}_{p})\big{]}$
(29)

where $L_{n^{\prime}}^{m^{\prime}-n^{\prime}}(z)$ is the generalized Laguerre
polynomial. Then, we have the matrix element of $H_{nm}(\hat{x},\hat{p})$ in
Fock representation

$\displaystyle\langle n^{\prime}|H_{nm}(\hat{x},\hat{p})|m^{\prime}\rangle$
$\displaystyle=$ $\displaystyle\frac{1}{2\pi}\int\int
dk_{x}dk_{p}f_{nm}(k_{x},k_{p})\langle
n^{\prime}|e^{i(k_{x}\hat{x}+k_{p}\hat{p})}|m^{\prime}\rangle$ (30)
$\displaystyle=$ $\displaystyle
e^{i\frac{\pi}{2}(m^{\prime}-n^{\prime})}\sqrt{\frac{n^{\prime}!}{m^{\prime}!}}\sqrt{\frac{n!}{m!}}\Big{(}i\sqrt{\frac{2}{\lambda}}\Big{)}^{m-n}\frac{\lambda}{\Gamma(1-m+n)}$
$\displaystyle\times\frac{1}{2\pi}\int_{0}^{2\pi}e^{i[(m-n)-(m^{\prime}-n^{\prime})]\tau}d\tau$
$\displaystyle\times\int_{0}^{\infty}dk\Big{(}\sqrt{\frac{\lambda}{2}}k\Big{)}^{m^{\prime}-n^{\prime}}k^{n-m+1}{{}_{1}F_{1}}(1+n;1-m+n;-\frac{\lambda}{2}k^{2})L_{n^{\prime}}^{m^{\prime}-n^{\prime}}\big{(}\frac{\lambda}{2}k^{2}\big{)}$
$\displaystyle=$
$\displaystyle\delta_{m-n,m^{\prime}-n^{\prime}}(-1)^{m-n}\sqrt{\frac{n^{\prime}!}{(m-n+n^{\prime})!}}\sqrt{\frac{n!}{m!}}\frac{\lambda}{\Gamma(1-m+n)}$
$\displaystyle\times\int_{0}^{\infty}kdk\
{{}_{1}F_{1}}(1+n;1-m+n;-\frac{\lambda}{2}k^{2})L_{n^{\prime}}^{m-n}\big{(}\frac{\lambda}{2}k^{2}\big{)}$
$\displaystyle=$
$\displaystyle\delta_{m-n,m^{\prime}-n^{\prime}}(-1)^{m-n}\sqrt{\frac{n^{\prime}!}{(m-n+n^{\prime})!}}\sqrt{\frac{n!}{m!}}\frac{2}{\Gamma(1-m+n)}$
$\displaystyle\times\int_{0}^{\infty}\tilde{k}d\tilde{k}\
{{}_{1}F_{1}}(1+n;1-m+n;-\tilde{k}^{2})L_{n^{\prime}}^{m-n}\big{(}\tilde{k}^{2}\big{)}\
\ \ (\mathrm{where}\ \ \tilde{k}=\sqrt{\frac{\lambda}{2}}k)$ $\displaystyle=$
$\displaystyle\delta_{n,n^{\prime}}\delta_{m,m^{\prime}}.$

This is the identity (26) we aim to prove. As a result, for a fixed $\lambda$,
the mapping between Floquet Hamiltonians and $Q$-functions is one to one.

Note that the exponentially suppression factor
$e^{-\frac{\lambda}{4}(k^{2}_{x}+k^{2}_{p})}$ in Eq. (29) cancels the same
exponentially increasing factor in Eq. (28). We will mention the sequence of
this point in the discussion of sharp-boundary elliptical well below.

##  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\tau,k_{p}=k\sin\tau$), and write the
Fourier expansion Eq. (13) as

$\displaystyle\hat{H}^{(T)}_{F}$ $\displaystyle=$
$\displaystyle\frac{1}{2\pi}\int_{0}^{2\pi}d\tau\int_{-\infty}^{+\infty}dk\frac{|k|}{2}f_{T}(k,\tau)e^{ik(\hat{x}\cos\tau+\hat{p}\sin\tau)}.\
\ \ \ $ (31)

Here, we have defined the Fourier component in the polar coordinate system via
$f_{T}(k,\tau)\equiv f_{T}(k_{x},k_{p})$, and allow for negative $k$ via the
relation $f_{T}(k,\tau)\equiv f^{*}_{T}(-k,\tau)$, cf. Eq.(22). From the
Fourier component of the target Hamiltonian $f_{T}(k,\tau)$, we design the
real space driving potential as follows

$\displaystyle V(x,t)$ $\displaystyle=$
$\displaystyle\int_{-\infty}^{+\infty}\frac{|k|}{2}f_{T}(k,\omega_{0}t)e^{ikx}dk,$
(32)

where the phase variable $\omega_{0}t$ in time domain plays the role of angle
$\tau$ in $(k_{x},k_{p})$ space. By setting driving period
$T_{d}=2\pi/(\omega_{0}q)$, the Hamiltonian in the rotating frame, cf. Eq. (3)
in the main text, becomes

$\displaystyle H(t)$ $\displaystyle=$
$\displaystyle\int_{-\infty}^{+\infty}\frac{1}{2}|k|f_{T}(k,\omega_{0}t)e^{ik(\hat{x}\cos\omega_{0}t+\hat{p}\sin\omega_{0}t)}dk.\
\ \ $ (33)

By time averaging the above Hamiltonian, cf. Eq. (4) in the main text, and
comparing the averaged result to Eq. (31), one can directly find that the
lowest-order Floquet-Magnus expansion (RWA) gives the target Hamiltonian
$H^{(T)}_{F}(\hat{x},\hat{p})$. We can also write the engineered driving
potential in real space as

$V(x,t)=\int_{0}^{+\infty}|kf(k,\omega_{0}t)|\cos[kx+\phi(k,t)]dk,$

where we have introduced phase $\phi(k,t)=\mathrm{Arg}[f_{T}(k,\omega_{0}t)]$
and used the property $f(-k,\omega_{0}t)=f^{*}(k,\omega_{0}t)$. Thus, the
driving potential can be engineered by superposing a series of cosine lattice
potentials with tunable amplitudes $|kf_{T}(k,\omega_{0}t)|$ and phases
$\phi(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

$\displaystyle\hat{H}^{(T)}_{F\gamma}=\frac{\beta}{|\alpha_{0}|^{2q}}e^{-\gamma\hat{a}^{\dagger}\hat{a}}(\hat{a}^{\dagger
q}-\alpha_{0}^{*q})(\hat{a}^{q}-\alpha_{0}^{q})e^{-\gamma\hat{a}^{\dagger}\hat{a}},$
(34)

where the factor $e^{-\gamma\hat{a}^{\dagger}\hat{a}}$ with $\gamma>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
$\alpha_{0}=1/\sqrt{2\lambda}$. Using the identity

$e^{-\gamma\hat{a}^{\dagger}\hat{a}}|\alpha\rangle=e^{-\frac{1}{2}(1-e^{-2\gamma})|\alpha|^{2}}|\alpha
e^{\gamma}\rangle,$

we have the Hamiltonian Q-function as follows

$\displaystyle H^{(T)}_{Q\gamma}(x,p)$ $\displaystyle=$
$\displaystyle\langle\alpha|\hat{H}^{(T)}_{F\gamma}|\alpha\rangle$ (35)
$\displaystyle=$
$\displaystyle\frac{\beta}{|\alpha_{0}e^{\gamma}|^{2q}}\exp(-\frac{x^{2}+p^{2}}{2\lambda\sigma^{2}_{\gamma}})\Big{|}\Big{(}\frac{x+ip}{\sqrt{2\lambda}}\Big{)}^{q}-\alpha^{q}_{0}e^{q\gamma}\Big{|}^{2}.\
\ $

Here, we have defined the parameter $\sigma_{\gamma}=1/\sqrt{1-e^{-2\gamma}}$.

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\phi,p=r\sin\phi)$ and
$(k_{x}=k\cos\tau,k_{p}=k\sin\tau)$. Plugging Eq. (35) into Eq. (21), we have

$\displaystyle f_{T}(k,\tau)$ $\displaystyle=$ $\displaystyle\frac{e^{\lambda
k^{2}/4}}{2\pi\beta}\int_{0}^{+\infty}rdr\int_{0}^{2\pi}d\phi
H^{(T)}_{Q\gamma}(r\cos\phi,r\sin\phi)e^{-ikr\cos(\phi-\tau)}$ (36)
$\displaystyle=$ $\displaystyle\frac{e^{\lambda
k^{2}/4}}{|\alpha_{0}e^{\gamma}|^{2q}}\int_{0}^{+\infty}rdre^{-\frac{r^{2}}{2\lambda\sigma^{2}_{\gamma}}}\Big{[}\frac{r^{2q}}{(2\lambda)^{q}}J_{0}(kr)-\frac{(ir\alpha_{0}e^{\gamma-i\tau})^{q}}{(2\lambda)^{\frac{q}{2}}}J_{-q}(kr)-\frac{(-ir\alpha^{*}_{0}e^{\gamma+i\tau})^{q}}{(2\lambda)^{\frac{q}{2}}}J_{q}(kr)+|\alpha_{0}e^{\gamma}|^{2q}J_{0}(kr)\Big{]}$
$\displaystyle=$ $\displaystyle\frac{\lambda
e^{\frac{\lambda}{4}k^{2}}\sigma^{2(q+1)}_{\gamma}}{|\alpha_{0}e^{\gamma}|^{2q}}\Big{[}n!{{}_{1}F_{1}}(1+q;1;-\frac{\lambda}{2}\sigma^{2}_{\gamma}k^{2})-\Big{(}-ie^{-i\tau}\alpha_{0}e^{\gamma}\sqrt{\frac{\lambda}{2}}\Big{)}^{q}k^{q}{{}_{1}F_{1}}(1+q;1+q;-\frac{\lambda}{2}\sigma^{2}_{\gamma}k^{2})$
$\displaystyle-\Big{(}-ie^{i\tau}\alpha_{0}^{*}e^{\gamma}\sqrt{\frac{\lambda}{2}}\Big{)}^{q}k^{q}{{}_{1}F_{1}}(1+q;1+q;-\frac{\lambda}{2}\sigma^{2}_{\gamma}k^{2})+\Big{|}\frac{\alpha_{0}e^{\gamma}}{\sigma_{\gamma}}\Big{|}^{2q}{{}_{1}F_{1}}(1;1;-\frac{\lambda}{2}\sigma^{2}_{\gamma}k^{2})\Big{]}.$

Here, $J_{q}(\bullet)$ is the Bessel function of $q$-th order and
${{}_{1}F_{1}}(a;b;\bullet)$ 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:

$\displaystyle\mathrm{Identity\ I:}\ $
$\displaystyle\int_{-\infty}^{+\infty}\frac{1}{2}|k|e^{\frac{\lambda}{4}k^{2}}k^{q}{{}_{1}F_{1}}(1+q;1+q;-\frac{\lambda}{2}\sigma^{2}_{\gamma}k^{2})e^{ikx}dk$
(37)
$\displaystyle=2^{q-1}(2\lambda\sigma^{2}_{\gamma}-\lambda)^{-\frac{q+2}{2}}\Big{[}q\Gamma(\frac{q}{2})(1+(-1)^{q}){{}_{1}F_{1}}(\frac{q+2}{2};\frac{1}{2};-\frac{x^{2}}{2\lambda\sigma^{2}_{\gamma}-\lambda})$
$\displaystyle+4ix(1-(-1)^{n})\Gamma(\frac{q+3}{2}){{}_{1}F_{1}}(\frac{q+3}{2};\frac{3}{2};-\frac{x^{2}}{2\lambda\sigma^{2}_{\gamma}-\lambda})\Big{]};$
$\displaystyle\mathrm{Identity\ II:}\ $
$\displaystyle\int_{-\infty}^{+\infty}\frac{1}{2}|k|e^{\frac{\lambda}{4}k^{2}}{{}_{1}F_{1}}(1;1;-\frac{\lambda}{2}\sigma^{2}_{\gamma}k^{2})e^{ikx}dk=\frac{2}{2\lambda\sigma^{2}_{\gamma}-\lambda}-\frac{4x}{(2\lambda\sigma^{2}_{\gamma}-\lambda)^{\frac{3}{2}}}D\big{(}\frac{x}{\sqrt{2\lambda\sigma^{2}_{\gamma}-\lambda}}\big{)};$
(38) $\displaystyle\mathrm{Identity\ III:}\ $
$\displaystyle\int_{-\infty}^{+\infty}\frac{1}{2}|k|e^{\frac{\lambda}{4}k^{2}}{{}_{1}F_{1}}(q+1;1;-\frac{\lambda}{2}\sigma^{2}_{\gamma}k^{2})e^{ikx}dk$
(39)
$\displaystyle=\int_{0}^{+\infty}ke^{\frac{\lambda}{4}k^{2}}{{}_{1}F_{1}}(q+1;1;-\frac{\lambda}{2}\sigma^{2}_{\gamma}k^{2})\cos(kx)dk$
$\displaystyle=\frac{1}{2}\sum_{m=0}^{+\infty}\frac{(-1)^{m}}{(2m)!}\int_{0}^{+\infty}e^{\frac{\lambda}{4}k^{2}}{{}_{1}F_{1}}(q+1;1;-\frac{\lambda}{2}\sigma^{2}_{\gamma}k^{2})(kx)^{2m}dk^{2}$
$\displaystyle=\frac{1}{2}\sum_{m=0}^{+\infty}\frac{(-1)^{m}}{(2m)!}x^{2m}\int_{0}^{+\infty}e^{\frac{\lambda}{4}z}{{}_{1}F_{1}}(q+1;1;-\frac{\lambda}{2}\sigma^{2}_{\gamma}z)z^{m}dz\
(\mathrm{here},\ z=k^{2})$
$\displaystyle=\frac{1}{2}\sum_{m=0}^{+\infty}\frac{(-1)^{m}}{(2m)!}x^{2m}\big{(}-\frac{4}{\lambda}\big{)}^{m+1}\Gamma(m+1){{}_{2}F_{1}}(m+1,q+1;1;2\sigma^{2}_{\gamma}).$

Here, $\Gamma(\bullet)$ is the Gamma function and
$D(z)=e^{-z^{2}}\int_{0}^{z}e^{t^{2}}dt$ is the Dawson function. Identity Eq.
(38) is the special case of identity Eq. (37) by setting $n=0$. In identity
Eq. (39), ${{}_{2}F_{1}}(a,b;c;z)$ is the hypergeometric function given by the
path integral in the complex $\zeta$-plane Srivastava and Choi (2012)

$\displaystyle{{}_{2}F_{1}}(a,b;c;z)=\frac{\Gamma(c)}{\Gamma(a)\Gamma(b)}\frac{1}{2\pi
i}\int_{-i\infty}^{+i\infty}\frac{\Gamma(a+\zeta)\Gamma(b+\zeta)\Gamma(-\zeta)}{\Gamma(c+\zeta)}(-z)^{\zeta}d\zeta.$
(40)

The above integral is valid for $|\arg(-z)|\leq\pi-\epsilon$
($0<\epsilon<\pi$) and $a,b\notin\mathbb{Z}^{-}_{0}$. For $|z|<1$, the
hypergeometric function can be written by the power series

$\displaystyle{{}_{2}F_{1}}(a,b;c;z)=\sum_{k=0}^{\infty}\frac{(a)_{k}(b)_{k}}{(c)_{k}}\frac{z^{k}}{k!},\
\ \ |z|<1$ (41)

with the Pochhammer symbol $(x)_{k}=\Gamma(x+k)/\Gamma(x)$. The analytical
continuation of ${{}_{2}F_{1}}(a,b;c;z)$ into the domain $|z|>1$ can be
realised from the following relationship Srivastava and Choi (2012)

$\displaystyle{{}_{2}F_{1}}(a,b;c;z)$ $\displaystyle=$
$\displaystyle\frac{\Gamma(c)\Gamma(b-a)}{\Gamma(b)\Gamma(c-a)}\big{(}-\frac{1}{z}\big{)}^{a}{{}_{2}F_{1}}(a,1-c+a;1-b+a;\frac{1}{z})$
(42)
$\displaystyle+\frac{\Gamma(c)\Gamma(a-b)}{\Gamma(a)\Gamma(c-b)}\big{(}-\frac{1}{z}\big{)}^{b}{{}_{2}F_{1}}(b,1-c+b;1-a+b;\frac{1}{z}).$

The above relationship is valid for $|\arg(-z)|\leq\pi-\epsilon$
($0<\epsilon<\pi$) and $a-b\notin\mathbb{Z}$. From Eqs. (32) and (36)-(39), we
obtain the analytical expression for the designed driving potential for finite
values of $\gamma>0$ and arbitrary complex number $\alpha_{0}$ as follows

$\displaystyle V_{\gamma}(x,t)$ $\displaystyle=$
$\displaystyle-\frac{(\lambda\sigma^{2}_{\gamma})^{q+1}\Gamma(q+1)}{|\sqrt{\lambda}\alpha_{0}e^{\gamma}|^{2q}}\sum_{m=0}^{+\infty}\frac{\Gamma(m+1)}{\Gamma(2m+1)}\frac{x^{2m}}{2}\big{(}\frac{4}{\lambda}\big{)}^{m+1}{{}_{2}F_{1}}(m+1,q+1;1;2\sigma^{2}_{\gamma})$
(43)
$\displaystyle-\frac{(\lambda\sigma^{2}_{\gamma})^{q+1}2^{q-1}}{|\sqrt{\lambda}\alpha_{0}e^{\gamma}|^{2q}(2\lambda\sigma^{2}_{\gamma}-\lambda)^{\frac{q+3}{2}}}\Big{[}\Big{(}-ie^{\gamma-i\tau}\alpha_{0}\sqrt{\frac{\lambda}{2}}\Big{)}^{q}+\Big{(}-ie^{\gamma+i\tau}\alpha^{*}_{0}\sqrt{\frac{\lambda}{2}}\Big{)}^{q}\Big{]}$
$\displaystyle\times\Big{[}(1+(-1)^{q})q\Gamma(\frac{q}{2})\sqrt{2\lambda\sigma^{2}_{\gamma}-\lambda}{{}_{1}F_{1}}(\frac{q+2}{2};\frac{1}{2};-\frac{x^{2}}{2\lambda\sigma^{2}_{\gamma}-\lambda})$
$\displaystyle+4ix(1-(-1)^{q})\Gamma(\frac{q+3}{2}){{}_{1}F_{1}}(\frac{q+3}{2};\frac{3}{2};-\frac{x^{2}}{2\lambda\sigma^{2}_{\gamma}-\lambda})\Big{]}$
$\displaystyle+\frac{2\sigma^{2}_{\gamma}}{2\sigma^{2}_{\gamma}-1}-\frac{1}{\sqrt{\lambda}}\frac{4x\sigma^{2}_{\gamma}}{(2\sigma^{2}_{\gamma}-1)^{\frac{3}{2}}}D\big{(}\frac{x}{\sqrt{2\lambda\sigma^{2}_{\gamma}-\lambda}}\big{)}.$

Next, we discuss how to calculate driving potential $V_{\gamma}(x,t)$ in the
limit of $\gamma\rightarrow 0$
($\sigma_{\gamma}=1/\sqrt{1-e^{-2\gamma}}\rightarrow+\infty$). Using Eq. (41),
(42) and the Euler’s reflection formula $\Gamma(z)\Gamma(1-z)=\pi/\sin(\pi
z)$, we obtain the following series expansion of the hypergeometric function
for $|z|>1$ and $q\in\mathbb{Z}^{+}_{0}$

$\displaystyle{{}_{2}F_{1}}(1+m,1+q;1;z)=\sum_{k=0}^{+\infty}\frac{\Gamma(m+k+1)^{2}(-1)^{m+q+1}}{\Gamma(1+m)\Gamma(1+q)\Gamma(m-q+k+1)}\frac{1}{k!}\big{(}\frac{1}{z}\big{)}^{k+m+1}.$
(44)

Note that the parameter $m$ actually can take the whole real values, i.e.,
$m\in\mathbb{R}$. Although the above expansion is not defined for $m$ at
negative integers $m\in\mathbb{Z}^{-}$, but the limit values do exist and can
be defined as the values of the expansion for $m\in\mathbb{Z}^{-}$ . By
plugging the above series expansion Eq. (44) and also the confluent
hypergeometric function
${{}_{1}F_{1}}(a;b;z)=\sum_{k=0}^{\infty}\frac{(a)_{k}}{(b)_{k}}\frac{z^{k}}{k!}$
into Eq. (43), we obtain the driving potential in the limit of
$\gamma\rightarrow 0$ ($\sigma_{\gamma}\rightarrow+\infty$)

$\displaystyle V_{\gamma\rightarrow 0}(x,t)$ $\displaystyle=$
$\displaystyle\frac{1}{|\sqrt{\lambda}\alpha_{0}|^{2q}}\sum_{m=0}^{+\infty}\sum_{k=0}^{+\infty}\frac{2^{2m-q}\lambda^{q-m}}{\Gamma(2m+1)}\frac{\Gamma(m+k+1)^{2}(-1)^{m+q}}{\Gamma(m-q+k+1)}\frac{x^{2m}}{k!}\big{(}\frac{1}{2\sigma^{2}_{\gamma}}\big{)}^{k+m-q}$
(49)
$\displaystyle-\frac{1}{|\sqrt{\lambda}\alpha_{0}|^{2q}}\Big{[}\Big{(}-ie^{\gamma-i\tau}\alpha_{0}\sqrt{\frac{\lambda}{2}}\Big{)}^{q}+\Big{(}-ie^{\gamma+i\tau}\alpha^{*}_{0}\sqrt{\frac{\lambda}{2}}\Big{)}^{q}\Big{]}$
$\displaystyle\times\Big{[}(1+(-1)^{q})\sum_{m=0}^{+\infty}\frac{\Gamma(1+m+\frac{q}{2})\Gamma(\frac{1}{2})}{\Gamma(\frac{1}{2}+m)}(-1)^{m}2^{\frac{q}{2}-m-1}(\lambda\sigma_{\gamma}^{2})^{\frac{q}{2}-m}x^{2m}$
$\displaystyle+i(1-(-1)^{q})\sum_{m=0}^{+\infty}\frac{\Gamma(\frac{3}{2}+m+\frac{q}{2})\Gamma(\frac{3}{2})}{\Gamma(\frac{3}{2}+m)}(-1)^{m}2^{\frac{q}{2}-m-\frac{1}{2}}(\lambda\sigma_{\gamma}^{2})^{\frac{q}{2}-m-\frac{1}{2}}x^{2m+1}\Big{]}$
$\displaystyle+1.$

In line (49), only terms that satisfy $k+m-q=0$ give nonzero contribution
otherwise
$\frac{1}{\Gamma(m+k-q+1)}\big{(}\frac{1}{2\sigma^{2}_{\gamma}}\big{)}^{k+m-q}=0$
in the limit of $\sigma_{\gamma}=+\infty$ (note that Gamma function
$\Gamma(z)=\infty$ for nonpositive integer argument $z\in\mathbb{Z}_{0}^{-}$).
In line (49), only terms with even integer $q$ and $m\leq 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^{\pm iq\tau}$ in line (49). Therefore, the only nontrivial contribution
comes from the term with $m=q/2$ in line (49). For the same reason, only
nontrivial contribution comes from the term with $m=(q-1)/2$ in line (49). By
neglecting these terms, we have the designed driving potential

$\displaystyle V(x,t)$ $\displaystyle=$
$\displaystyle\frac{1}{|\sqrt{2\lambda}\alpha_{0}|^{2q}}\sum_{m=0}^{q}(-1)^{m+q}\frac{(2^{m}q!)^{2}}{(2m)!(q-m)!}\lambda^{q-m}x^{2m}$
(50)
$\displaystyle-\frac{(-i)^{q}}{|\sqrt{2\lambda}\alpha_{0}|^{2q}}\Big{[}\Big{(}e^{-i\tau}\alpha_{0}\sqrt{2\lambda}\Big{)}^{q}+\Big{(}e^{i\tau}\alpha^{*}_{0}\sqrt{2\lambda}\Big{)}^{q}\Big{]}$
$\displaystyle\times\Big{(}[1+(-1)^{q}]\frac{\Gamma(1+q)\Gamma(\frac{1}{2})}{\Gamma(\frac{1}{2}+\frac{q}{2})}\frac{(-1)^{\frac{q}{2}}}{2}+i[1-(-1)^{q}]\frac{\Gamma(q+1)\Gamma(\frac{3}{2})}{\Gamma(\frac{q}{2}+1)}(-1)^{\frac{q-1}{2}}\Big{)}x^{q}$
$\displaystyle+1.$

By taking the value of $\alpha_{0}=1/\sqrt{2\lambda}$, we obtain the driving
potential given by Eq. (10) shown in the main text, i.e.,

$\displaystyle V(x,t)$ $\displaystyle=$
$\displaystyle\sum_{m=0}^{q}\frac{(2^{m}q!)^{2}(-1)^{q+m}}{(2m)!(q-m)!}\lambda^{q-m}x^{2m}-\sqrt{\pi}q!\Bigg{[}\frac{1+(-1)^{q}}{\Gamma(\frac{q}{2}+\frac{1}{2})}+\frac{1-(-1)^{q}}{\Gamma(\frac{q}{2}+1)}\Bigg{]}\cos(q\tau)x^{q}+1$
(51) $\displaystyle=$
$\displaystyle\sum_{m=0}^{q}B_{q,m}\lambda^{q-m}x^{2m}-C_{q}\cos(q\omega_{0}t)x^{q}+1.\
\ \ \ \ \ $

Here, we have defined the coefficients
$B_{q,m}=\frac{(2^{m}q!)^{2}(-1)^{q+m}}{(2m)!(q-m)!}$ and
$C_{q}=\sqrt{\pi}q!\Big{[}\frac{1+(-1)^{q}}{\Gamma(\frac{q}{2}+\frac{1}{2})}+\frac{1-(-1)^{q}}{\Gamma(\frac{q}{2}+1)}\Big{]}=\frac{2\sqrt{\pi}q!}{\Gamma[(2q+3-(-1)^{q})/4]}$.

##  Appendix E Ground state property of rotational lattice Hamiltonian

In the main text, we have claimed that the groundstate manifold of present
rotational lattice Hamiltonian is $q$-dimensional and is spanned by $q$
$q$-legged cat states. This is because the $q$ coherent states
$|\alpha_{m}\rangle=|e^{im\frac{2\pi}{q}}/\sqrt{2\lambda}\rangle$ with
$m=0,1,\cdots,q-1$ are exact zero-energy eigenstates and quantum fluctuations
do not introduce any tunneling between these coherent states. Mathematically,
the $q$-fold symmetry of Hamiltonian in phase space is described by

$\displaystyle\hat{R}_{q}\hat{H}^{(T)}_{F}\hat{R}_{q}^{\dagger}=\hat{H}^{(T)}_{F},$
(52)

where $\hat{R}_{q}\equiv e^{i\hat{a}^{\dagger}\hat{a}\frac{2\pi}{q}}$ is the
discrete rotational operation Guo et al. (2013); Grimsmo et al. (2020).
According to the Bloch theorem extended in phase space Guo et al. (2013), the
eigentsates of Hamiltonian can be written in form of

$\displaystyle|\psi_{l,m}\rangle=\frac{1}{\sqrt{q}}\sum_{p=0}^{q-1}e^{imp\frac{2\pi}{q}}\hat{R}_{q}^{p}|\phi_{l}\rangle,$
(53)

where index $l$ labels the Bloch bands, $m$ is called quasinumber and
$|\phi_{l}\rangle$ is the Wannier state of $l$-th Bloch band. The Wannier
state $|\phi_{l}\rangle$ for the lowest band can be chosen by the coherent
state $|\alpha_{0}\rangle$, and the constructed Bloch states are exact
eigenstates that are exactly degenerate (flat band). The single-photon loss
only alternates the Bloch eigenstates inside the band:

$\displaystyle\hat{a}|\psi_{l,m}\rangle$ $\displaystyle=$
$\displaystyle\frac{1}{\sqrt{q}}\sum_{p=0}^{q-1}e^{imp\frac{2\pi}{q}}\hat{a}\hat{R}_{q}^{p}|\alpha_{0}\rangle$
(54) $\displaystyle=$
$\displaystyle\frac{1}{\sqrt{q}}\sum_{p=0}^{q-1}e^{imp\frac{2\pi}{q}}\hat{R}_{q}^{p}\hat{R}_{q}^{\dagger
p}\hat{a}\hat{R}_{q}^{p}|\alpha_{0}\rangle$ $\displaystyle=$
$\displaystyle\frac{\alpha_{0}}{\sqrt{q}}\sum_{p=0}^{q-1}e^{i(m-1)p\frac{2\pi}{q}}\hat{R}_{q}^{p}|\alpha_{0}\rangle$
$\displaystyle\propto$ $\displaystyle|\psi_{l,m-1}\rangle.$

In the third line, we have used the property
$\hat{R}_{q}^{\dagger}\hat{a}\hat{R}_{q}=\hat{a}e^{-i\frac{2\pi}{q}}$ and
$\hat{a}|\alpha_{0}\rangle=\alpha_{0}|\alpha_{0}\rangle$. Thus, photon decay
does not leak quantum information out of the code subspaces Gottesman et al.
(2001); Puri et al. (2019); Rymarz et al. (2021); Conrad (2021), allowing for
fault tolerant syndrome detection that exponentially suppresses phase flips
and does not cause error backaction on the encoded system Puri et al. (2019).

##  Appendix F Sharp-boundary elliptical well in phase space

In this section, we calculate analytical expression of the engineered driving
potential for the elliptical Hamiltonian in phase space. We start from the
general formula of the engineered driving potential as follows

$\displaystyle
V(x,t)=\frac{1}{2}\int_{-\infty}^{+\infty}|k|f_{T}(k,t)e^{ikx}dk.$ (55)

According to the identity
$x^{-2}=-\frac{1}{2}\int_{-\infty}^{+\infty}|k|e^{ikx}dk$ and defining the
following function

$\displaystyle
v(x,t)\equiv\frac{1}{2\pi}\int_{-\infty}^{+\infty}f_{T}(k,t)e^{ikx}dk,$ (56)

the driving potential $V(x,t)$ is the convolution of $-x^{-2}$ and $u(x,t)$,
i.e.,

$\displaystyle V(x,t)$ $\displaystyle=$ $\displaystyle-x^{-2}*v(x,t)$ (57)
$\displaystyle=$
$\displaystyle-\int_{-\infty}^{+\infty}\frac{1}{z^{2}}v(x-z,t)dz$
$\displaystyle=$ $\displaystyle-\lim_{\epsilon\to
0^{+}}\int_{-\infty}^{+\infty}\Re\Big{[}\frac{1}{(z-i\epsilon)^{2}}\Big{]}v(x-z,t)dz$
$\displaystyle=$ $\displaystyle-\lim_{\epsilon\to
0^{+}}\int_{-\infty}^{+\infty}\frac{z^{2}-\epsilon^{2}}{(z^{2}+\epsilon^{2})^{2}}v(x-z,t)dz.$

Here, we should first replace the convolution function $1/z^{2}$ by
${(z^{2}-\epsilon^{2})}/{(z^{2}+\epsilon^{2})^{2}}$ to get converged integral,
and then take the limit of $\epsilon\to 0$ to obtain $V(x,t)$.

In the $(x,p)$ phase space, we set a new coordinate $(x^{\prime},p^{\prime})$
system which is rotated by an angle $\tau$ between $x$ and $x^{\prime}$ axes
given by the following orthogonal transformation

$\displaystyle x=x^{\prime}\cos\tau-p^{\prime}\sin\tau;\ \ \
p=x^{\prime}\sin\tau+p^{\prime}\cos\tau.$ (58)

We can express the target Hamiltonian Q-function with rotated coordinates by
$H^{{}^{\prime}(T)}_{Q}(x^{\prime},p^{\prime})\equiv H^{(T)}_{Q}(x,p)$. Then,
we project the Hamiltonian on the $x^{\prime}$ axis by the so-called _Radon
transformation_

$\displaystyle\mathcal{R}_{\tau}[H^{{}^{\prime}(T)}_{Q}](x^{\prime})=\int_{-\infty}^{+\infty}H^{{}^{\prime}(T)}_{Q}(x^{\prime},p^{\prime})dp^{\prime}.$
(59)

The 1D Fourier transformation of the above projected function is given by

$\displaystyle\mathcal{F}[\mathcal{R}_{\tau}[H^{{}^{\prime}(T)}_{Q}]](k_{x^{\prime}})=\int\mathcal{R}_{\tau}[H^{{}^{\prime}(T)}_{Q}](x^{\prime})e^{-ik_{x^{\prime}}x^{\prime}}dx^{\prime}.$
(60)

The 2D Fourier transformation of
$H^{{}^{\prime}(T)}_{Q}(x^{\prime},p^{\prime})$ is

$\displaystyle\mathcal{F}[H^{{}^{\prime}(T)}_{Q}](k_{x^{\prime}},k_{p^{\prime}})=\int\int
H^{{}^{\prime}(T)}_{Q}(x^{\prime},p^{\prime})e^{-i(k_{x^{\prime}}x^{\prime}+k_{p^{\prime}}p^{\prime})}dx^{\prime}dp^{\prime}.$
(61)

By plugging Eq. (59) into Eq. (60) and comparing with Eq. (61), we have

$\displaystyle\mathcal{F}[\mathcal{R}_{\tau}[H^{{}^{\prime}(T)}_{Q}]](k_{x^{\prime}})=\mathcal{F}[H^{{}^{\prime}(T)}_{Q}](k_{x^{\prime}},k_{p^{\prime}}=0).$
(62)

This is the so-called projection-slice theorem Bracewell (1956); Ng (2005).
Comparing to Eq. (21) and using the orthogonal transformation (58), we have
$\mathcal{F}[\mathcal{R}_{\tau}[H^{{}^{\prime}(T)}_{Q}]](k)=2\pi
f_{T}(k,\tau)e^{-\frac{\lambda}{4}k^{2}}$ (here we have set
$k_{x^{\prime}}=k$) and thus

$\displaystyle\mathcal{R}_{\tau}[H^{{}^{\prime}(T)}_{Q}](x)=\int^{+\infty}_{-\infty}f_{T}(k,\tau)e^{-\frac{\lambda}{4}k^{2}}e^{ikx}dk.\
$ (63)

In the classical limit $\lambda=0$, according to Eqs. (56) and (57), we have
$v(x,\tau)=\frac{1}{2\pi}\mathcal{R}_{\tau}[H^{{}^{\prime}(T)}_{Q}](x)$ and
thus

$\displaystyle V(x,t)$ $\displaystyle=$
$\displaystyle-\frac{1}{2\pi}\frac{1}{x^{2}}*\mathcal{R}_{\tau=t}[H^{{}^{\prime}(T)}_{Q}](x).$
(64)

Now, we apply the convolutional form (57) and Radon transformation (59) to
calculate the engineered driving potential for the elliptical well in phase
space in the classical limit $\lambda=0$. The boundary of the elliptical
Hamiltonian Q-function in phase space is given by

$\frac{(x^{\prime}\cos\tau-p^{\prime}\sin\tau)^{2}}{a^{2}}+\frac{(x^{\prime}\sin\tau+p^{\prime}\cos\tau)^{2}}{b^{2}}=1.$

Using the transformation (58), we have the following

$\Big{(}\frac{\sin\tau^{2}}{a^{2}}+\frac{\cos\tau^{2}}{b^{2}}\Big{)}p^{\prime
2}-\sin(2\tau)\Big{(}\frac{1}{a^{2}}-\frac{1}{b^{2}}\Big{)}x^{\prime}p^{\prime}+\Big{(}\frac{\cos\tau^{2}}{a^{2}}+\frac{\sin\tau^{2}}{b^{2}}\Big{)}x^{\prime
2}-1=0.$

Using the two solutions $p^{\prime}_{1}$ and $p^{\prime}_{2}$ of the above
equation, the length across the ellipse is

$\displaystyle|p^{\prime}_{1}-p^{\prime}_{2}|$ $\displaystyle=$
$\displaystyle\sqrt{(p^{\prime}_{1}+p^{\prime}_{2})^{2}-4p^{\prime}_{1}p^{\prime}_{2}}$
(65) $\displaystyle=$ $\displaystyle\sqrt{\frac{1}{A^{2}}x^{\prime
2}\sin^{2}2\tau\Big{(}\frac{1}{a^{2}}-\frac{1}{b^{2}}\Big{)}^{2}-\frac{4}{A}\Big{[}x^{\prime
2}\Big{(}\frac{\cos\tau^{2}}{a^{2}}+\frac{\sin\tau^{2}}{b^{2}}\Big{)}-1\Big{]}}$
$\displaystyle=$ $\displaystyle\frac{2}{\sqrt{A}}\sqrt{1-Bx^{\prime 2}}$

where

$\displaystyle
A(\tau)\equiv\frac{\sin\tau^{2}}{a^{2}}+\frac{\cos\tau^{2}}{b^{2}},\ \ \
B(\tau)\equiv\Big{(}\frac{\cos\tau^{2}}{a^{2}}+\frac{\sin\tau^{2}}{b^{2}}\Big{)}-\frac{1}{4A(\tau)}\sin^{2}2\tau\Big{(}\frac{1}{a^{2}}-\frac{1}{b^{2}}\Big{)}^{2}$
(66)

From Eq. (59), we have the Radon transformation

$\displaystyle\mathcal{R}_{\tau}[H^{{}^{\prime}(T)}_{Q}](x^{\prime})=-|p^{\prime}_{1}-p^{\prime}_{2}|=-\frac{2}{\sqrt{A(\tau)}}\sqrt{1-B(\tau)x^{\prime
2}}\ \ \ \mathrm{for}\ \ \ |x^{\prime}|<\frac{1}{\sqrt{B(\tau)}}.\ \ \ \ \ \ \
\ $ (67)

From Eq. (64), the driving potential is given by

$\displaystyle V(x,t)$ $\displaystyle=$
$\displaystyle-\frac{1}{2\pi}x^{-2}*\mathcal{R}_{t}[H^{{}^{\prime}(T)}_{Q}](x)$
(68) $\displaystyle=$ $\displaystyle-\frac{1}{2\pi}\lim_{\epsilon\to
0^{+}}\int_{-\infty}^{+\infty}\frac{z^{2}-\epsilon^{2}}{(z^{2}+\epsilon^{2})^{2}}\mathcal{R}_{t}[H^{{}^{\prime}(T)}_{Q}](x-z)dz$
$\displaystyle=$
$\displaystyle\frac{1}{2\pi}2\sqrt{\frac{B}{A}}\lim_{\epsilon\to
0^{+}}\int_{-\infty}^{+\infty}\frac{z^{2}-\epsilon^{2}}{(z^{2}+\epsilon^{2})^{2}}\sqrt{\big{(}\frac{1}{\sqrt{B}}-x+z\big{)}\big{(}\frac{1}{\sqrt{B}}+x-z\big{)}}dz$
$\displaystyle=$ $\displaystyle-\sqrt{\frac{B}{A}}\lim_{\epsilon\to
0^{+}}\Big{(}1-\Re\big{[}\frac{\epsilon-ix}{\sqrt{B^{-1}+(\epsilon-
ix)^{2}}}\big{]}\Big{)}$

In the limit $\epsilon=0$, we have

$\displaystyle
V(x,t)=\begin{cases}-\sqrt{\frac{B(t)}{A(t)}},&|x|<\frac{1}{\sqrt{B(t)}}\\\
-\sqrt{\frac{B(t)}{A(t)}}\Big{(}1-\frac{|x|}{\sqrt{x^{2}-\frac{1}{B(t)}}}\Big{)},&|x|\geq\frac{1}{\sqrt{B(t)}}.\end{cases}$
(69)

This is the driving potential that can generate classical elliptical potential
with sharp boundary in phase space.

Note that the value of $V(x,t)$ is divergent at $x=\pm{1}/{\sqrt{B(t)}}$
according to Eq. (69). In fact, the Fourier coefficient Eq. (21) is always
divergent for $k\rightarrow\infty$ in the quantum regime $\lambda>0$. For
physical result, we can add an exponentially suppressing factor to the Fourier
coefficient:

$\displaystyle f_{T}(k_{x}.k_{p})\rightarrow
f_{T}(k_{x}.k_{p})e^{-\frac{1}{2}\sigma^{2}(k_{x}^{2}+k_{p}^{2})}\ \
\mathrm{with}\ \ \sigma>\sqrt{{\lambda}/{2}}.$ (70)

According to the convolution theorem Arfken and Weber (2005), we equivalently
modify the Hamiltonian Q-function with a convolution operation, i.e.,

$\displaystyle H^{(T)}_{Q,\sigma}=\int\int
dx^{\prime}dp^{\prime}g(x^{\prime},p^{\prime})H^{(T)}_{Q}(x-x^{\prime},p-p^{\prime})\
\ \mathrm{with}\ \
g(x,p)=\frac{1}{2\pi\sigma^{2}}e^{-\frac{x^{2}+p^{2}}{2\sigma^{2}}}.$ (71)

The kernel function $g(x,p)$ smooths the sharp boundary of the elliptical
well. Correspondingly, the driving potential is also modified with a
convolution, i.e.,

$\displaystyle V_{\sigma}(x,t)=\int\
dx^{\prime}h(x^{\prime})V(x-x^{\prime},t)\ \ \mathrm{with}\ \
h(x)=\frac{1}{\sqrt{\pi(2\sigma^{2}-\lambda)}}e^{-\frac{x^{2}}{2\sigma^{2}-\lambda}}.$
(72)

We point out that, in the Fock representation, the exponentially suppression
factor $e^{-\frac{\lambda}{4}(k^{2}_{x}+k^{2}_{p})}$ in Eq. (29) cancels the
same exponentially increasing factor in Eq. (28). As a result, a sharp well
with boundary narrower than quantum fluctuations ($\sigma<\sqrt{\lambda/2}$)
does exist. However, this scenario cannot be realised by our present method.

##  Appendix G Moiré superlattice

In this section, we discuss how to engineer a Moiré superlattice in phase
space that is formed by two honeycomb phase space lattices overlaid with a
relative twist angle and confined in a finite region with radius $R$,

$\displaystyle
H^{(T)}_{Q}(x,p)=\begin{cases}H_{\theta=0}+H_{\theta=\theta_{0}},&\sqrt{x^{2}+p^{2}}\leq
R\\\ 0,&\sqrt{x^{2}+p^{2}}>R.\end{cases}$ (73)

Here,
$H_{\theta}(x,p)\equiv-\prod_{n=1}^{3}\sin^{2}\big{[}\frac{1}{2}\textbf{{v}}_{n}\cdot\textbf{{z}}(\theta)\big{]}+\frac{3}{32}$
is the honeycomb lattice in phase space Guo et al. (2022). We have defined the
vector
${\textbf{{z}}}(\theta)\equiv(x\cos\theta+p\sin\theta,-x\sin\theta+p\cos\theta)$,
and three ancillary vectors

$\displaystyle\textbf{{v}}_{1}=(\frac{2\sqrt{3}}{3},0),\ \
\textbf{{v}}_{2}=(-\frac{\sqrt{3}}{3},1),\ \
\textbf{{v}}_{3}=(-\frac{\sqrt{3}}{3},-1).$ (74)

In Fig. 3(a) in the main text, we plot the resulting Moiré superlattice
Q-function with twisted angle $\theta_{0}=10^{\circ}$.

We calculate and plot the NcFT coefficient $f_{T}(k_{x},k_{p})$ in Fig. 3(b)
in the main text, which composes of discrete peaks reflecting the discrete
translational symmetry of target Hamiltonian in phase space. The centers of
these peaks take place at
$(k^{q}_{x}=k_{q}\cos\tau_{q},k^{q}_{p}=k_{q}\sin\tau_{q})$ where
$\tau_{q}=q{\pi}/{6},\ q{\pi}/{6}+\theta_{0}$ and $k_{q}={2}/{\sqrt{3}},\ 2,\
{4}/{\sqrt{3}}$ with $q\in\mathbb{Z}$. The finite width of peaks comes from
the boundary condition of Moiré superlattice. In fact, we have the analytical
expression from Eq. (21)

$\displaystyle
f_{T}(k_{x},k_{p})=\sum_{q}\frac{A_{q}}{2\pi}\frac{J_{1}\Big{(}R\sqrt{(k_{x}-k^{q}_{x})^{2}+(k_{p}-k^{q}_{p})^{2}}\Big{)}}{R^{-1}\sqrt{(k_{x}-k^{q}_{x})^{2}+(k_{p}-k^{q}_{p})^{2}}},\
\ $ (75)

where $A_{q}={\pi}/{16},\ -{\pi}/{16},\ {\pi}/{32}$ for $k_{q}={2}{\sqrt{3}},\
2,\ {4}{\sqrt{3}}$ respectively. Due to the long-distance asymptotic behavior
of Bessel function

$\displaystyle J_{1}(k)\approx\sqrt{\frac{2}{\pi k}}\cos(k-\frac{\pi}{4})\ \
\mathrm{for}\ \ |k|\gg 1,$ (76)

we have added an exponentially suppressing factor
$e^{-\frac{\lambda}{4}(k_{x}^{2}+k_{p}^{2})}$ to obtain convergent NcFT
coefficient. As a result, the Hamiltonian Q-function smoothed with a
convolution kernel function
$g(x,p)=\frac{1}{\pi\lambda}e^{-\frac{x^{2}+p^{2}}{\lambda}}$, cf. Eq. (71).

##  Appendix H Hamiltonian operator and Q-function symmetry

In this section, we prove that the Hamiltonian operator and its Q-function has
the same symmetry in phase space. We first discuss the Hamiltonian operator in
the discrete rotational operation $\hat{R}_{q}\equiv
e^{i\hat{a}^{\dagger}\hat{a}\frac{2\pi}{q}}$. From the Fourier form of
Hamiltonian operator Eq. (13), we have

$\displaystyle\hat{R}_{q}H(\hat{x},\hat{p})\hat{R}_{q}^{\dagger}$
$\displaystyle=$ $\displaystyle\frac{\beta}{2\pi}\int\int
dk_{x}dk_{p}f_{T}(k_{x},k_{p})\hat{R}_{q}\exp[i(k_{x}\hat{x}+k_{p}\hat{p})]\hat{R}_{q}^{\dagger}$
(77) $\displaystyle=$ $\displaystyle\frac{\beta}{2\pi}\int\int
dk_{x}dk_{p}f_{T}(k_{x},k_{p})\exp\Big{[}i\big{(}[k_{x}\cos(\frac{2\pi}{q})-k_{p}\sin(\frac{2\pi}{q})]\hat{x}+[k_{x}\sin(\frac{2\pi}{q})+k_{p}\cos(\frac{2\pi}{q})]\hat{p}\big{)}\Big{]}$
$\displaystyle=$ $\displaystyle\frac{\beta}{2\pi}\int\int
dk^{\prime}_{x}dk^{\prime}_{p}f_{T}(k^{\prime}_{x},k^{\prime}_{p})\exp[i(k^{\prime}_{x}\hat{x}+k^{\prime}_{p}\hat{p})].$

Here, we have used the property

$\displaystyle\left\\{\begin{array}[]{lll}\hat{R}_{q}\hat{x}\hat{R}_{q}^{\dagger}=\hat{x}\cos(\frac{2\pi}{q})+\hat{p}\sin(\frac{2\pi}{q})\\\
\hat{R}_{q}\hat{p}\hat{R}_{q}^{\dagger}=-\hat{x}\sin(\frac{2\pi}{q})+\hat{p}\cos(\frac{2\pi}{q}),\\\
\end{array}\right.$ (80)

and transformed the integral coordinates by

$\displaystyle\left\\{\begin{array}[]{lll}k^{\prime}_{x}=k_{x}\cos(\frac{2\pi}{q})-k_{p}\sin(\frac{2\pi}{q})\\\
k^{\prime}_{p}=k_{x}\sin(\frac{2\pi}{q})+k_{p}\cos(\frac{2\pi}{q})\\\
\end{array}\right.$ (83)

with the property $dk_{x}dk_{p}=dk^{\prime}_{x}dk^{\prime}_{p}$ as the
transformation is orthogonal.

Next, we will prove the NcFT coefficient of rotational lattice satisfies
$f(k^{\prime}_{x},k^{\prime}_{p})=f(k_{x},k_{p})$. According to (21), we have

$\displaystyle f_{T}(k^{\prime}_{x},k^{\prime}_{p})$ $\displaystyle=$
$\displaystyle\frac{e^{\frac{\lambda}{4}(k^{\prime 2}_{x}+k^{\prime
2}_{p})}}{2\pi\beta}\int\int
dxdpH^{(T)}_{Q}(x,p)e^{-i(k^{\prime}_{x}x+k^{\prime}_{p}p)}$ (84)
$\displaystyle=$
$\displaystyle\frac{e^{\frac{\lambda}{4}(k^{2}_{x}+k^{2}_{p})}}{2\pi\beta}\int\int
dx^{\prime}dp^{\prime}H^{(T)}_{Q}(x^{\prime},p^{\prime})e^{-i(k_{x}x^{\prime}+k_{p}p^{\prime})},$

where we have used transformation Eq. (83) and made the orthogonal
transformation in the phase space plane

$\displaystyle\left\\{\begin{array}[]{lll}x^{\prime}=x\cos(\frac{2\pi}{q})+p\sin(\frac{2\pi}{q})\\\
p^{\prime}=-x\sin(\frac{2\pi}{q})+p\cos(\frac{2\pi}{q}).\\\
\end{array}\right.$ (87)

Therefore, if the Hamiltonian Q-function satisfied the discrete rotational
symmetry $H^{(T)}_{Q}(x^{\prime},p^{\prime})=H^{(T)}_{Q}(x,p)$, the NcFT
coefficient of rotational lattice satisfies
$f(k_{x},k_{p})=f(k^{\prime}_{x},k^{\prime}_{p})$. Then, comparing Eq. (13)
and Eq. (77), we have
$\hat{R}_{q}H(\hat{x},\hat{p})\hat{R}_{q}^{\dagger}=H(\hat{x},\hat{p})$ and
vice versa .

The same discussion can also be applied for the translational symmetry
described by the displacement operator

$\displaystyle\hat{D}_{\alpha_{0}}H(\hat{x},\hat{p})\hat{D}_{\alpha_{0}}^{\dagger}$
$\displaystyle=$ $\displaystyle\frac{\beta}{2\pi}\int\int
dk_{x}dk_{p}f_{T}(k_{x},k_{p})\hat{D}_{\alpha_{0}}\exp[i(k_{x}\hat{x}+k_{p}\hat{p})]\hat{D}_{\alpha_{0}}^{\dagger}$
(88) $\displaystyle=$ $\displaystyle\frac{\beta}{2\pi}\int\int
dk_{x}dk_{p}f_{T}(k_{x},k_{p})e^{-i(k_{x}x_{0}+k_{p}p_{0})}\exp[i(k_{x}\hat{x}+k_{p}\hat{p})],$

where we have used the property
$\hat{D}_{\alpha_{0}}\hat{a}\hat{D}_{\alpha_{0}}^{\dagger}=\hat{a}-\alpha_{0}$
and thus

$\displaystyle\left\\{\begin{array}[]{lll}\hat{D}_{\alpha_{0}}\hat{x}\hat{D}_{\alpha_{0}}^{\dagger}=\hat{x}-x_{0}\
\mathrm{with}\ x_{0}=\langle\alpha_{0}|\hat{x}|\alpha_{0}\rangle\\\
\hat{D}_{\alpha_{0}}\hat{p}\hat{D}_{\alpha_{0}}^{\dagger}=\hat{p}-p_{0}\
\mathrm{with}\ p_{0}=\langle\alpha_{0}|\hat{p}|\alpha_{0}\rangle.\\\
\end{array}\right.$ (91)

From Eq. (21), we have

$\displaystyle f_{T}(k_{x},k_{p})e^{-i(k_{x}x_{0}+k_{p}p_{0})}$
$\displaystyle=$
$\displaystyle\frac{e^{\frac{\lambda}{4}(k^{2}_{x}+k^{2}_{p})}}{2\pi\beta}\int\int
dxdpH^{(T)}_{Q}(x,p)e^{-i[k_{x}(x+x_{0})+k_{p}(p+p_{0})}$ (92)
$\displaystyle=$
$\displaystyle\frac{e^{\frac{\lambda}{4}(k^{2}_{x}+k^{2}_{p})}}{2\pi\beta}\int\int
dxdpH^{(T)}_{Q}(x+x_{0},p+p_{0})e^{-i(k_{x}x+k_{p}p)}.$

Therefore, if the Hamiltonian Q-function satisfied the translational symmetry
$H^{(T)}_{Q}(x+x_{0},p+p_{0})=H^{(T)}_{Q}(x,p)$, the NcFT coefficient of
rotational lattice satisfies
$f_{T}(k_{x},k_{p})e^{-i(k_{x}x_{0}+k_{p}p_{0})}=f(k_{x},k_{p})$. Then,
comparing Eq. (13) and Eq. (88), we have
$\hat{D}_{\alpha_{0}}H(\hat{x},\hat{p})\hat{D}_{\alpha_{0}}^{\dagger}=H(\hat{x},\hat{p})$
and vice versa .

The same discussion can also be applied for other phase-space symmetries like
the mirror symmetry, i.e., if the Hamiltonian Q-function satisfies
$H^{(T)}_{Q}(x_{0},p)=H^{(T)}_{Q}(\pm x,\pm p)$, the Hamiltonian operator has
$H(\pm\hat{x},\pm\hat{p})=H(\pm\hat{x},\pm\hat{p})$ and vice versa.

Lastly, we point out that the conclusion is also true for the smoothed
Hamiltonian Q-function with a convolution operation, i.e.,

$\displaystyle H^{(T)}_{Q,\sigma}=\int\int
dx^{\prime}dp^{\prime}g(x^{\prime},p^{\prime})H^{(T)}_{Q}(x-x^{\prime},p-p^{\prime})$
(93)

as long as the kernel function is rotationally symmetric, e.g., the standard
Gaussian kernel
$g(x,p)=\frac{1}{2\pi\sigma^{2}}\exp(-\frac{x^{2}+p^{2}}{2\sigma^{2}})$.

## References

  * Gerry and Knight (2004) C. Gerry and P. Knight, _Nonclassical light_ (Cambridge University Press, Cambridge, 2004), pp. 150–194, Introductory Quantum Optics, ISBN 9780521527354\.
  * Strekalov and Leuchs (2019) D. V. Strekalov and G. Leuchs, _Nonlinear Interactions and Non-classical Light_ (Springer International Publishing, Cham, 2019), pp. 51–101, ISBN 978-3-319-98402-5.
  * Kubala et al. (2015a) B. Kubala, V. Gramich, and J. Ankerhold, Physica Scripta T165, 014029 (2015a), URL https://doi.org/10.1088/0031-8949/2015/t165/014029.
  * Braunstein and van Loock (2005) S. L. Braunstein and P. van Loock, Rev. Mod. Phys. 77, 513 (2005), URL https://link.aps.org/doi/10.1103/RevModPhys.77.513.
  * Pan et al. (2012) J.-W. Pan, Z.-B. Chen, C.-Y. Lu, H. Weinfurter, A. Zeilinger, and M. Żukowski, Rev. Mod. Phys. 84, 777 (2012), URL https://link.aps.org/doi/10.1103/RevModPhys.84.777.
  * Yan et al. (2021) Z.-H. Yan, J.-L. Qin, Z.-Z. Qin, X.-L. Su, X.-J. Jia, C.-D. Xie, and K.-C. Peng, Fundamental Research 1, 43 (2021), ISSN 2667-3258, URL https://www.sciencedirect.com/science/article/pii/S2667325820300054.
  * Leghtas et al. (2013) Z. Leghtas, G. Kirchmair, B. Vlastakis, R. J. Schoelkopf, M. H. Devoret, and M. Mirrahimi, Phys. Rev. Lett. 111, 120501 (2013), URL https://link.aps.org/doi/10.1103/PhysRevLett.111.120501.
  * Mirrahimi et al. (2014) M. Mirrahimi, Z. Leghtas, V. V. Albert, S. Touzard, R. J. Schoelkopf, L. Jiang, and M. H. Devoret, New Journal of Physics 16, 045014 (2014), URL https://doi.org/10.1088/1367-2630/16/4/045014.
  * Heeres et al. (2017) R. W. Heeres, P. Reinhold, N. Ofek, L. Frunzio, L. Jiang, M. H. Devoret, and R. J. Schoelkopf, Nature Communications 8, 94 (2017), ISSN 2041-1723, URL https://doi.org/10.1038/s41467-017-00045-1.
  * Rosenblum et al. (2018) S. Rosenblum, P. Reinhold, M. Mirrahimi, L. Jiang, L. Frunzio, and R. J. Schoelkopf, Science 361, 266 (2018), eprint https://www.science.org/doi/pdf/10.1126/science.aat3996, URL https://www.science.org/doi/abs/10.1126/science.aat3996.
  * Flühmann et al. (2019) C. Flühmann, T. L. Nguyen, M. Marinelli, V. Negnevitsky, K. Mehta, and J. P. Home, Nature 566, 513 (2019), ISSN 1476-4687, URL https://doi.org/10.1038/s41586-019-0960-6.
  * Hu et al. (2019) L. Hu, Y. Ma, W. Cai, X. Mu, Y. Xu, W. Wang, Y. Wu, H. Wang, Y. P. Song, C.-L. Zou, et al., Nature Physics 15, 503 (2019), ISSN 1745-2481, URL https://doi.org/10.1038/s41567-018-0414-3.
  * Campagne-Ibarcq et al. (2020) P. Campagne-Ibarcq, A. Eickbusch, S. Touzard, E. Zalys-Geller, N. E. Frattini, V. V. Sivak, P. Reinhold, S. Puri, S. Shankar, R. J. Schoelkopf, et al., Nature 584, 368 (2020), ISSN 1476-4687, URL https://doi.org/10.1038/s41586-020-2603-3.
  * Gertler et al. (2021) J. M. Gertler, B. Baker, J. Li, S. Shirol, J. Koch, and C. Wang, Nature 590, 243 (2021), ISSN 1476-4687, URL https://doi.org/10.1038/s41586-021-03257-0.
  * Cochrane et al. (1999) P. T. Cochrane, G. J. Milburn, and W. J. Munro, Phys. Rev. A 59, 2631 (1999), URL https://link.aps.org/doi/10.1103/PhysRevA.59.2631.
  * Gottesman et al. (2001) D. Gottesman, A. Kitaev, and J. Preskill, Phys. Rev. A 64, 012310 (2001), URL https://link.aps.org/doi/10.1103/PhysRevA.64.012310.
  * Travaglione and Milburn (2002) B. C. Travaglione and G. J. Milburn, Phys. Rev. A 66, 052322 (2002), URL https://link.aps.org/doi/10.1103/PhysRevA.66.052322.
  * Michael et al. (2016) M. H. Michael, M. Silveri, R. T. Brierley, V. V. Albert, J. Salmilehto, L. Jiang, and S. M. Girvin, Phys. Rev. X 6, 031006 (2016), URL https://link.aps.org/doi/10.1103/PhysRevX.6.031006.
  * Albert et al. (2018) V. V. Albert, K. Noh, K. Duivenvoorden, D. J. Young, R. T. Brierley, P. Reinhold, C. Vuillot, L. Li, C. Shen, S. M. Girvin, et al., Phys. Rev. A 97, 032346 (2018), URL https://link.aps.org/doi/10.1103/PhysRevA.97.032346.
  * Grimsmo et al. (2020) A. L. Grimsmo, J. Combes, and B. Q. Baragiola, Phys. Rev. X 10, 011058 (2020), URL https://link.aps.org/doi/10.1103/PhysRevX.10.011058.
  * Tzitrin et al. (2020) I. Tzitrin, J. E. Bourassa, N. C. Menicucci, and K. K. Sabapathy, Phys. Rev. A 101, 032315 (2020), URL https://link.aps.org/doi/10.1103/PhysRevA.101.032315.
  * Terhal et al. (2020) B. M. Terhal, J. Conrad, and C. Vuillot, Quantum Science and Technology 5, 043001 (2020), URL https://doi.org/10.1088/2058-9565/ab98a5.
  * Joshi et al. (2021) A. Joshi, K. Noh, and Y. Y. Gao, Quantum Science and Technology 6, 033001 (2021), URL https://doi.org/10.1088/2058-9565/abe989.
  * Cai et al. (2021) W. Cai, Y. Ma, W. Wang, C.-L. Zou, and L. Sun, Fundamental Research 1, 50 (2021), ISSN 2667-3258, URL https://www.sciencedirect.com/science/article/pii/S2667325820300145.
  * Krastanov et al. (2015) S. Krastanov, V. V. Albert, C. Shen, C.-L. Zou, R. W. Heeres, B. Vlastakis, R. J. Schoelkopf, and L. Jiang, Phys. Rev. A 92, 040303 (2015), URL https://link.aps.org/doi/10.1103/PhysRevA.92.040303.
  * Fösel et al. (2020) T. Fösel, S. Krastanov, F. Marquardt, and L. Jiang, _Efficient cavity control with snap gates_ (2020), URL https://arxiv.org/abs/2004.14256.
  * Kudra et al. (2022) M. Kudra, M. Kervinen, I. Strandberg, S. Ahmed, M. Scigliuzzo, A. Osman, D. P. Lozano, M. O. Tholén, R. Borgani, D. B. Haviland, et al., PRX Quantum 3, 030301 (2022), URL https://link.aps.org/doi/10.1103/PRXQuantum.3.030301.
  * Puri et al. (2019) S. Puri, A. Grimm, P. Campagne-Ibarcq, A. Eickbusch, K. Noh, G. Roberts, L. Jiang, M. Mirrahimi, M. H. Devoret, and S. M. Girvin, Phys. Rev. X 9, 041009 (2019), URL https://link.aps.org/doi/10.1103/PhysRevX.9.041009.
  * Rymarz et al. (2021) M. Rymarz, S. Bosco, A. Ciani, and D. P. DiVincenzo, Phys. Rev. X 11, 011032 (2021), URL https://link.aps.org/doi/10.1103/PhysRevX.11.011032.
  * Conrad (2021) J. Conrad, Phys. Rev. A 103, 022404 (2021), URL https://link.aps.org/doi/10.1103/PhysRevA.103.022404.
  * Zaslavskii et al. (1986) G. M. Zaslavskii, M. Y. Zakharov, R. Z. Sagdeev, D., Usikov, and A. A. Chernikov, Zh. Eksp. Teor. Fiz. 91, 500 (1986).
  * Berman et al. (1991) G. P. Berman, V. Y. Rubaev, and G. M. Zaslavsky, Nonlinearity 4, 543 (1991), URL https://dx.doi.org/10.1088/0951-7715/4/2/015.
  * Leboeuf et al. (1990) P. Leboeuf, J. Kurchan, M. Feingold, and D. P. Arovas, Phys. Rev. Lett. 65, 3076 (1990), URL https://link.aps.org/doi/10.1103/PhysRevLett.65.3076.
  * Leboeuf et al. (1992) P. Leboeuf, J. Kurchan, M. Feingold, and D. P. Arovas, Chaos: An Interdisciplinary Journal of Nonlinear Science 2, 125 (1992), eprint https://doi.org/10.1063/1.165915, URL https://doi.org/10.1063/1.165915.
  * Billam and Gardiner (2009) T. P. Billam and S. A. Gardiner, Phys. Rev. A 80, 023414 (2009), URL https://link.aps.org/doi/10.1103/PhysRevA.80.023414.
  * Guo et al. (2013) L. Guo, M. Marthaler, and G. Schön, Phys. Rev. Lett. 111, 205303 (2013), URL https://link.aps.org/doi/10.1103/PhysRevLett.111.205303.
  * Zhang et al. (2017) Y. Zhang, J. Gosner, S. M. Girvin, J. Ankerhold, and M. I. Dykman, Phys. Rev. A 96, 052124 (2017), URL https://link.aps.org/doi/10.1103/PhysRevA.96.052124.
  * Liang et al. (2018) P. Liang, M. Marthaler, and L. Guo, New Journal of Physics 20, 023043 (2018), URL https://doi.org/10.1088/1367-2630/aaa7c3.
  * Lörch et al. (2019) N. Lörch, Y. Zhang, C. Bruder, and M. I. Dykman, Phys. Rev. Research 1, 023023 (2019), URL https://link.aps.org/doi/10.1103/PhysRevResearch.1.023023.
  * Guo et al. (2022) L. Guo, V. Peano, and F. Marquardt, Phys. Rev. B 105, 094301 (2022), URL https://link.aps.org/doi/10.1103/PhysRevB.105.094301.
  * Rahav et al. (2003) S. Rahav, I. Gilary, and S. Fishman, Phys. Rev. A 68, 013820 (2003), URL https://link.aps.org/doi/10.1103/PhysRevA.68.013820.
  * Goldman and Dalibard (2014) N. Goldman and J. Dalibard, Phys. Rev. X 4, 031027 (2014), URL https://link.aps.org/doi/10.1103/PhysRevX.4.031027.
  * Bukov et al. (2015a) M. Bukov, L. D’Alessio, and A. Polkovnikov, Advances in Physics 64, 139 (2015a), eprint https://doi.org/10.1080/00018732.2015.1055918, URL https://doi.org/10.1080/00018732.2015.1055918.
  * Casas et al. (2001) F. Casas, J. A. Oteo, and J. Ros, Journal of Physics A: Mathematical and General 34, 3379 (2001), URL https://doi.org/10.1088/0305-4470/34/16/305.
  * Blanes et al. (2009) S. Blanes, F. Casas, J. Oteo, and J. Ros, Physics Reports 470, 151 (2009), ISSN 0370-1573, URL https://www.sciencedirect.com/science/article/pii/S0370157308004092.
  * Eckardt and Anisimovas (2015) A. Eckardt and E. Anisimovas, New Journal of Physics 17, 093039 (2015), URL https://doi.org/10.1088/1367-2630/17/9/093039.
  * Mikami et al. (2016) T. Mikami, S. Kitamura, K. Yasuda, N. Tsuji, T. Oka, and H. Aoki, Phys. Rev. B 93, 144307 (2016), URL https://link.aps.org/doi/10.1103/PhysRevB.93.144307.
  * Bukov et al. (2015b) M. Bukov, L. D’Alessio, and A. Polkovnikov, Advances in Physics 64, 139 (2015b), eprint https://doi.org/10.1080/00018732.2015.1055918, URL https://doi.org/10.1080/00018732.2015.1055918.
  * Rudner and Lindner (2020) M. S. Rudner and N. H. Lindner, Nature Reviews Physics 2, 229 (2020), ISSN 2522-5820, URL https://doi.org/10.1038/s42254-020-0170-z.
  * Scully and Zubairy (1997) M. Scully and M. Zubairy, _Quantum Optics_ , Quantum Optics (Cambridge University Press, 1997), ISBN 9780521435956, URL https://books.google.de/books?id=20ISsQCKKmQC.
  * Marthaler and Dykman (2006) M. Marthaler and M. I. Dykman, Phys. Rev. A 73, 042108 (2006), URL https://link.aps.org/doi/10.1103/PhysRevA.73.042108.
  * Peano et al. (2012) V. Peano, M. Marthaler, and M. I. Dykman, Phys. Rev. Lett. 109, 090401 (2012), URL https://link.aps.org/doi/10.1103/PhysRevLett.109.090401.
  * Dykman (2012) M. Dykman, _Fluctuating Nonlinear Oscillators: From Nanomechanics to Quantum Superconducting Circuits_ (Oxford University Press, 2012), ISBN 9780199691388, URL https://doi.org/10.1093/acprof:oso/9780199691388.001.0001.
  * Bachtold et al. (2022) A. Bachtold, J. Moser, and M. I. Dykman, Rev. Mod. Phys. 94, 045005 (2022), URL https://link.aps.org/doi/10.1103/RevModPhys.94.045005.
  * Bistritzer and MacDonald (2011) R. Bistritzer and A. H. MacDonald, Proceedings of the National Academy of Sciences 108, 12233 (2011), eprint https://www.pnas.org/doi/pdf/10.1073/pnas.1108174108, URL https://www.pnas.org/doi/abs/10.1073/pnas.1108174108.
  * Cao et al. (2018a) Y. Cao, V. Fatemi, A. Demir, S. Fang, S. L. Tomarken, J. Y. Luo, J. D. Sanchez-Yamagishi, K. Watanabe, T. Taniguchi, E. Kaxiras, et al., Nature 556, 80 (2018a), ISSN 1476-4687, URL https://doi.org/10.1038/nature26154.
  * Cao et al. (2018b) Y. Cao, V. Fatemi, S. Fang, K. Watanabe, T. Taniguchi, E. Kaxiras, and P. Jarillo-Herrero, Nature 556, 43 (2018b), ISSN 1476-4687, URL https://doi.org/10.1038/nature26160.
  * Thomson et al. (2018) A. Thomson, S. Chatterjee, S. Sachdev, and M. S. Scheurer, Phys. Rev. B 98, 075109 (2018), URL https://link.aps.org/doi/10.1103/PhysRevB.98.075109.
  * Guo et al. (2016) L. Guo, M. Liu, and M. Marthaler, Phys. Rev. A 93, 053616 (2016), URL https://link.aps.org/doi/10.1103/PhysRevA.93.053616.
  * Guo (2021) L. Guo, in _Phase Space Crystals_ (IOP Publishing, 2021), 2053-2563, pp. 4–1 to 4–28, ISBN 978-0-7503-3563-8, URL https://dx.doi.org/10.1088/978-0-7503-3563-8ch4.
  * Sacha (2015) K. Sacha, Scientific Reports 5, 10787 (2015), ISSN 2045-2322, URL https://doi.org/10.1038/srep10787.
  * Mierzejewski et al. (2017) M. Mierzejewski, K. Giergiel, and K. Sacha, Phys. Rev. B 96, 140201 (2017), URL https://link.aps.org/doi/10.1103/PhysRevB.96.140201.
  * Giergiel et al. (2018) K. Giergiel, A. Miroszewski, and K. Sacha, Phys. Rev. Lett. 120, 140401 (2018), URL https://link.aps.org/doi/10.1103/PhysRevLett.120.140401.
  * Sacha (2020) K. Sacha, _Condensed Matter Physics in the Time Dimension_ (Springer International Publishing, Cham, 2020), pp. 173–235, ISBN 978-3-030-52523-1, URL https://doi.org/10.1007/978-3-030-52523-1_5.
  * Guo and Liang (2020) L. Guo and P. Liang, New Journal of Physics 22, 075003 (2020), URL https://doi.org/10.1088/1367-2630/ab9d54.
  * Hannaford and Sacha (2022) P. Hannaford and K. Sacha, AAPPS Bulletin 32, 12 (2022), ISSN 2309-4710, URL https://doi.org/10.1007/s43673-022-00041-8.
  * Moritz et al. (2003) H. Moritz, T. Stöferle, M. Köhl, and T. Esslinger, Phys. Rev. Lett. 91, 250402 (2003), URL https://link.aps.org/doi/10.1103/PhysRevLett.91.250402.
  * Hadzibabic et al. (2004) Z. Hadzibabic, S. Stock, B. Battelier, V. Bretin, and J. Dalibard, Phys. Rev. Lett. 93, 180403 (2004), URL https://link.aps.org/doi/10.1103/PhysRevLett.93.180403.
  * Chen et al. (2014) F. Chen, J. Li, A. D. Armour, E. Brahimi, J. Stettenheim, A. J. Sirois, R. W. Simmonds, M. P. Blencowe, and A. J. Rimberg, Phys. Rev. B 90, 020506 (2014), URL https://link.aps.org/doi/10.1103/PhysRevB.90.020506.
  * Hofheinz et al. (2011) M. Hofheinz, F. Portier, Q. Baudouin, P. Joyez, D. Vion, P. Bertet, P. Roche, and D. Esteve, Phys. Rev. Lett. 106, 217005 (2011), URL https://link.aps.org/doi/10.1103/PhysRevLett.106.217005.
  * Chen et al. (2011) F. Chen, A. J. Sirois, R. W. Simmonds, and A. J. Rimberg, Applied Physics Letters 98, 132509 (2011).
  * Armour et al. (2013) A. D. Armour, M. P. Blencowe, E. Brahimi, and A. J. Rimberg, Phys. Rev. Lett. 111, 247001 (2013), URL https://link.aps.org/doi/10.1103/PhysRevLett.111.247001.
  * Gramich et al. (2013) V. Gramich, B. Kubala, S. Rohrer, and J. Ankerhold, Phys. Rev. Lett. 111, 247002 (2013), URL https://link.aps.org/doi/10.1103/PhysRevLett.111.247002.
  * Leppäkangas et al. (2013) J. Leppäkangas, G. Johansson, M. Marthaler, and M. Fogelström, Phys. Rev. Lett. 110, 267004 (2013), URL https://link.aps.org/doi/10.1103/PhysRevLett.110.267004.
  * Leppäkangas et al. (2015) J. Leppäkangas, M. Fogelström, A. Grimm, M. Hofheinz, M. Marthaler, and G. Johansson, Phys. Rev. Lett. 115, 027004 (2015), URL https://link.aps.org/doi/10.1103/PhysRevLett.115.027004.
  * Leppäkangas et al. (2016) J. Leppäkangas, M. Fogelström, M. Marthaler, and G. Johansson, Phys. Rev. B 93, 014506 (2016), URL https://link.aps.org/doi/10.1103/PhysRevB.93.014506.
  * Armour et al. (2015) A. D. Armour, B. Kubala, and J. Ankerhold, Phys. Rev. B 91, 184508 (2015), URL https://link.aps.org/doi/10.1103/PhysRevB.91.184508.
  * Trif and Simon (2015) M. Trif and P. Simon, Phys. Rev. B 92, 014503 (2015), URL https://link.aps.org/doi/10.1103/PhysRevB.92.014503.
  * Kubala et al. (2015b) B. Kubala, V. Gramich, and J. Ankerhold, Physica Scripta T165, 014029 (2015b), URL https://doi.org/10.1088/0031-8949/2015/t165/014029.
  * Hofer et al. (2016) P. P. Hofer, J.-R. Souquet, and A. A. Clerk, Phys. Rev. B 93, 041418 (2016), URL https://link.aps.org/doi/10.1103/PhysRevB.93.041418.
  * Dambach et al. (2017) S. Dambach, B. Kubala, and J. Ankerhold, New Journal of Physics 19, 023027 (2017), URL https://doi.org/10.1088/1367-2630/aa5bb6.
  * Lang and Armour (2021) B. Lang and A. D. Armour, New Journal of Physics 23, 033021 (2021), URL https://doi.org/10.1088/1367-2630/abe483.
  * Lang et al. (2022) B. Lang, G. F. Morley, and A. D. Armour, _Discrete time translation symmetry breaking in a josephson junction laser_ (2022), URL https://arxiv.org/abs/2208.03191.
  * Guo and Marthaler (2016) L. Guo and M. Marthaler, New Journal of Physics 18, 023006 (2016), URL https://doi.org/10.1088/1367-2630/18/2/023006.
  * Srivastava and Choi (2012) H. Srivastava and J. Choi, in _Zeta and q-Zeta Functions and Associated Series and Integrals_ , edited by H. Srivastava and J. Choi (Elsevier, London, 2012), pp. 1–140, ISBN 978-0-12-385218-2, URL https://www.sciencedirect.com/science/article/pii/B9780123852182000013.
  * Bracewell (1956) R. N. Bracewell, Australian Journal of Physics 9, 198 (1956), URL https://doi.org/10.1071/PH560198.
  * Ng (2005) R. Ng, ACM Trans. Graph. 24, 735–744 (2005), ISSN 0730-0301, URL https://doi.org/10.1145/1073204.1073256.
  * Arfken and Weber (2005) G. B. Arfken and H. J. Weber, _Mathematical Methods for Physicists_ (Academic Press; 6th edition, 2005), ISBN 978-0120598762.