# Decoherence Limits the Cost to Simulate an Anharmonic Oscillator

Tzula B. Propp  Center for Quantum Information & Control, University of New
Mexico, Albuquerque, NM 87131, USA  Sayonee Ray  Center for Quantum
Information & Control, University of New Mexico, Albuquerque, NM 87131, USA
Department of Physics and Astronomy, University of Waterloo, Ontario, N2L 3G1,
Canada  Current affiliation: IonQ Inc, 4505 Campus Dr, College Park, MD 20740,
USA  John B. DeBrota  Center for Quantum Information & Control, University of
New Mexico, Albuquerque, NM 87131, USA  Tameem Albash  Center for Quantum
Information & Control, University of New Mexico, Albuquerque, NM 87131, USA
Department of Electrical and Computer Engineering, University of New Mexico,
Albuquerque, NM 87131, USA  Department of Physics and Astronomy, University of
New Mexico, Albuquerque, NM 87131, USA  Ivan Deutsch  Center for Quantum
Information & Control, University of New Mexico, Albuquerque, NM 87131, USA
Department of Physics and Astronomy, University of New Mexico, Albuquerque, NM
87131, USA

###### Abstract

We study how decoherence increases the efficiency with which we can simulate
the quantum dynamics of an anharmonic oscillator, governed by the Kerr effect.
As decoherence washes out the fine-grained subPlanck structure associated with
phase-space quantum interference in the closed quantum system, the open
quantum dynamics can be more efficiently simulated using a coarse-grained
finite-difference numerical integration. We tie this to the way in which
decoherence recovers the semiclassical truncated Wigner approximation (TWA),
which strongly differs from the exact closed-system dynamics at times when
quantum interference leads to cat states and more general superpositions of
coherent states. The regression to semiclassical dynamics become more
pronounced as the initial amplitude of the oscillator grows, with implications
for the quantum advantage that might be accessible in noisy devices.

##  I Introduction

The macroscopic world is largely described by classical statistical physics
even though the underlying fundamental description is quantum mechanical. As
emphasized in the seminal work of Zurek Zurek (1981, 1991); Zurek and Paz
(1994); Anglin  _et al._ (1997); Zurek (1998); Paz and Zurek (1999);
Karkuszewski  _et al._ (2002); Zurek (2001), decoherence helps us understand
the transition from the quantum-to-classical world: classical states are more
robust to the coupling to the environment whereas highly nonclassical states,
such as macroscopic superposition “cat states” Dodonov  _et al._ (1974), are
fragile in the face of decoherence Leggett  _et al._ (1987); Caldeira and
Leggett (1983); Kim and Bužek (1992); Zurek (2001).

This observation not only explains why quantum effects are largely unobserved
in the classical world, but it also represents the fundamental challenge of
large scale quantum information processing as decoherence limits the quantum
complexity one can harness at the macroscopic scale. While in principle one
can tame decoherence through quantum error correction Shor (1995); Gottesman
_et al._ (2001); Albert and Faist , we do not yet have the means to do so
fault-tolerantly Shor (1996); Aharonov and Ben-Or (1997); Kitaev (1997); Knill
_et al._ (1998) except at very small scales et al. (2023).

Nevertheless, one may hope to achieve a meaningful quantum advantage without
fault tolerant error correction. A key question, thus, is how much complexity
one can expect to harness in a noisy-intermediate scale quantum (NISQ)
Preskill (2018) device given that decoherence largely washes out the most
nonclassical features.

One way to quantify the quantum-to-classical transition is using the phase
space representation of quantum mechanics Curtright  _et al._ (2013). The
nonclassical features of quantum states are reflected in the negativity of the
Wigner quasiprobability distribution Cahill and Glauber (1969); Hudson (1974);
Kenfack and Życzkowski (2004). This is associated with the fine-grained
subPlanck-scale structure that arises in quantum dynamics beyond what is
described by classical flow in phase space Zurek (2001). These properties
challenge the classical simulation of macroscopic quantum dynamics; direct
numerical propagation of the quasiprobability function scales poorly when
trying to capture ultra-fine features Heller (1976), and applying techniques
such as the methods of characteristics Weedbrook  _et al._ (2012) or Monte
Carlo sampling Sellier  _et al._ (2015) is difficult when the quasiprobability
function develops substantial negativity Welland and Ferry (2020).
Decoherence, however, washes out subPlanck structure and the associated
negativity by introducing a time scale over which the coherence of a quantum
state is lost Habib  _et al._ (1998). This transition occurs more rapidly as
the system becomes more macroscopic Zurek (2001). Thus, the phase space
representation exhibits the trade off between robustness and complexity — more
subPlanck structure and negativity lead to faster decoherence, while those
features that are robust to decoherence are essentially classical and more
efficiently simulatable.

In this work we consider how the phase space representation quantifies the
robustness versus complexity trade off in the context of quantum simulation.
Instead of focusing on negativity itself (which precludes semiclassical
sampling), we will study how decoherence simplifies numerical simulation by
erasing fine-grained structure. We consider a canonical toy model – an
anharmonic oscillator such as a single bosonic mode evolving under a Kerr
nonlinearity. While a trivial model, it exhibits a wide range of well-studied
nonclassical features, such as squeezing at short times Milburn (1986);
Milburn and Holmes (1986), the collapse and revival of quantum oscillations
Polkovnikov (2010), and the generation of superpositions of coherent states
Miranowicz  _et al._ (1990); Tara  _et al._ (1993); van Enk (2003, 2005);
Stobińska  _et al._ (2008). These effects have been studied in a variety of
platforms including the seminal proposal in quantum optics Yurke and Stoler
(1986) and studies in trapped ions Stobińska  _et al._ (2011) as well as the
realizations in atomic BEC Greiner  _et al._ (2002) and the pioneering
experiments in circuit QED Kirchmair  _et al._ (2013) which now form the basis
of quantum error-correcting codes that have been realized Grimm  _et al._
(2020). While there are closed form solutions for this integrable model Tara
_et al._ (1993); van Enk (2003); Sudheesh  _et al._ (2004), including for the
open quantum system Milburn and Holmes (1986); Peinová and Luk (1990);
Chaturvedi and Srinivasan (1991a, b); McDonald and Clerk (2022, 2023), in
practice determining the expectation values of some observables can be
numerically intensive when the mean number of bosons is deep in the
macroscopic regime. We emphasize that formally nothing is computationally
complex in this model, as all quantities can be extracted with algorithms that
scale polynomially with the mean boson number. Nonetheless, there remains a
tradeoff between robustness and ease of numerical simulation. Our goal is to
study how decoherence makes this system more “classical” and in doing so
determine which observables correspond to properties that are fragile to
decoherence, and which observables can be more efficiently simulated for the
open quantum system.

The remainder of this article is structured as follows. In Sec. II we derive
the dynamics of the quantum state under pure unitary evolution, both in
discrete and continuous variable descriptions. In Sec. III we expand on the
continuous variable dynamics to include open system effects using the Fokker-
Planck equations. In Sec. IV we study the time evolution of moments of the
Wigner function which capture fine-grained information and discuss the cost of
simulating the open quantum system. We conclude with a summary and outlook in
Sec. V.

##  II Closed system dynamics

We consider an anharmonic oscillator for a single bosonic mode with a Kerr
nonlinearity, governed by the Hamiltonian

$\hat{H}=\frac{\kappa}{2}\hat{a}^{{\dagger}2}\hat{a}^{2}=\frac{\kappa}{2}\hat{n}(\hat{n}-1).$
(1)

Here and throughout $\hbar=1$. This Hamiltonian arises, e.g., in quantum
optics in the presence of an intensity-dependent index of refraction (Kerr
effect) Walls and Milburn (2008), and in atom optics for bosonic atoms
undergoing cold collisions Steel  _et al._ (1998). For concreteness, we will
consider the bosons here to be photons. The dynamics of this system has been
well studied, and we review the salient results here.

Figure 1: The Wigner function for an initial coherent state with amplitude
$\alpha_{0}=4$ at $t=0$ evolving under the Kerr interaction. At short times we
see characteristic squeezing by at longer times we see negativity in the
Wigner function and the formation of intricate subPlanck-scale structure. At
times $\kappa t_{M,N}=(M/N)2\pi$, $M$ and $N$ coprime, the state is a
superposition of $N$ coherent states, Eq. (II). For $N\gg\alpha_{0}$, the
coherent states are not distinguishable and we observe structures like the
short time “banana” with minimal negativity in the top row. At later times,
distinguishable superpositions form. These include high-order kitten states
which may appear multiple times (e.g. the $N=9$, $M=2$ kitten state in the
center of the third row) as well as the Schrödinger cat state at
$t=\pi/\kappa$ (bottom right corner). Right before a distinguishable kitten
state forms, we observe shuriken-like Wigner functions (bottom left corner).
In between $t=\pi/\kappa$ and $t=2\pi/\kappa$, the evolution is symmetric to
the evolution between $t=0$ and $t=\pi/\kappa$ as the state progresses through
the same states in reverse (not pictured for brevity). An animation depicting
the continuous dynamics is given in Propp (2023a).

As photon number is conserved here, the Heisenberg equation of motion for the
annihilation operator $d\hat{a}/dt=-i\kappa\hat{a}^{\dagger}\hat{a}\hat{a}$ is
integrable, $\hat{a}(t)=e^{-i\kappa t\hat{a}^{\dagger}\hat{a}}\hat{a}(0)$. In
the classical (mean field) limit, $\hat{a}(t)\rightarrow\alpha_{c}(t)$,
$\alpha_{c}(t)=e^{-i\kappa t|\alpha_{c}|^{2}}\alpha_{c}(0)$, describing the
rotation of the phasor at an angular rate proportional to the amplitude
squared (the classical Kerr effect). An arbitrary pure state evolves according
to the trivial solution

$\left|\psi(t)\right\rangle=\sum_{n=0}^{\infty}c_{n}e^{-i\kappa
tn(n-1)/2}\left|n\right\rangle.$ (2)

We will be particularly interested in the case that the initial state is a
coherent state $\left|\alpha_{0}\right\rangle$, in which case
$c_{n}=e^{-\frac{|\alpha_{0}|^{2}|}{2}}\alpha_{0}^{n}/\sqrt{n!}$.

The evolution of the state is periodic returning to its initial condition when
$\kappa t=M\,2\pi$, for integer $M$. There are additional revivals due to
quantum interference at times $\kappa t_{M,N}=(M/N)2\pi$, where $M$ and $N$
are co-prime integers. To see this consider the discrete functions,
$f_{e}[n]\equiv e^{-i\frac{M\pi n^{2}}{N}}$ and $f_{o}[n]\equiv
e^{-i\frac{M\pi n(n-1)}{N}}$. These are periodic with
$f_{e(o)}[n]=f_{e(o)}[n+N]$ for $N$ even (odd) and thus can be expanded as
discrete Fourier transforms,

$\displaystyle f_{e(o)}[n]$ $\displaystyle=$
$\displaystyle\frac{1}{\sqrt{N}}\sum_{k=0}^{N-1}c_{e(o)}[k]e^{i\frac{2\pi
nk}{N}},$ (3) $\displaystyle c_{e(o)}[k]$ $\displaystyle=$
$\displaystyle\frac{1}{\sqrt{N}}\sum_{n=0}^{N-1}f_{e(o)}[n]e^{-i\frac{2\pi
nk}{N}},$ (4)

where the coefficients $c_{e(o)}[k]$ implicitly depend on both $M$ and $N$. It
follows that the time evolution operator at times $t_{M,N}$,
$\hat{U}(t_{M,N})=e^{-i\frac{M\pi}{N}\hat{n}(\hat{n}-1)}$, is a superposition
of phase space rotation operators

$\displaystyle\hat{U}(t_{M,N})$ $\displaystyle=$
$\displaystyle\frac{1}{\sqrt{N}}\sum_{k=0}^{N-1}c_{e}[k]e^{i\frac{(2k+1)\pi}{N}\hat{n}},\;\;N\;\mathrm{even},$
(5) $\displaystyle=$
$\displaystyle\frac{1}{\sqrt{N}}\sum_{k=0}^{N-1}c_{o}[k]e^{i\frac{2k\pi}{N}\hat{n}},\;\;N\;\mathrm{odd}.$
(6)

Therefore, at the discrete times $t_{M,N}$, the state revives to coherent
superposition of $N$ copies of the initial state, distributed in phase space
on the circle with radius $|\alpha_{0}|$ at multiples of $2\pi/N$. We can find
the expansion coefficients by noting that

$\displaystyle c_{e}[k]$ $\displaystyle=$
$\displaystyle\frac{1}{\sqrt{N}}\sum_{n=0}^{N-1}f_{e}[n+Mk]e^{i\frac{\pi
Mk^{2}}{N}}$ (7) $\displaystyle=$ $\displaystyle c_{e}[0]e^{i\frac{\pi
Mk^{2}}{N}},$ $\displaystyle c_{o}[k]$ $\displaystyle=$
$\displaystyle\frac{1}{\sqrt{N}}\sum_{n=0}^{N-1}f_{o}[n+Mk]e^{i\frac{\pi
Mk(k+1)}{N}}$ (8) $\displaystyle=$ $\displaystyle c_{o}[0]e^{i\frac{\pi
Mk(k+1)}{N}},$

where the second equality follows from the periodicity of $f_{e(o)}[n]$.

We thus find for an initial coherent state $\left|\alpha_{0}\right\rangle$.
Using Parseval’s theorem, $|c_{e(o)}[0]|^{2}=1$, and neglecting the overall
phase $c_{e(o)}[0]$, the state at the revival time $t_{M,N}$ is

$\displaystyle\left|\psi(t_{M,N})\right\rangle$ $\displaystyle=$
$\displaystyle\begin{cases}\frac{1}{\sqrt{N}}\sum_{k=0}^{N-1}e^{i\frac{\pi
Mk^{2}}{N}}\left|\alpha_{0}e^{i\frac{(2k+M)\pi}{N}}\right\rangle&\,N\;\mathrm{even},\\\
\frac{1}{\sqrt{N}}\sum_{k=0}^{N-1}e^{i\frac{\pi
Mk(k+1)}{N}}\left|\alpha_{0}e^{i\frac{2k\pi}{N}}\right\rangle&\,N\;\mathrm{odd}.\end{cases}$

For the case $M=1$, $N=2$, we obtain
$\left|\psi(t_{1,2})\right\rangle=\frac{1}{\sqrt{2}}(\left|i\alpha_{0}\right\rangle+i\left|-i\alpha_{0}\right\rangle)$,
a “cat state.” For general $N$ and $M$, the state is a superposition of $N$
coherent states as originally analyzed by Bialynicka-Birula Bialynicka-Birula
(1968), which we will denote “kitten states” 111The kitten states and the cat
state generated by the Kerr interaction all have this remarkable feature: each
phase is such that the overall normalization constant is independent of the
coherent state amplitude $\alpha_{0}$.. These states form a dense set and
exist at all $t_{M,N}$, but only come to the fore when the coherent states are
distinguishably separated.

We are particularly interested in studying the dynamics of the system using
phase space representations. The equation of motion for the Wigner function,
$W(\alpha,\alpha^{*},t)$, with a quadratic nonlinearity is

$\displaystyle\frac{\partial W}{\partial t}$ $\displaystyle=$
$\displaystyle\\{H_{W},W\\}_{MP}$ $\displaystyle=$
$\displaystyle-i\left(\frac{\partial H_{W}}{\partial\alpha}\frac{\partial
W}{\partial\alpha^{*}}-\frac{\partial H_{W}}{\partial\alpha^{*}}\frac{\partial
W}{\partial\alpha}\right)$
$\displaystyle-\frac{i}{8}\left(\frac{\partial^{3}H_{W}}{\partial^{2}\alpha^{*}\partial\alpha}\frac{\partial^{3}W}{\partial^{2}\alpha\partial\alpha^{*}}-\frac{\partial^{3}H_{W}}{\partial^{2}\alpha\partial\alpha^{*}}\frac{\partial^{3}W}{\partial^{2}\alpha^{*}\partial\alpha}\right),$

where
$H_{W}(\alpha,\alpha^{*})=\frac{\kappa}{2}|\alpha|^{2}(|\alpha|^{2}-2)+\frac{\kappa}{4}$
is the Weyl-symbol of the Hamiltonian and $\\{H_{W},W\\}_{MP}$ is the Moyal
bracket Polkovnikov (2013) . The first line in Eq. (II) is
$\\{H_{W},W\\}_{PB}$, the Poisson bracket, corresponding to the classical flow
on phase space generated by $H_{W}$. The second line represents the
nonclassical dynamics, which can lead to negativity of the Wigner function.
These terms become nonnegligible when the Wigner function develops subPlanck-
scale structure as this corresponds to a rapidly varying function whose higher
order derivatives are at least order one. Neglecting these quantum dynamics is
known as the truncated Wigner approximation (TWA) Polkovnikov (2013), which
well describes the dynamics for times short times compared to the Ehrenfest
time (that is, when corrections to the Poisson bracket become nonnegligible).

Substituting $H_{W}$ for the Kerr Hamiltonian into Eq. (II), we have

$\displaystyle\frac{\partial W}{\partial t}$ $\displaystyle=$
$\displaystyle-i\kappa(|\alpha|^{2}-1)\left(\alpha^{*}\frac{\partial
W}{\partial\alpha^{*}}-\alpha\frac{\partial W}{\partial\alpha}\right)$ (11)
$\displaystyle-\frac{i\kappa}{4}\left(\alpha\frac{\partial^{3}W}{\partial^{2}\alpha\partial\alpha^{*}}-\alpha^{*}\frac{\partial^{3}W}{\partial^{2}\alpha^{*}\partial\alpha}\right).$

The first line is the TWA representing classical flow on phase space. Note
$-i(\alpha^{*}\partial_{\alpha^{*}}-\alpha\partial_{\alpha})=X\partial_{P}-P\partial_{X}$
is the rotation operator on phase space, where the quadratures are defined by
$\alpha=(X+iP)/\sqrt{2}$. The TWA has the expected form of a Kerr effect, that
is, a rotation at angular rate $\kappa|\alpha|^{2}$ for $|\alpha|$ large. The
evolution according to the TWA can be solved by the method of characteristics,
$W_{\mathrm{TWA}}(\alpha,\alpha^{*},t)=W_{\mathrm{TWA}}(\alpha(-t),\alpha^{*}(-t),0)$,
where $\alpha(t)=e^{-i\kappa t(|\alpha|^{2}-1)}\alpha$ is the classical flow.
The TWA evolution of an initial coherent state represented as a Gaussian
wavepacket is thus,
$W_{\mathrm{TWA}}(\alpha,\alpha^{*},t)=\frac{1}{\pi}\exp\\{-2|e^{i\kappa
t(|\alpha|^{2}-1)}\alpha-\alpha_{0}|^{2}\\}$, as shown in the top row of Fig.
2. For short times the nonlinear rotation “shears” the distribution, leading
to a squeezed Gaussian. At longer times the distribution becomes stretched and
becomes highly nonGaussian, but remains a positive probability distribution in
the TWA. At these longer times the corrections to the TWA become
nonneglegible. Figure 1 shows the Wigner function evolution for
$\alpha_{0}=4$, calculated using the numerical methods described in Appendix
A. For very short times, the TWA evolution approximately matches the exact
evolution, but negativity soon develops and revivals occur at the expected
times, yielding kitten states and subPlanck structure.

While in principle the exact state is available in any representation, in
practice extracting measured values of observables becomes numerically
intensive when $\alpha_{0}$ is sufficiently large, which may describe
realistic mean photon numbers, e.g.,
$\langle\hat{n}\rangle=|\alpha_{0}|^{2}=10^{6}$. We will specifically focus on
symmetrically-ordered correlation functions, which are moments of the Wigner
function,

$\langle\\{\hat{X}^{n}\hat{P}^{m}\\}_{\mathrm{sym}}\rangle(t)=\int
dXdP\;X^{n}P^{m}\;W(X,P,t),$ (12)

where the symmetrically-ordered product of $\hat{X}$ and $\hat{P}$ is to be
normalized by $n+m$, e.g.,
$\\{\hat{X}\hat{P}\\}_{\mathrm{sym}}=\frac{1}{2}(\hat{X}\hat{P}+\hat{P}\hat{X})$.
Such correlation functions arise, e.g., in considering moments of a quadrature
$\hat{X}_{\theta}=\cos\theta\hat{X}+\sin\theta\hat{P}$ optically measured in
homodyne detection, $\langle\hat{X}^{n}_{\theta}\rangle(t)=\int
dX_{\theta}\mathscr{P}_{\theta}(X_{\theta},t)X^{n}_{\theta}$.

For the closed system, the probability distribution which determines these
moments $\langle\hat{X}^{n}_{\theta}\rangle(t)$ is the square of the wave
function as a function of the quadrature $X_{\theta}$ eigenvalue, which has an
analytic form

$\displaystyle\mathscr{P}(X_{\theta},t)=|\langle
X_{\theta}|\psi(t)\rangle|^{2}$ (13)
$\displaystyle=\left|\sum_{n=0}^{\infty}e^{-\frac{|\alpha_{0}|^{2}}{2}}\frac{\alpha_{0}^{n}}{\sqrt{n!}}e^{-i\left(\frac{\kappa
t}{2}+\theta\right)}e^{-i\frac{\kappa
tn^{2}}{2}}u_{n}(X_{\theta})\right|^{2},$

where we have defined
$u_{n}(X_{\theta})=A_{n}H_{n}(X_{\theta})e^{-\frac{X_{\theta}^{2}}{2}}$, with
$A_{n}=\left(\sqrt{\pi}2^{n}n!\right)^{-1/2}$ and $H_{n}(X_{\theta})$ the
$n$th Hermite polynomial. The relevant number of terms to evaluate this formal
expression grows as $\sim|\alpha_{0}|^{4}$. More generally, this probability
distribution can be obtained by marginalizing the Wigner function

$\mathscr{P}(X_{\theta},t)=\int dP_{\theta}W(X_{\theta},P_{\theta},t)$ (14)

where we define $P_{\theta}=X_{\theta+\pi/2}$.

The higher order moments of the marginal capture the fine-grained structure in
the Wigner function, and generically this is hardest to simulate in the large
$|\alpha_{0}|$ limit, when nonclassical subPlanck-scale structure develops.
However, as discussed above, this fine grained structure is washed out by
decoherence. Thus, our goal is to determine how much of this complexity
remains in the open quantum system and how the reduction of this complexity
leads to more efficient simulations.

The Kerr anharmonic oscillator is a useful system for benchmarking the
nonclassical dynamics and comparing numerical simulations to the exact
analytic solution. The conservation of photon number implies normally ordered
correlation functions have a simple closed form for an initial coherent state,

$\langle\hat{a}^{\dagger\mu}\hat{a}^{\nu}\rangle(t)\equiv\langle\alpha_{0}|\hat{a}^{{\dagger}\mu}(t)\hat{a}^{\nu}(t)|\alpha_{0}\rangle,$
(15)

where $\mu$ and $\nu$ are nonnegative integers. Using
$\hat{a}^{\nu}(t)\left|\alpha_{0}\right\rangle=\alpha_{0}^{\nu}e^{-i\frac{\kappa
t}{2}\nu(\nu-1)}\left|\alpha_{0}e^{-i\nu\kappa t}\right\rangle$ and
$\langle\alpha|\beta\rangle=e^{-\frac{1}{2}|\alpha-\beta|^{2}}e^{-\frac{1}{2}(\alpha\beta^{*}-\alpha^{*}\beta)}$
one finds

$\displaystyle\langle\hat{a}^{\dagger\mu}\hat{a}^{\nu}\rangle(t)=(\alpha_{0}^{*})^{\mu}(\alpha_{0})^{\nu}e^{i\frac{\kappa
t}{2}[\mu(\mu-1)-\nu(\nu-1)]}$
$\displaystyle\times\;e^{-|\alpha_{0}|^{2}(1-\cos[(\mu-\nu)\kappa
t])}e^{i|\alpha_{0}|^{2}\sin[(\mu-\nu)\kappa t]},$ (16)

which tends towards zero exponentially with $|\alpha_{0}|^{2}$ except at
special times when $\kappa t=M(\mu-\nu)\pi$, with $M$ an integer.
Understanding quantum dynamics at these “recurrences” will be the focus of
Sec. IV.

In contrast, no such simple form exists for direct calculation of correlation
functions associated with the symmetrically ordered expectation values.
Symmetrically ordered expectation values are related to the normally ordered
ones according to

$\displaystyle\left\langle\\{\hat{a}^{\dagger\mu}\hat{a}^{\nu}\\}_{\mathrm{sym}}\right\rangle$
$\displaystyle=$ $\displaystyle\sum_{\gamma=0}^{{\rm
min}(\mu,\nu)}\frac{\mu!\nu!}{\gamma!(\mu-\gamma)!(\nu-\gamma)!2^{\gamma}}\left\langle\hat{a}^{\dagger\mu-\gamma}\,\hat{a}^{\nu-\gamma}\right\rangle.$

For very large order moments with large $|\alpha_{0}|$, accurate calculation
of such sums becomes numerically intensive.

##  III Open System dynamics

We consider the simplest model of decoherence for the single mode - the damped
simple harmonic oscillator due to photon loss, whereby the state evolves
according to the master equation

$\frac{\partial\hat{\rho}}{\partial
t}=-i[\hat{H},\hat{\rho}]+\mathscr{L}[\hat{\rho}],$ (18)

where

$\mathscr{L}[\hat{\rho}]=-\frac{\gamma}{2}\left(\hat{a}^{\dagger}\hat{a}\hat{\rho}+\hat{\rho}\hat{a}^{\dagger}\hat{a}\right)+\gamma\hat{a}\rho\hat{a}^{\dagger}$
(19)

is the Lindbladian describing photon loss in a zero temperature reservoir. In
the absence of Hamiltonian evolution, for a general superposition of coherent
states, $\left|\psi\right\rangle=\sum_{k}c_{k}\left|\alpha_{k}\right\rangle$
we have the closed form solution Walls and Milburn (2008)

$\displaystyle
e^{\mathscr{L}t}[\left|\psi\right\rangle\left\langle\psi\right|]$
$\displaystyle=$
$\displaystyle\sum_{k}|c_{k}|^{2}\left|\alpha_{k}e^{-\frac{\gamma
t}{2}}\right\rangle\left\langle\alpha_{k}e^{-\frac{\gamma t}{2}}\right|$
$\displaystyle+\sum_{k\neq
k^{\prime}}c_{k}c^{*}_{k^{\prime}}\left\langle\alpha_{k^{\prime}}|\alpha_{k}\right\rangle^{1-e^{-\gamma
t}}\left|\alpha_{k^{\prime}}e^{-\frac{\gamma
t}{2}}\right\rangle\left\langle\alpha_{k}e^{-\frac{\gamma t}{2}}\right|.$

The amplitude of coherent states in the mixture decay as expected, but the
coherences in the superposition of coherent states decay more rapidly, at a
rate that depends on their overlap (distinguishability). In particular, for a
cat state, for $t\ll 1/\gamma$,

$\left\langle-\alpha_{0}|\alpha_{0}\right\rangle^{1-e^{-\gamma
t}}\approx\left\langle-\alpha_{0}|\alpha_{0}\right\rangle^{\gamma
t}=e^{-2|\alpha_{0}|^{2}\gamma t},$ (21)

and this coherence decays at the rate $2|\alpha_{0}|^{2}\gamma$. For
macroscopic $\alpha_{0}$ this implies that the coherence between the coherent
states is lost essentially instantaneously compared to the time of substantial
energy loss.

The Wigner function for the cat state thus evolves under the damping channel
as

$\displaystyle W(X,P,t)$ $\displaystyle=$
$\displaystyle\frac{1}{2}W_{0}(X,P-X_{0}e^{-\frac{\gamma t}{2}})$
$\displaystyle+\frac{1}{2}W_{0}(X,P+X_{0}e^{-\frac{\gamma
t}{2}})+e^{-X_{0}^{2}\gamma t}\sin(2XX_{0})W_{0}(X,P),$

where $W_{0}(X,P)=e^{-(X^{2}+P^{2})}/\sqrt{\pi}$ is the Wigner function of the
vacuum. The subPlanck-scale structure decays rapidly, at the rate
$X_{0}^{2}\gamma=2\left\langle\hat{n}\right\rangle\gamma$. This can also be
seen directly from the phase space dynamics. Recall the Weyl representation of
the Lindbladian

$\displaystyle\mathscr{L}[\hat{\rho}]_{W}$ $\displaystyle=$
$\displaystyle\frac{\gamma}{2}\left(\frac{\partial}{\partial\alpha}\alpha+\frac{\partial}{\partial\alpha^{*}}\alpha^{*}\right)W+\frac{\gamma}{2}\frac{\partial^{2}W}{\partial\alpha\partial\alpha^{*}}$
$\displaystyle=\frac{\gamma}{2}\left(\frac{\partial}{\partial
X}X+\frac{\partial}{\partial
P}P\right)W+\frac{\gamma}{4}\left(\frac{\partial^{2}W}{\partial
X^{2}}+\frac{\partial^{2}W}{\partial P^{2}}\right),$

which is a Fokker-Planck equation. The first order derivative “drift” terms
generate the decay of energy of the coherent state. The Laplacian diffusion
terms lead to the rapid washing out of subPlanck scale structure, and thus
decoherence. Indeed, we see the action of diffusion on the interference term
in Eq. (III) gives, to leading order in $X_{0}$,
$\frac{\gamma}{4}(\frac{\partial^{2}}{\partial
X^{2}}+\frac{\partial^{2}}{\partial P^{2}})\sin(2XX_{0})W_{0}(X,P)\approx-
X_{0}^{2}\gamma\sin(2XX_{0})W_{0}(X,P)$, which shows that diffusion leads to
decays of coherence with the same rate we obtained by different methods. The
finer the subPlanck scale structure, the faster decoherence washes it out
Zurek (2001) – a manifestation of the robustness vs. complexity tradeoff.

Figure 2: The Wigner functions of the closed system (top row) and open system
(bottom two rows), evolving from an initial coherent state with
$\alpha_{0}=4$, are plotted at times corresponding to the first six-kitten
state (first column), the cat state (middle column), and the second six-kitten
state (last column). In each inset, both the exact quantum state and the state
calculated using the truncated Wigner approximation (TWA) are shown on the
left and right, respectively. For nonzero coupling to the environment, the
exact state more closely resembles that given by the TWA, with the resemblance
increasing with both decoherence strength $\gamma$ (vertical axis) and time
(horizontal axis).

The equation of motion for the Wigner function under the concurrent action of
the Kerr Hamiltonian and photon loss is given by $\frac{\partial W}{\partial
t}=\\{H,W\\}_{MB}+\mathscr{L}[\hat{\rho}]_{W}$. While a formal solution exists
for the Husimi representation Milburn and Holmes (1986), no such solution
exists for the Wigner function, and generally one must resort to numerical
integration. As first studied by Stobinska et al. Stobińska  _et al._ (2008),
given the invariance of the Kerr interaction under rotation in phase space,
this is best done in polar coordinates by expressing the complex amplitude
$\alpha\equiv re^{i\phi}$, giving the equation of motion

$\displaystyle\frac{\partial}{\partial t}W(r,\phi,t)$ $\displaystyle=$
$\displaystyle\\!\left[\kappa(r^{2}-1)\frac{\partial}{\partial\phi}-\frac{\kappa}{16}\left(\frac{1}{r}\frac{\partial}{\partial
r}+\frac{\partial^{2}}{\partial
r^{2}}+\frac{1}{r^{2}}\frac{\partial^{2}}{\partial\phi^{2}}\right)\frac{\partial}{\partial\phi}\right]W(r,\phi,t)$
(24) $\displaystyle+$
$\displaystyle\left[\gamma\left(1+\frac{r}{2}\frac{\partial}{\partial
r}\right)+\frac{\gamma}{8}\left(\frac{1}{r}\frac{\partial}{\partial
r}+\frac{\partial^{2}}{\partial
r^{2}}+\frac{1}{r^{2}}\frac{\partial^{2}}{\partial\phi^{2}}\right)\right]W(r,\phi,t).$

The first line is the closed-system evolution, where the first term is the
expected classical rotation in phase phase depending on the amplitude squared,
$r^{2}$ (the TWA evolution), and the second term is the nonclassical
Hamiltonian flow. The second line is the Fokker-Planck equation where the
first term is the damping (drift radially inward to the origin), and the
second term is diffusion.

Closed form solutions for the Wigner function are generally not available.
Recently McDonald and Clerk developed a powerful method based on
diagonalization of the Lindblad superoperator, which follows from fundamental
dissipative symmetry present in all quadratic bosonic Lindbladians McDonald
and Clerk (2023). Using this they were able to obtain a closed form propagator
in phase space, and a solution for the Wigner function in terms of a compact
but infinite sum. We seek to probe the role of dissipation in reducing
computational complexity associated with subPlanck-scale structure. For this
we look to the the approach of Stobinska et al., who used a finite difference
method in order to numerically integrate this partial differential equation
(PDE) Stobińska  _et al._ (2008). This method discretized the Wigner funcion
in phase space but the numerics were limited to small $\alpha_{0}$ as the
required grid size would otherwise grow too large. This is a reflection of the
need to capture the fine grained structure that develops in the Wigner
function. However, as noted, we expect the fine-grained structure to be
limited by diffusion (decoherence), and thus a coarse grained approximation to
the Wigner function should give a good approximation for the open quantum
system. Similarly, since decoherence is washing out the features generated by
nonclassical flow, we expect the TWA to give a better representation of the
state in the presence of decoherence as seen in Fig. 2. We study both of these
quantitatively in the section to follow.

##  IV Quadrature Moments

We study in detail the behavior of expectation values of powers of the phase
quadrature operator

$\langle\hat{X}^{n}_{\theta}\rangle(t)=\sum_{m=0}^{n}{n\choose
m}\cos^{n-m}\theta\sin^{m}\theta\langle\\{\hat{X}^{n-m}\hat{P}^{m}\\}_{\mathrm{sym}}\rangle(t).$
(25)

These describe statistics of homodyne measurements, with higher-order moments
corresponding to finer-scale features of the Wigner function, and with their
evolution governed by the quantum dynamics. As discussed above, the
distinction between the quantum and classical dynamics is due to quantum
interference between different photons numbers. For an initial distribution of
field amplitudes, as for an initial coherent state, under classical dynamics
the quadrature moments will collapse as a function of time due to the
nonlinear phase shifts, whereas under closed-system unitary dynamics these
moments also exhibit a series of revivals at recurrence times due to quantum
interference, as discussed above. Away from these times, the expectation value
is well-described by the semiclassical dynamics given by the TWA Polkovnikov
(2013).

We seek to understand the deviation of the exact quantum dynamics from the
semiclassical dynamics given by the TWA. In this section, we first derive the
existence and times $t^{(n)}_{M,N}$ of these recurrences for the $n^{th}$
phase-quadrature moment. Then, we show why the deviations from the TWA occur
_solely_ at times near these recurrences (with the effect becoming more
prominent in the large-$\alpha_{0}$ limit), despite the dense set of highly-
quantum kitten states throughout the evolution. Finally, we study the quantum
behavior of $\langle\hat{X}^{n}_{\theta}\rangle(t)$ itself in terms of its
deviation from semiclassical dynamics.

The existence of recurrence times can be seen from the expression for normally
ordered correlation functions for the closed system given an initial coherent
state with amplitude $\alpha_{0}$ in Eq. (II). We observe that the magnitude
of $\langle\hat{a}^{{\dagger}\mu}\hat{a}^{\nu}\rangle$ has periodicity (that
is, recurrence) at times

$\displaystyle T_{\mu\nu}=\frac{2\pi}{|\mu-\nu|\kappa}.$ (26)

$T_{\mu\mu}$ is undefined, which is consistent with the Kerr interaction being
number preserving; there is never a recurrence for powers of the number
operator because it does not change during the evolution. Since the
expectation value of any normally ordered operator can be expressed as a sum
of expectation values of the form of Eq. (II), any normally ordered operator
will also exhibit recurrences at these times. Furthermore, the period
$T_{\mu\nu}$ of an expectation value
$\left\langle\hat{a}^{\dagger\mu}\hat{a}^{\nu}\right\rangle(t)$ is unchanged
by operator ordering since the recurrence depends only on the difference
$\mu-\nu$ and this difference is preserved by the commutation relations.
Because of this, a symmetrized expectation value
$\left\langle\\{\hat{a}^{\dagger\mu}\hat{a}^{\nu}\\}_{\rm
sym}\right\rangle(t)$ also has the same period $T_{\mu\nu}$.

The dependence of the recurrence times on the difference $\mu-\nu$ and not
$\mu$ and $\nu$ separately is also true for the open quantum system governed
by Eq. (18). The open system recurrences also exhibit independence of operator
ordering, so that normally ordered and symmetrically ordered operators with
the same difference $\mu-\nu$ have the same recurrence times. We derive these
results below with further details given in Appendix B.

The Lindblad master equation Eq. (18) expressed in the number basis is given
by

$\displaystyle\frac{d}{dt}\rho_{nm}(t)=-\left(\frac{\gamma}{2}(n+m)+i\frac{\kappa}{2}\left((n^{2}-n)-(m^{2}-m)\right)\right)\rho_{nm}(t)+\gamma\sqrt{(n+1)(m+1)}\rho_{n+1,\,m+1}(t).$
(27)

We can write a formal solution for the density matrix elements as a
generalized discrete Fourier series

$\displaystyle\rho_{nm}(t)=\sum_{n^{\prime},m^{\prime}}A_{n^{\prime}m^{\prime}}^{nm}e^{-\tilde{\gamma}_{n^{\prime}m^{\prime}}t}\delta_{n^{\prime}-m^{\prime},n-m},$
(28)

where we define a complex function

$\displaystyle\tilde{\gamma}_{n^{\prime}m^{\prime}}=\frac{\gamma}{2}(n^{\prime}+m^{\prime})+i\frac{\kappa}{2}\left((n^{\prime
2}-n^{\prime})-(m^{\prime 2}-m^{\prime})\right).$ (29)

For the closed system, $\tilde{\gamma}_{n^{\prime}m^{\prime}}$ is purely
imaginary, and the Fourier coefficients $A_{n^{\prime}m^{\prime}}^{nm}$ are
all zero except for $n^{\prime}=m,m^{\prime}=m$ with
$A_{nm}^{nm}=\rho_{nm}(0)$. The sum thus collapses to a single term. For the
open system, all Fourier coefficients with $n^{\prime}<n$ and $m^{\prime}<m$
are zero, and the coefficients $A_{n^{\prime}m^{\prime}}^{nm}$ can be solved
recursively by truncating the Fock space (for details see Appendix B).

In all cases the form of Eq. (28) implies that each set of elements
$\rho_{nm}(t)$ of the density matrix on a diagonal “stripe,” with $n-m$ a
fixed constant, evolves independently. This follows directly for the “weak
symmetry” associated with the Lindbladian on the dissipative Kerr oscillator,
recently identified by McDonald and Clerk McDonald and Clerk (2022). We can
use this fact to study the nature of the recurrences in the open quantum
system. It follows from Eq. (28) that

$\displaystyle\left\langle\hat{a}^{\dagger\mu}\hat{a}^{\nu}\right\rangle(t)$
$\displaystyle=\sum\limits_{nm}P_{nm}^{\mu\nu}\rho_{nm}(t)\delta_{n-m,\mu-\nu}$
$\displaystyle=\sum_{nm}R_{nm}^{\mu\nu}e^{-\tilde{\gamma}_{nm}t}\delta_{n-m,\mu-\nu},$
(30)

where $P_{nm}^{\mu\nu}$ is a positive real combinatorial factor and where we
have defined new generalized Fourier coefficients

$\displaystyle
R_{nm}^{\mu\nu}=\sum_{n^{\prime}}P_{n-m+n^{\prime},n^{\prime}}^{\mu\nu}A_{nm}^{n-m+n^{\prime},n^{\prime}}.$
(31)

From this form of the expectation value, we make the following three
conclusions:

  1. 1.

The recurrence times of operators depend only on the difference, $\mu-\nu$,
and not the powers $\mu$ and $\nu$ directly.

  2. 2.

The recurrence times are unchanged by operator ordering. For symmetricly
ordered and normally ordered expectation values the combinatorial factor
$P_{nm}^{\mu\nu}$ will be different, but otherwise the expressions for the
expectation value in Eq. (IV) are the same.

  3. 3.

Open system effects will not change the recurrence times as the correlation
function depends on the same frequency components that contribute in the
closed quantum system.

Thus, the periodicity of the expectation value is independent of operator
ordering, and of being for an open or closed system, up to rescaling of those
coefficients by the real part of $\tilde{\gamma}_{nm}$ and a change of the
value of the combinatorial prefactor; both of these leave the recurrence time
unchanged.222Note that this argument fails for $\nu=\mu$ since these depend on
the diagonal of the density matrix, where $\tilde{\gamma}_{jk}$ is completely
real. However, this is to be expected, as these expectation values are simply
powers of the number operator and exhibit no periodicity/recurrence in the
first place. This behavior extends to any anharmonic, single-mode bosonic
system at finite temperature as we show in Appendix B.

Figure 3: Deviation of the expectation value
$\langle\hat{X}_{\theta}^{6}\rangle(t)$ from its value calculated under the
TWA is plotted for the closed quantum system evolving under the Kerr
interaction (that is,
$\delta\langle\hat{X}_{\theta}^{6}\rangle(t)\equiv\langle\hat{X}_{\theta}^{6}\rangle(t)-\langle\hat{X}_{\theta}^{6}\rangle_{\rm
TWA}(t)$). The initial state is a coherent state with $\alpha_{0}=6$, and we
normalize the expectation value by $\alpha_{0}^{6}$. Different values of
$\theta$ correspond to different phase quadratures and give the same times of
recurrence. The $\theta$ are chosen to correspond to angles present in the $6$
kitten state, which form at $\kappa t=\pi/3$. Plots of the Wigner function at
times $\kappa t=\pi/18,\,2\pi/9,\pi/3$, and $2\pi/5$ are included as insets.
Notably, there is a dense set of times throughout evolution where the Wigner
function is a kitten state with highly quantum features, yet the TWA only
fails to give the correct expectation value for $\hat{X}_{\theta}^{6}$ near
the beginning of evolution and near the first recurrence time $t^{(6)}_{1,6}$.
At the times $\kappa t=2\pi/9$ and $\kappa t=2\pi/5$, we do not observe a
deviation from the TWA due to a mismatch between the symmetry of the
measurement operator $\hat{X}_{\theta}^{6}$ (which has the $6$-fold symmetry
of the cyclic group $\mathbb{Z}_{6}$) and the symmetries of the $5$ and $9$
kitten states ($\mathbb{Z}_{5}$ and $\mathbb{Z}_{9}$, respectively). Near the
recurrence time at $\kappa t=\pi/3$, there are additional contributions from
high order kittens such that the terms of order
$\mathscr{O}(\left\langle\alpha_{0}|\alpha_{0}\,e^{i\frac{2\pi}{N}}\right\rangle)$
in the last line of Eq. (IV) become non-negligible, which are exponentially
suppressed in the large-$\alpha_{0}$ limit.

To understand why away from these recurrences the system is well-described by
the TWA, consider the normally ordered expectation value
$\left\langle\psi(t_{M,N})\right|\hat{a}^{\dagger\mu}\hat{a}^{\nu}\left|\psi(t_{M,N})\right\rangle$,
where $\left|\psi(t_{M,N})\right\rangle$ is the $N$ kitten state at time
$t_{M,N}$ defined in Eq. (II). This expectation value is only sensitive to the
number of kittens $N$ when $\mu-\nu=p\,N$ with $p$ a non-zero integer. To see
this, note from its form that
$\hat{a}^{N}\left|\psi(t_{M,N})\right\rangle=\pm\alpha_{0}^{N}\left|\psi(t_{M,N})\right\rangle$,
with the minus sign occurring for $N$-even and $M$-odd only333This is due to
an overall rotation in the location of kitten states when $N$ is even and $M$
is odd.. This property follows from the considerations. Up to normalization
and an overall phase, the effect of $\hat{a}$ on an $N$-kitten state is to
permute the phases of the coherent states in the superposition, forming a
cyclic group ($\mathbb{Z}_{N}$) Nash (2016), so that all phases return to
their original values after $N$ applications Propp (2023b). Because of this
property, these states were originally denoted “generalized coherent states.”
Bialynicka-Birula (1968). The same is true for the action of
$\hat{a}^{\dagger}$ on $\left\langle\psi(t_{M,N})\right|$. Only when $\hat{a}$
and $\hat{a}^{\dagger}$ are applied a number of times whose difference is an
integer multiple of $N$ do the phases fail to average out to zero, giving rise
to recurrence at the times predicted by the period formula in Eq. (26). The
recurrence becomes more sharply peaked at this time only in the
large-$\alpha_{0}$ limit, when the coherent states are approximately
orthogonal.

Figure 4: Averaged-deviation from the TWA ($\overline{\delta}^{(n)}(t)$) for
$\alpha_{0}=3$ (left) and $\alpha_{0}=6$ (right) as defined in Eq. (35). In
both cases, the effect of increased decoherence is to suppress quantum
deviations from the TWA. Increasing the amplitude of the initial coherent
state both magnifies the effects of decoherence and reduces the initial short-
time deviation from the TWA, and the oscillations become more sharply peaked
around the recurrence times $t^{(n)}_{M,N}$ Polkovnikov (2010, 2003).

We derive these properties for general kitten states (not only those generated
by the Kerr interaction) in Appendix C. The behavior of the expectation value
is given by

$\displaystyle\left\langle\psi(t_{M,N})\right|\hat{a}^{\dagger\mu}\hat{a}^{\nu}\left|\psi(t_{M,N})\right\rangle=\begin{cases}\big{\langle}\hat{n}\big{\rangle}^{\mu}&\,\mu=\nu,\\\
\alpha_{0}^{*\mu-\nu}\big{\langle}\hat{n}\big{\rangle}^{\nu}&\,\mu-\nu=p\,N,\,p\in\mathbb{Z}^{\pm},\\\
0+\mathscr{O}\bigr{(}\big{\langle}\alpha_{0}\big{|}\alpha_{0}e^{i\frac{2\pi}{N}}\big{\rangle}\bigr{)}&\mu-\nu=p\,N,\,p\notin\mathbb{Z},\end{cases}$

where $\mathbb{Z}^{\pm}$ is the set of non-zero integers, and $\mathbb{Z}$ is
the set of all integers including zero. Note that, for normally ordered
expectation values of a _mixture_ of $N$ coherent states, the same expression
holds except that the last case becomes exactly zero for all $\alpha_{0}$; it
is the superposition states that create oscillations around the recurrence
times for finite $\alpha_{0}$.

We thus find that the effects of quantum coherence, as seen in the correlation
functions, depend solely on the difference $\mu-\nu$, which is preserved by
operator reordering. In particular, for the symmetrically ordered moments
under consideration here, when this difference appropriately aligns with
recurrence times we expect strong deviations from the TWA. At other times, the
lack of alignment leads to cancellation making phase relationships
approximately irrelevant and we expect the TWA to be a good approximation.
Furthermore, as $\alpha_{0}$ increases, quantum deviations from the TWA at
times near (but not exactly at) recurrences $T_{\mu\nu}$ become exponentially
suppressed as the coherent states in the final case of Eq. (IV) become
increasingly orthogonal. From our earlier expression for the normally ordered
expectation values in the closed system, Eq. (II), we also see that this
suppression is exponential in time, with a dependence of
$e^{-|\alpha_{0}|^{2}(t-T_{\mu\nu})^{2}}$. The quantum deviations from the TWA
become more suppressed as $|\alpha_{0}|$ increases, including the initial
short-time deviation from the TWA Polkovnikov (2010, 2003) and at the
recurrence times, which become increasingly sharp.

We now have the full set of tools needed to consider the expectation value of
a moment of a quadrature operator $\langle\hat{X}_{\theta}^{n}\rangle$. Any
such operator is decomposable into a sum of symmetrically ordered operators

$\displaystyle\langle\hat{X}_{\theta}^{n}\rangle=\sum_{i}^{n}C_{i,n,\theta}\langle\\{\hat{a}^{\dagger
n-i}\,\hat{a}^{i}\\}_{\mathrm{sym}}\rangle,$ (33)

where the index $i$ runs along the even integers for $n$-even, and along the
odd integers for $n$-odd. As we established earlier, for every symmetrically
ordered expectation value
$\langle\\{\hat{a}^{\dagger\mu}\hat{a}^{\nu}\\}_{\mathrm{sym}}\rangle$, there
is a set of recurrence times $pT_{\mu\nu}$ (with $p\in\mathbb{Z}^{+}$) where
the expectation value is expected to deviate from semiclassical dynamics. For
expectation values of the operator $\hat{X}_{\theta}^{n}$ given in Eq. (33),
there are a set of times $t^{(n)}_{M,N}$ where at least one of the terms in
the expansion of $\hat{X}_{\theta}^{n}$ will exhibit recurrence, leading to an
overall deviation from the TWA. These times are defined as the recurrence
times that satisfy $\kappa t^{(n)}_{M,N}=2\pi M/N$ where
$M,N\in\mathbb{Z}^{+}$ and $N$ is even(odd) for $n$ even(odd), and in both
cases $N\leq n$. The case with $N=0$ comes from terms with $\mu=\nu$, that is,
operators that are proportional to powers of the number operator, and these
never exhibit recurrence under the Kerr interaction444Note that the number of
recurrences between $t=0$ and $t=\frac{2\pi}{\kappa}$ for an operator
$\hat{X}_{\theta}^{n}$ can be compactly written in terms of the Euler totient
summary function, and thus scales quadratically:
$\sim\frac{1}{\xi(2)}n^{2}+\mathscr{O}(n{\rm log}n)$, where
$\xi(2)=\frac{2\pi}{6}$ is the Riemann zeta function evaluated at $2$. In this
way, the high spatial frequency features of the state characterized by higher-
order moments are connected to high temporal frequency dynamics, exhibiting
recurrence more frequently.. Note that the set of recurrence times
$t^{(n)}_{M,N}$ for any particular operator $\hat{X}_{\theta}^{n}$ is a strict
subset of the times $t_{M,N}$ when kitten states form; for the latter
$t_{M,N}$, there is no restriction on $M$ and $N$ other than they are positive
integers and co-prime (hence, forming a dense set), while for the former
$t^{(n)}_{M,N}$, we require $N\leq n$ and $N$ to be even (odd) if $n$ is even
(odd)555We do not require co-prime for the recurrence times as we did for the
Kitten times, as we do not need to avoid double counting. Indeed, when there
are multiple values of $N$ and $M$ corresponding to the same $t^{(n)}_{M,N}$,
we observe a larger deviation from the TWA at the recurrence time due to
contributions from multiple terms in the expansion in Eq. (33)..

Figure 5: Error induced by finite grid effects on the calculation of
$\langle{\hat{X}_{\theta}^{6}}\rangle(t)$ measured relative to $\Delta r_{\rm
min}=10^{-2}\frac{\pi}{\alpha_{0}^{2}}$, for $\alpha_{0}$ varying from $3$ to
$6$, and averaging over $20$ values of $\theta$. We define the error in Eq.
(36), and offset all data by machine precision $10^{-24}$. The maximum time
$t_{max}$ is chosen to be $1.3\frac{\pi}{3\kappa}$ so that, for an initial
coherent state at $t=0=t_{\rm min}$, the full behavior of the first recurrence
for $\left\langle\hat{X}_{\theta}^{6}\right\rangle(t)$ illustrated in Fig. 3
is captured. We comment that, for $\alpha_{0}=3$, the ratio of the error
induced by discretization for the closed system to the largest $\gamma$
simulation at $\Delta r=0.1$ is three orders of magnitude, whereas for
$\alpha_{0}=6$ the same ratio also at $\Delta r=0.1$ is four orders of
magnitude.

Consider now the deviation of the exact evolution from that given by
semiclassical TWA,

$\displaystyle\delta\langle\hat{X}^{n}_{\theta}\rangle(t)\equiv\langle\hat{X}^{n}_{\theta}\rangle(t)-\langle\hat{X}^{n}_{\theta}\rangle_{\rm
TWA}(t).$ (34)

Whether the deviation from the TWA is maximal precisely at $t^{(n)}_{M,N}$
depends on the particular phase-quadrature as well as the choice of
$\alpha_{0}$ (due to the fast $\alpha_{0}$-dependent oscillation in Eq. (II)).
We thus define an average deviation from the TWA for our expectation values of
interest

$\displaystyle\overline{\delta}^{(n)}(t)\equiv\frac{1}{|\alpha_{0}|^{n}}\int\frac{d\theta}{2\pi}\bigr{|}\delta\langle\hat{X}^{n}_{\theta}\rangle(t)|,$
(35)

which we normalize by $|\alpha_{0}|^{n}$ so we can meaningfully compare
expectation values with differing $\alpha_{0}$.

As shown in Fig. 4, at times far from any recurrences $t^{(n)}_{M,N}$ the
behavior of the expectation values are well approximated by the TWA (as
predicted from the last two cases of Eq. (IV)). This approximation becomes
more exact in the large-$\alpha_{0}$ limit, as seen with the increasing
sharpness of the recurrence with increasing $\alpha_{0}$. In this limit the
coherent states in the kitten-states become increasingly orthogonal, as
detailed above, which also serves to explain the decrease in initial
deviations from the TWA at short times as $\alpha_{0}$ increases, as has been
detailed elsewhere Polkovnikov (2010, 2003). In Fig. 4 we also observe that
the TWA improves as $\gamma$ increases, and for the same decoherence rate, the
TWA improves for larger $\alpha_{0}$ as expected in the macroscopic limit.

Given the average deviation of the exact solution from the TWA (Eq. (35)), we
now define our quantitative metric of interest. We numerically integrate
Eq.(24) to find the Wigner function and calculate expectation values based a
finite difference method, as described in detail in Appendix A. Given the
symmetry, we decompose the Wigner function into a Fourier series in the
angular variable and discretize the radial variable with resolution $\Delta
r$. The coarseness of this grid determines efficiency of the numerical
integration. Thus, we define a cumulative relative error induced by the phase
space discretization in the finite difference method,

$\displaystyle\epsilon^{(n)}_{\tau}(\Delta
r)\equiv\frac{\int_{\tau}dt\left|\overline{\delta}^{(n)}_{\Delta
r}(t)-\overline{\delta}^{(n)}_{dr}(t)\right|}{\int_{\tau}dt\left|\overline{\delta}^{(n)}_{dr}(t)\right|}.$
(36)

Here $dr$ is a suitably fine-grained radial grid size that well approximates
the exact solution, $\tau$ is a time-window of interest. This is plotted in
Fig. 5 for the moment $n=6$ with a variety of decoherence rates $\gamma$ and a
time window $\tau$ that includes the first recurrence
$t^{(6)}_{1,6}=\pi/(3\kappa)$. This metric has the advantage that, if
increasing the coarse-grained grid size $\Delta r$ has the effect of over-
estimating $\overline{\delta}^{(n)}_{\tau}$ at some times and under-estimating
$\overline{\delta}^{(n)}_{\tau}$ at other times, these effects will not cancel
out.

For the cumulative error, in all cases the error resulting from coarse-
graining the grid is suppressed for larger decoherence strengths. That this
becomes more prominent for larger $\alpha_{0}$ is made evident by carefully
comparing the top and bottom plots in Fig. 5: the separation between the lines
corresponding to the closed and open system increases with $\alpha_{0}$ at
intermediate $\Delta r$, especially for the smaller values of $\gamma$ where
the subPlanck-scale features are not yet completely washed out. Since these
fine-grained features are generated more prominently for larger $\alpha_{0}$,
we expect the error to increase with $\alpha_{0}$, which it indeed does for
the closed system. However, the error induced in calculations of the open
system increases more slowly with grid size $\Delta r$. In the limit of
asymptotically-strong decoherence, we would expect the error to be completely
flat at intermediate $\Delta r$ since no subPlanck features arise.666At
sufficiently high $\Delta r$ there will always be discretization errors due to
failure to capture even the semiclassical dynamics given by the TWA, including
the open system effect of energy loss.

##  V Summary and Outlook

In this work we have studied the effect of decoherence on the efficiency of
simulating quantum dynamics in the case of an anharmonic Kerr oscillator, a
standard paradigm in quantum optics. As the decoherence facilitates the
quantum-to-classical transition and washes out quantum complexity, we expect
more efficient representations are possible for the open quantum system. To
quantify this, we studied quantum phase space dynamics represented by the
Wigner function. Whereas the closed-system quantum dynamics leads to fine-
grained subPlanck scale structure, decoherence acts to wash this out, making
the semiclassical description according to the TWA more accurate. This
semiclassical dynamics is more efficiently simulated by a coarse-grained
finite-difference integration of the Wigner function, when including the
Fokker-Planck terms associated with decoherence.

For the specific cases of the Kerr oscillator we showed how semiclassical and
quantum dynamics diverge due to quantum interference and the generation of
superpositions of coherent states. Without the fine structure refeeding
coherence to the system, evolution of the open system does not produce these
high-order kitten states. The result of this is a reduction of the quantum
deviations from the TWA (Fig. 4) and a suppression of numerical errors induced
by simulation of the system on a coarser discrete grid for numerical
integration (Fig. 5).

The tendency towards expectation value dynamics given by the TWA is well-known
in the literature Polkovnikov (2003). Here, we have explained this in terms of
the breaking of the $\mathbb{Z}_{n}$ symmetry of the kitten states by
decoherence concurrent with Hamiltonian evolution. It remains an open question
whether the source of this reduction in computational cost is a reduction in
signal-to-noise ratio (such as is the case for depolarizing and dephasing
noise models for qubits), or a nontrivial transition from quantum dynamics to
classical dynamics. We conjecture that, for the anharmonic oscillator at zero
temperature, the open system is not a noisy simulation of the closed system
but a different (and easier) problem entirely.

Our work makes quantitative the intuition long established in the early work
on decoherence. Decoherence reduces the complexity of the quantum state in
phase space, and it does so more prominently for more macroscopic initial
states. This translates into more efficient simulation of the open quantum
system as the system size grows. While we have studied this for the toy
problem of a single anharmonic oscillator, its implications for obtaining a
quantum advantage in NISQ devices is another open question. How does the
rapidity of the regression to semiclassical dynamics with increasing system
size generalize for mulitpartite systems such as multiple coupled nonlinear
oscillators or qubits, and what is an efficient representation of the open
quantum systems that would allow for efficient classical simulations?

The current work gives some hints in this direction. Phase space
representations may provide for scalable efficient simulations of open quantum
systems when decoherence is sufficiently large. Coarse graining and well-
chosen finite discretization, informed by semiclassical dynamics is one
potential method, as seen in recent work Roda-Llordes  _et al._ (2023).
Another potential method is mapping the phase space dynamics to an underlying
set of stochastic Langevin equations. While such a method cannot efficiently
capture all complex quantum dynamics, it may do so in the presence of
decoherence Deuar  _et al._ (2021). In future work we will study this to
better understand when the quantum-to-classical transition is also a
transition in computational complexity.

##  VI Acknowledgements

TzBP gratefully acknowledges helpful conversations with Jun Takahashi,
Changhao Yi, and Chris Jackson.

The authors would like to thank the UNM Center for Advanced Research
Computing, supported in part by the National Science Foundation, for providing
the high performance computing resources used in this work. This work was
supported by National Science Foundation Grant No. PHY-2116246 and Grant No.
2037755, and is based upon work partially supported by the U.S. Department of
Energy, Office of Science, National Quantum Information Science Research
Centers, Quantum Systems Accelerator. This material is based upon work
supported by the Air Force Office of Scientific Research under award number
FA9550-22-1-0498. Any opinions, findings, and conclusions or recommendations
expressed in this material are those of the author(s) and do not necessarily
reflect the views of the United States Air Force.

We acknowledge the indigenous peoples of the Pueblo of Sandia as the original
inhabitants, stewards, and protectors of the lands on which the University of
New Mexico now sits.

## References

  * Zurek (1981) W. H. Zurek, “Pointer basis of quantum apparatus: Into what mixture does the wave packet collapse?” Phys. Rev. D 24, 1516–1525 (1981).
  * Zurek (1991) W. H. Zurek, “Decoherence and the transition from quantum to classical,” Physics Today 44, 36–44 (1991).
  * Zurek and Paz (1994) W. H. Zurek and J. P. Paz, “Decoherence, chaos, and the second law,” Physical Review Letters 72, 2508–2511 (1994).
  * Anglin  _et al._ (1997) J. R. Anglin, J. P. Paz, and W. H. Zurek, “Deconstructing decoherence,” Physical Review A 55, 4041–4053 (1997).
  * Zurek (1998) W. H. Zurek, “Decoherence, chaos, quantum-classical correspondence, and the algorithmic arrow of time,” Physica Scripta T76, 186 (1998).
  * Paz and Zurek (1999) J. P. Paz and W. H. Zurek, “Quantum limit of decoherence: Environment induced superselection of energy eigenstates,” Physical Review Letters 82, 5181–5185 (1999).
  * Karkuszewski  _et al._ (2002) Z. P. Karkuszewski, C. Jarzynski, and W. H. Zurek, “Quantum chaotic environments, the butterfly effect, and decoherence,” Physical Review Letters 89, 170405 (2002).
  * Zurek (2001) W. H. Zurek, “Sub-planck structure in phase space and its relevance for quantum decoherence,” Nature 412, 712–717 (2001).
  * Dodonov  _et al._ (1974) V. Dodonov, I. Malkin, and V. Man’ko, “Even and odd coherent states and excitations of a singular oscillator,” Physica 72, 597–615 (1974).
  * Leggett  _et al._ (1987) A. J. Leggett, S. Chakravarty, A. T. Dorsey, M. P. A. Fisher, A. Garg, and W. Zwerger, “Dynamics of the dissipative two-state system,” Reviews of Modern Physics 59, 1–85 (1987).
  * Caldeira and Leggett (1983) A. Caldeira and A. Leggett, “Quantum tunnelling in a dissipative system,” Annals of Physics 149, 374–456 (1983).
  * Kim and Bužek (1992) M. S. Kim and V. Bužek, “Schrödinger-cat states at finite temperature: Influence of a finite-temperature heat bath on quantum interferences,” Physical Review A 46, 4239–4251 (1992).
  * Shor (1995) P. W. Shor, “Scheme for reducing decoherence in quantum computer memory,” Phys. Rev. A 52, R2493–R2496 (1995).
  * Gottesman  _et al._ (2001) D. Gottesman, A. Kitaev, and J. Preskill, “Encoding a qubit in an oscillator,” Physical Review A 64 (2001), 10.1103/physreva.64.012310.
  * (15) V. V. Albert and P. Faist, “The error correction zoo,”  https://errorcorrectionzoo.org/.
  * Shor (1996) P. W. Shor, “Fault-tolerant quantum computation,” in  _Proceedings of 37th Conference on Foundations of Computer Science_ (IEEE, 1996) pp. 56–65.
  * Aharonov and Ben-Or (1997) D. Aharonov and M. Ben-Or, “Fault-tolerant quantum computation with constant error,” in  _Proceedings of the Twenty-Ninth Annual ACM Symposium on Theory of Computing_, STOC ’97 (Association for Computing Machinery, New York, NY, USA, 1997) pp. 176–188.
  * Kitaev (1997) A. Y. Kitaev, “Quantum computations: algorithms and error correction,” Uspekhi Matematicheskikh Nauk 52, 53–112 (1997).
  * Knill  _et al._ (1998) E. Knill, R. Laflamme, and W. H. Zurek, “Resilient quantum computation: error models and thresholds,” Proceedings of the Royal Society of London. Series A: Mathematical, Physical and Engineering Sciences 454, 365–384 (1998).
  * et al. (2023) R. A. et al., “Suppressing quantum errors by scaling a surface code logical qubit,” Nature 614, 676–681 (2023).
  * Preskill (2018) J. Preskill, “Quantum computing in the NISQ era and beyond,” Quantum 2, 79 (2018).
  * Curtright  _et al._ (2013) T. L. Curtright, D. B. Fairlie, and C. K. Zachos,  _A Concise Treatise on Quantum Mechanics in Phase Space_ (WORLD SCIENTIFIC, 2013).
  * Cahill and Glauber (1969) K. E. Cahill and R. J. Glauber, “Density operators and quasiprobability distributions,” Physical Review 177, 1882–1902 (1969).
  * Hudson (1974) R. Hudson, “When is the Wigner quasi-probability density non-negative?” Reports on Mathematical Physics 6, 249–252 (1974).
  * Kenfack and Życzkowski (2004) A. Kenfack and K. Życzkowski, “Negativity of the wigner function as an indicator of non-classicality,” Journal of Optics B: Quantum and Semiclassical Optics 6, 396–404 (2004).
  * Heller (1976) E. J. Heller, “Wigner phase space method: Analysis for semiclassical applications,” The Journal of Chemical Physics 65, 1289–1298 (1976).
  * Weedbrook  _et al._ (2012) C. Weedbrook, S. Pirandola, R. García-Patrón, N. J. Cerf, T. C. Ralph, J. H. Shapiro, and S. Lloyd, “Gaussian quantum information,” Reviews of Modern Physics 84, 621–669 (2012).
  * Sellier  _et al._ (2015) J. Sellier, M. Nedjalkov, and I. Dimov, “An introduction to applied quantum mechanics in the Wigner Monte Carlo formalism,” Physics Reports 577, 1–34 (2015).
  * Welland and Ferry (2020) I. Welland and D. K. Ferry, “Wavepacket phase-space quantum Monte Carlo method,” Journal of Computational Electronics 20, 267–273 (2020).
  * Habib  _et al._ (1998) S. Habib, K. Shizume, and W. H. Zurek, “Decoherence, chaos, and the correspondence principle,” Physical Review Letters 80, 4361–4365 (1998).
  * Milburn (1986) G. J. Milburn, “Quantum and classical Liouville dynamics of the anharmonic oscillator,” Phys. Rev. A 33, 674–685 (1986).
  * Milburn and Holmes (1986) G. J. Milburn and C. A. Holmes, “Dissipative quantum and classical Liouville mechanics of the anharmonic oscillator,” Physical Review Letters 56, 2237–2240 (1986).
  * Polkovnikov (2010) A. Polkovnikov, “Phase space representation of quantum dynamics,” Annals of Physics 325, 1790–1852 (2010).
  * Miranowicz  _et al._ (1990) A. Miranowicz, R. Tanaś, and S. Kielich, “Generation of discrete superpositions of coherent states in the anharmonic oscillator model,” Quantum Optics: Journal of the European Optical Society Part B 2, 253–265 (1990).
  * Tara  _et al._ (1993) K. Tara, G. S. Agarwal, and S. Chaturvedi, “Production of Schrödinger macroscopic quantum-superposition states in a Kerr medium,” Phys. Rev. A 47, 5024–5029 (1993).
  * van Enk (2003) S. J. van Enk, “Entanglement capabilities in infinite dimensions: Multidimensional entangled coherent states,” Phys. Rev. Lett. 91, 017902 (2003).
  * van Enk (2005) S. J. van Enk, “Decoherence of multidimensional entangled coherent states,” Phys. Rev. A 72, 022308 (2005).
  * Stobińska  _et al._ (2008) M. Stobińska, G. J. Milburn, and K. Wódkiewicz, “Wigner function evolution of quantum states in the presence of self-Kerr interaction,” Phys. Rev. A 78, 013810 (2008).
  * Yurke and Stoler (1986) B. Yurke and D. Stoler, “Generating quantum mechanical superpositions of macroscopically distinguishable states via amplitude dispersion,” Phys. Rev. Lett. 57, 13–16 (1986).
  * Stobińska  _et al._ (2011) M. Stobińska, A. S. Villar, and G. Leuchs, “Generation of Kerr non-Gaussian motional states of trapped ions,” EPL (Europhysics Letters) 94, 54002 (2011).
  * Greiner  _et al._ (2002) M. Greiner, O. Mandel, T. W. Hänsch, and I. Bloch, “Collapse and revival of the matter wave field of a Bose–Einstein condensate,” Nature 419, 51–54 (2002).
  * Kirchmair  _et al._ (2013) G. Kirchmair, B. Vlastakis, Z. Leghtas, S. E. Nigg, H. Paik, E. Ginossar, M. Mirrahimi, L. Frunzio, S. M. Girvin, and R. J. Schoelkopf, “Observation of quantum state collapse and revival due to the single-photon Kerr effect,” Nature 495, 205–209 (2013).
  * Grimm  _et al._ (2020) A. Grimm, N. E. Frattini, S. Puri, S. O. Mundhada, S. Touzard, M. Mirrahimi, S. M. Girvin, S. Shankar, and M. H. Devoret, “Stabilization and operation of a Kerr-cat qubit,” Nature 584, 205–209 (2020).
  * Sudheesh  _et al._ (2004) C. Sudheesh, S. Lakshmibala, and V. Balakrishnan, “Manifestations of wave packet revivals in the moments of observables,” Physics Letters A 329, 14–21 (2004).
  * Peinová and Luk (1990) V. Peinová and A. Luk, “Exact quantum statistics of a nonlinear dissipative oscillator evolving from an arbitrary state,” Phys. Rev. A 41, 414–420 (1990).
  * Chaturvedi and Srinivasan (1991a) S. Chaturvedi and V. Srinivasan, “Class of exactly solvable master equations describing coupled nonlinear oscillators,” Phys. Rev. A 43, 4054–4057 (1991a).
  * Chaturvedi and Srinivasan (1991b) S. Chaturvedi and V. Srinivasan, “Solution of the master equation for an attenuated or amplified nonlinear oscillator with an arbitrary initial condition,” J. Mod. Opt. 38, 777 (1991b).
  * McDonald and Clerk (2022) A. McDonald and A. A. Clerk, “Exact solutions of interacting dissipative systems via weak symmetries,” Phys. Rev. Lett. 128, 033602 (2022).
  * McDonald and Clerk (2023) A. McDonald and A. A. Clerk, “Third quantization of open quantum systems: new dissipative symmetries and connections to phase-space and Keldysh field theory formulations,” arXiv preprint arXiv:2302.14047 (2023), 10.48550/arXiv.2302.14047.
  * Walls and Milburn (2008) D. F. Walls and G. F. Milburn,  _Quantum Optics_ , 2nd ed. (Springer, Berlin, Germany, 2008).
  * Steel  _et al._ (1998) M. J. Steel, M. K. Olsen, L. I. Plimak, P. D. Drummond, S. M. Tan, M. J. Collett, D. F. Walls, and R. Graham, “Dynamical quantum noise in trapped Bose-Einstein condensates,” Phys. Rev. A 58, 4824–4835 (1998).
  * Propp (2023a) Tz. B. Propp, “Video of the Wigner function evolving in a Kerr medium generated using the Fock-space method described in Appendix B for the open and (approximately) closed system. Here, we depict the rotation counter-clockwise as is often done in the literature Stobińska  _et al._ (2011), as opposed to the (physical) clockwise rotation depicted elsewhere in the main text.” https://github.com/tzpropp/KerrVideos/blob/main/WignerMovieFull.gif (2023a).
  * Bialynicka-Birula (1968) Z. Bialynicka-Birula, “Properties of the generalized coherent state,” Phys. Rev. 173, 1207–1209 (1968).
  * Polkovnikov (2013) A. Polkovnikov, “Phase space representation of quantum dynamics. Lecture notes. Boulder summer school,”  (2013).
  * Nash (2016) D. Nash,  _A friendly introduction to group theory_ (2016).
  * Propp (2023b) Tz. B. Propp, “Video illustrating effect of the operator $\hat{a}^{6}$ on the phases in a 6-kitten superposition state.” https://github.com/tzpropp/KerrVideos/blob/main/PhaseGif.gif (2023b).
  * Polkovnikov (2003) A. Polkovnikov, “Quantum corrections to the dynamics of interacting bosons: Beyond the truncated Wigner approximation,” Physical Review A 68 (2003), 10.1103/PhysRevA.68.053604.
  * Roda-Llordes  _et al._ (2023) M. Roda-Llordes, D. Candoli, P. T. Grochowski, A. Riera-Campeny, T. Agrenius, J. J. García-Ripoll, C. Gonzalez-Ballestero, and O. Romero-Isart, “Numerical simulation of large-scale nonlinear open quantum mechanics,”  (2023).
  * Deuar  _et al._ (2021) P. Deuar, A. Ferrier, M. Matuszewski, G. Orso, and M. H. Szymańska, “Fully quantum scalable description of driven-dissipative lattice models,” PRX Quantum 2, 010319 (2021).
  * Bank  _et al._ (1985) R. Bank, W. Coughran, W. Fichtner, E. Grosse, D. Rose, and R. Smith, “Transient simulation of silicon devices and circuits,” IEEE Transactions on Computer-Aided Design of Integrated Circuits and Systems 4, 436–451 (1985).
  * Propp and van Enk (2019) Tz. B. Propp and S. J. van Enk, “Quantum networks for single photon detection,” Phys. Rev. A 100, 033836 (2019).

##  Appendix A Summary of Numerical Approach

Our aim is to numerically solve the (linear) partial differential equation
(PDE) for the Wigner function in a dissipative medium with a self-Kerr
interaction in cylindrical coordinates Stobińska  _et al._ (2008)

$\displaystyle\partial_{\tau}W(\tau,r,\varphi)$ $\displaystyle=$
$\displaystyle\left\\{\underbrace{(r^{2}-1)\partial_{\varphi}}_{\text{mean
field
rotation}}-\frac{1}{16}\left(\underbrace{\frac{1}{r}\partial_{r}}_{\text{Gaussian
shearing}}+\underbrace{\partial_{r}^{2}+\frac{1}{r^{2}}\partial_{\varphi}^{2}}_{\text{non-
Gaussian rotation}}\right)\partial_{\varphi}\right.$ (37)
$\displaystyle\left.+\underbrace{\xi+\frac{\xi}{2}r\partial_{r}}_{\text{drift}}+\underbrace{\frac{\xi}{8}\left(\partial_{r}^{2}+\frac{2}{r}\partial_{r}+\frac{1}{r^{2}}\partial_{\varphi}^{2}\right)}_{\text{diffusion}}\right\\}W(\tau,r,\varphi)\
,$

where we have defined dimensionless time $\tau=\kappa t$ and dimensionless
decay $\xi=\frac{\gamma}{\kappa}$ with boundary conditions

$W(\tau,0,0)=\frac{2}{\pi}e^{-2|\alpha_{0}|^{2}e^{-\tau\xi}}\
,\quad\lim_{r\to\infty}W(\tau,r,\varphi)=0\ ,$ (38)

and initial condition

$W(0,r,\varphi)=\frac{2}{\pi}e^{-2|\alpha_{0}-re^{i\varphi}|^{2}}\ .$ (39)

Without loss of generality assuming $\alpha_{0}$ to be real, we have

$\displaystyle W(0,r,\varphi)$ $\displaystyle=$
$\displaystyle\frac{2}{\pi}e^{-2(\alpha_{0}^{2}+r^{2})}e^{4\alpha
r\cos\varphi}$
$\displaystyle=\frac{2}{\pi}e^{-2(\alpha_{0}^{2}+r^{2})}\left[I_{0}(4r\alpha_{0})+2\sum_{k=1}^{\infty}I_{k}(4r\alpha)\cos(k\varphi)\right]\
,$

where $I_{k}(x)$ are the modified Bessel functions of the first kind. Because
the coefficients of Eq. (37) are independent of $\varphi$ and the Wigner
function must be periodic in $\varphi$, we choose to expand the Wigner
function in a sine/cosine series

$\displaystyle W(\tau,r,\varphi)$ $\displaystyle=$ $\displaystyle
a_{0}(\tau,r)+\sum_{k=1}^{\infty}\left(a_{k}(\tau,r)\cos(k\varphi)\right.$
(41) $\displaystyle\left.+b_{k}(\tau,r)\sin(k\varphi)\right)\ .$

The normalization condition can now be expressed as

$\int_{0}^{2\pi}d\varphi\int_{0}^{\infty}dr\ rW(\tau,r,\varphi)=2\pi\int dr\
ra_{0}(\tau,r)=1\ .$ (42)

Plugging the expansion of the Wigner function into Eq. (37), we get a coupled
system of ordinary differential equations (ODE’s) to solve for each $k$ value

$\displaystyle\partial_{\tau}a_{0}(\tau,r)$ $\displaystyle=$
$\displaystyle\left\\{\xi+\frac{\xi}{2}\left(r+\frac{1}{4r}\right)\partial_{r}+\frac{\xi}{8}\partial_{r}^{2}\right\\}a_{0}(\tau,r)\
,$ (43) $\displaystyle\partial_{\tau}a_{k}(\tau,r)$ $\displaystyle=$
$\displaystyle\left\\{k(r^{2}-1)-\frac{k}{16}\left(\frac{1}{r}\partial_{r}+\partial_{r}^{2}-\frac{k^{2}}{r^{2}}\right)\right\\}b_{k}(\tau,r)$
(44)
$\displaystyle+\left\\{\xi+\frac{\xi}{2}\left(r+\frac{1}{4r}\right)\partial_{r}+\frac{\xi}{8}\left(\partial_{r}^{2}-\frac{k^{2}}{r^{2}}\right)\right\\}a_{k}(\tau,r)\
,$ $\displaystyle\partial_{\tau}b_{k}(\tau,r)$ $\displaystyle=$
$\displaystyle\left\\{-k(r^{2}-1)+\frac{k}{16}\left(\frac{1}{r}\partial_{r}+\partial_{r}^{2}-\frac{k^{2}}{r^{2}}\right)\right\\}a_{k}(\tau,r)$
(45)
$\displaystyle+\left\\{\xi+\frac{\xi}{2}\left(r+\frac{1}{4r}\right)\partial_{r}+\frac{\xi}{8}\left(\partial_{r}^{2}-\frac{k^{2}}{r^{2}}\right)\right\\}b_{k}(\tau,r)\
.$

The equations for different $k$’s do not couple, but for each $k$ the pair
$a_{k}(\tau,r)$ and $b_{k}(\tau,r)$ are governed by coupled ODE’s. The
boundary condition at $r=0$ is given by

$a_{0}(\tau,0)=\frac{2}{\pi}e^{-2|\alpha_{0}|^{2}e^{-\tau\xi}}\ ,\quad
a_{k}(\tau,0)=0\ ,\quad b_{k}(\tau,0)=0\ ,$ (46)

and the boundary condition at $r\to\infty$ is given by

$\displaystyle\lim_{r\to\infty}a_{0}(\tau,r)$ $\displaystyle=0\ ,$ (47a)
$\displaystyle\lim_{r\to\infty}a_{k}(\tau,r)$ $\displaystyle=0\ ,\quad k>0\ ,$
(47b) $\displaystyle\lim_{r\to\infty}b_{k}(\tau,r)$ $\displaystyle=0\ ,\quad
k>0\ ,$ (47c)

The initial condition is given by (again, assuming $\alpha_{0}$ is real)

$\displaystyle a_{0}(0,r)$
$\displaystyle=\frac{1}{\pi^{2}}e^{-2(\alpha_{0}^{2}+r^{2})}\int_{-\pi}^{\pi}d\varphi\
e^{4r\alpha_{0}\cos(\varphi)}$
$\displaystyle=\frac{2}{\pi}e^{-2(\alpha_{0}^{2}+r^{2})}I_{0}(4r\alpha_{0})\
,$ (48a) $\displaystyle a_{k}(0,r)$
$\displaystyle=\frac{2}{\pi^{2}}e^{-2(\alpha_{0}^{2}+r^{2})}\int_{-\pi}^{\pi}d\varphi\
e^{4r\alpha_{0}\cos(\varphi)}\cos(k\varphi)$
$\displaystyle=\frac{4}{\pi}e^{-2(\alpha_{0}^{2}+r^{2})}I_{k}(4r\alpha_{0})\
,\quad k>0\ ,$ (48b) $\displaystyle b_{k}(0,r)$ $\displaystyle=0\ ,\quad k>0\
.$ (48c)

From the initial condition, we can find the leading behavior for the functions
near $r=0$

$\displaystyle a_{0}(0,r)$
$\displaystyle=\tilde{a}_{0}+r^{2}\tilde{\tilde{a}}_{0}+\dots\ ,$ (49a)
$\displaystyle a_{k}(0,r)$
$\displaystyle=\tilde{a}_{k}r^{k}+\tilde{\tilde{a}}_{k}r^{k+2}+\dots\ .$ (49b)

We find that this provides a consistent expansion for the functions at $r=0$
even for $\tau>0$

$\displaystyle a_{0}(\tau,r)$
$\displaystyle=\tilde{a}_{0}(\tau)+r^{2}\tilde{\tilde{a}}_{0}(\tau)+\dots\ ,$
(50a) $\displaystyle a_{k}(\tau,r)$
$\displaystyle=\tilde{a}_{k}(\tau)r^{k}+\tilde{\tilde{a}}_{k}r^{k+2}(\tau)+\dots\
,$ (50b) $\displaystyle b_{k}(\tau,r)$
$\displaystyle=\tilde{b}_{k}(\tau)r^{k}+\tilde{\tilde{b}}_{k}r^{k+2}(\tau)+\dots\
.$ (50c)

The ODE’s are solved by discretizing the radial direction into $n$ equal steps

$r_{j}=j\Delta r\ ,\quad j=1,\dots n\ ,$ (51)

with $\Delta r=r_{\mathrm{max}}/n$, where $r_{\mathrm{max}}$ being some
sufficiently large $r$ value. In our simulations, we take
$r_{\mathrm{max}}=2\alpha_{0}$. We will use the index $j=0$ to denote the
boundary condition at $r=0$. The derivatives discretized up to fourth order
are given by

$\displaystyle\partial_{r}f(\tau,r)$ $\displaystyle\rightarrow$
$\displaystyle\frac{f_{j-2}(\tau)-8f_{j-1}(\tau)+8f_{j+1}(\tau)-f_{j+2}(\tau)}{12\Delta
r}\ ,$ (52) $\displaystyle\partial_{r}^{2}f(\tau,r)$
$\displaystyle\rightarrow$
$\displaystyle\frac{-f_{j-2}(\tau)+16f_{j-1}(\tau)-30f_{j}(\tau)+16f_{j+1}(\tau)-f_{j+2}(\tau)}{12\Delta
r^{2}}\ ,$ (53)

where $f_{j}(\tau)\equiv f(\tau,r_{j})$. Let us define the vector
$\vec{g}_{k}(\tau)=\left(a_{k,1}(\tau),\dots,a_{k,n}(\tau),b_{k,1}(\tau),\dots
b_{k,n}(\tau)\right)^{T}$. We will then write our ODE’s as a system of linear
equations to be solved

$\displaystyle\partial_{t}\vec{g}_{0}(\tau)$
$\displaystyle=A_{0}\vec{g}_{0}(\tau)+\vec{b}_{0}(\tau)\ ,$ (54a)
$\displaystyle\partial_{t}\vec{g}_{k}(\tau)$
$\displaystyle=A_{k}\vec{g}_{k}(\tau)\ ,k>0\ .$ (54b)

The vector $\vec{b}_{0}(\tau)$ arises from the discretization near $r=0$: when
terms such as $a_{0,0}(\tau)$ appear when evaluating the ODE at $j=1$ and
$j=2$, we must replace them with the initial condition

$a_{0,0}(\tau)=\frac{2}{\pi}e^{-2\alpha^{2}}\ ,$ (55)

which then represent constants in the system of equations. (Recall that
$a_{k,0}(\tau)=0$ so this term does not appear for $k>0$.) The terms
$a_{k,-1}(\tau)$ and $b_{k,-1}(\tau)$ that arise at $j=1$ must be replaced by

$a_{k,-1}(\tau)=(-1)^{k}a_{k,1}(\tau)\ ,\quad
b_{k,-1}(\tau)=(-1)^{k}b_{k,1}(\tau)\ .$ (56)

This choice is to enforce the even/odd behavior from Eq. (50). We note that
this term only appears from the term $\frac{\xi}{2}r\partial_{r}$, and it is
cancelled otherwise. We enforce the boundary condition at infinity by setting
to zero terms beyond $r_{\mathrm{max}}$, i.e.
$a_{k,n+1}(\tau)=a_{k,n+2}(\tau)=b_{k,n+1}(\tau)=b_{k,n+2}(\tau)=0$.

To solve the system of linear equations in Eq. (54), we use TR-BDF2 Bank  _et
al._ (1985), an implicit Runge-Kutta method. We use an absolute error
tolerance of $10^{-9}$ and a relative error tolerance of $10^{-6}$. Our
maximum temporal step size is chosen to be $\pi/(100\alpha_{0}^{2})$. We
simulate pairs of $(a_{k},b_{k})$ up to a maximum $k$ value of 60.

After solving for the time-evolved Wigner function up to some time $\tau$, we
can plot the Wigner function in terms of $(x,p)$ for fixed $\kappa t$ as was
done in the main text by plotting the function
$a_{0}(\tau,r)+\sum_{k=1}^{\infty}\left(a_{k}(\tau,r)\cos(k\varphi)+b_{k}(\tau,r)\sin(k\varphi)\right)$.
Our $\varphi$ grid is discretized in steps of $\pi/200$.

We can also calculate the expectation value of an observable $\hat{A}$ as (we
drop the time dependence for now)

$\langle\hat{A}\rangle=\int_{-\infty}^{\infty}dx\int_{-\infty}^{\infty}dp\
W(x,p)\tilde{A}(x,p)\ ,$ (57)

where $\tilde{A}(x,p)$ is the Weyl transform of the operator $\tilde{A}$

$\displaystyle\tilde{A}(x,p)$ $\displaystyle=$ $\displaystyle\int
dye^{-ipy/\hbar}\langle x+\frac{y}{2}|\tilde{A}|x-\frac{y}{2}\rangle$ (58)
$\displaystyle=$ $\displaystyle\int due^{ixu/\hbar}\langle
p+\frac{u}{2}|\tilde{A}|p-\frac{u}{2}\rangle\ .$

In particular, the expectation value of arbitrary powers of $\hat{X}$ and
$\hat{P}$ is given by

$\langle\hat{X}^{m}\hat{P}^{n}\rangle=\int_{-\infty}^{\infty}dx\int_{-\infty}^{\infty}dp\
x^{m}p^{n}W(x,p)\ ,$ (59)

Using our solution for the Wigner function in polar coordinates, this takes
the form

$\displaystyle\langle\hat{X}^{m}\hat{P}^{n}\rangle$ $\displaystyle=$
$\displaystyle\sum_{k=0}^{\infty}\left[\int_{0}^{2\pi}d\varphi\
\cos^{m}\varphi\sin^{n}\varphi\cos(k\varphi)\int_{0}^{\infty}dr\
r^{m+n+1}a_{k}(r)\right.$ (60) $\displaystyle\left.+\int_{0}^{2\pi}d\varphi\
\cos^{m}\varphi\sin^{n}\varphi\sin(k\varphi)\int_{0}^{\infty}dr\
r^{m+n+1}b_{k}(r)\right]\ .$

The integral over $\varphi$ can be calculated (numerically) by noting that

$\displaystyle\int_{0}^{2\pi}d\varphi\
\cos^{m}\varphi\sin^{n}\varphi\cos(k\varphi)$
$\displaystyle=\frac{2\pi}{2^{m+n+1}i^{n}}\sum_{k_{m}=0}^{m}\sum_{k_{n}=0}^{n}\binom{m}{k_{m}}\binom{n}{k_{n}}(-1)^{k_{n}}\left(\delta_{m+n+k-2k_{m}-2k_{n},0}+\delta_{m+n-k-2k_{m}-2k_{n},0}\right)\
,$ (61a) $\displaystyle\int_{0}^{2\pi}d\varphi\
\cos^{m}\varphi\sin^{n}\varphi\sin(k\varphi)$
$\displaystyle=\frac{2\pi}{2^{m+n+1}i^{n+1}}\sum_{k_{m}=0}^{m}\sum_{k_{n}=0}^{n}\binom{m}{k_{m}}\binom{n}{k_{n}}(-1)^{k_{n}}\left(\delta_{m+n+k-2k_{m}-2k_{n},0}-\delta_{m+n-k-2k_{m}-2k_{n},0}\right)\
.$ (61b)

In this form, it is clear that we need only to perform radial integrals of the
form $\int_{0}^{\infty}dr\ r^{p+1}a_{k}(r),b_{k}(r)$.

##  Appendix B Derivation of Recurrence Time Formulae and Zero-Temperature
Fock-Space Solution for the Open Quantum System

In Eq. (28), we observe that the density matrix for a bosonic mode with an
anharmonic, number-preserving Hamiltonian has the formal solution

$\displaystyle\rho_{jk}(t)=\sum_{j^{\prime},k^{\prime}}A_{j^{\prime}k^{\prime}}^{(j,k)}e^{-\tilde{\gamma}_{j^{\prime}k^{\prime}}t}\delta_{j^{\prime}-k^{\prime},j-k},$
(62)

where $\tilde{\gamma}_{j^{\prime}k^{\prime}}$ is a complex coefficient
specific to each density matrix element, depending on both the Hamiltonian
operator and Lindbladian superoperator’s elements in the Fock basis. Such an
equation is a solution to the family of equations

$\displaystyle\frac{d}{dt}\rho_{jk}(t)=-\tilde{\gamma}_{jk}\rho_{jk}(t)+\sum_{n}\left(g^{(n)\,+}_{jk}\rho_{j+n,k+n}(t)+g^{(n)\,-}_{jk}\rho_{j-n,k-n}(t)\right),$
(63)

which are generated by the Lindblad master equation with multiple arbitrary
$n$-photon-loss/gain channels, described by a coupling $g^{(n)\,\pm}_{jk}$
respectively. Recast in this way, Eq. (63) has the form of a set of coupled
systems where population can leak out either to other systems or to an
environment.777Explicitly, we work with a time domain generalization of Eq.
(28) in Propp and van Enk (2019), making a reindexing substitution
$\rho_{jk}\rightarrow c_{i}$, which recasts the evolution of a single stripe
of the density matrix, as depicted in Fig. 6, with respect to an anharmonic
Hamiltonian and open system effects as energy moving between simple harmonic
oscillators which are coupled to each other and to an unmonitored environment.
Using techniques for such systems Propp and van Enk (2019), we can solve each
diagonal stripe of the density matrix separately, as depicted in Fig. 6 for a
small number of basis elements.

Figure 6: For open quantum systems with anharmonic, number-preserving
Hamiltonians, each diagonal stripe of the density matrix in the number basis
(colored) evolves independently.

If we specialize to the Kerr Hamiltonian defined in Eq. (1) and zero-
temperature with amplitude damping (that is, single photon loss only) giving
rise to the Lindbladian in Eq. (19), the complex decay rates
$\tilde{\gamma}_{jk}$ have the form

$\displaystyle\tilde{\gamma}_{jk}=\frac{\gamma}{2}(j+k)+i\frac{\kappa}{2}\left((j^{2}-j)-(k^{2}-k)\right),$
(64)

and the couplings are given by

$\displaystyle g^{(1)\,+}_{jk}=\gamma\sqrt{(j+1)(k+1)},$ (65)

with all other $g^{(n)\,\pm}_{jk}=0$.

Returning to the more general case of an anharmonic number-preserving
oscillator with arbitrary $n$-photon loss and gain channels, we will now show
that the two results derived in the main text hold under this more general
case: recurrence times are independent of both operator ordering and open
system effects. Here we define $s=+1$ to be normally ordered, $s=0$ to be
symmetrically ordered, and $s=-1$ to be anti-normally ordered.

We can decompose the expectation value of an $s$-ordered product of powers of
the creation and annihilation operators in terms of a single stripe of the
density matrix

$\displaystyle\left\langle\hat{a}^{\dagger\mu}\hat{a}^{\nu}\right\rangle_{s}(t)=\sum\limits_{jk}P_{jk}^{\mu\nu
s}\rho_{jk}(t)\delta_{j-k,\mu-\nu},$ (66)

with $P_{jk}^{(\mu\nu s)}$ a positive and real combinatorial factor.
Substituting in the solution for density matrix elements given in Eq. (62), we
find

$\displaystyle\left\langle\hat{a}^{\dagger\mu}\hat{a}^{\nu}\right\rangle_{s}(t)=\sum_{j^{\prime}k^{\prime}}R_{j^{\prime}k^{\prime}}^{\mu\nu
s}e^{-\tilde{\gamma}_{j^{\prime}k^{\prime}}t}\delta_{j^{\prime}-k^{\prime},\mu-\nu},$
(67)

where we have defined new Fourier coefficients for the expectation value

$\displaystyle R_{j^{\prime}k^{\prime}}^{\mu\nu
s}=\sum_{j}P_{j^{\prime}-k^{\prime}+j,j}^{\mu\nu
s}A_{j^{\prime}k^{\prime}}^{j^{\prime}-k^{\prime}+j,j}.$ (68)

where we reindexed $j\rightarrow j^{\prime}-k^{\prime}+j$ and contracted the
Kronecker product.

Note that it is the sums defining the Fourier coefficients in Eq. (68) that
collapse to a single term for the closed system, and not the sum in the
expectation value in Eq. (67). Thus, there is a set of frequencies
$\omega_{j}$ present in the time-dependent expectation value, where
$\omega_{j}={\rm Im}[\tilde{\gamma}_{j,j+\mu-\nu}]$ for $\mu>\nu$ and
$\omega_{j}={\rm Im}[\tilde{\gamma}_{j+\nu-\mu,j}]$ for $\mu<\nu$ (the two
cases ensure that the indices on $\tilde{\gamma}_{j,k}$ remain nonnegative).
For any anharmonic number-preserving Hamiltonian, all these frequencies are
distinct.

To see this, consider a Hamiltonian decomposed into the form

$\displaystyle H_{\rm anh}=\sum_{k}b_{k}\hat{n}^{k}.$ (69)

We see that, for an expectation value
$\left\langle\hat{a}^{\dagger\mu}\hat{a}^{\nu}\right\rangle_{s}(t)$ (assuming
$\mu<\nu$) the frequencies for the closed system will have the form

$\displaystyle\omega_{j}=\sum_{k}b_{k}(j^{k}-(j+\nu-\mu)^{k}),$ (70)

which are negative in this case. For the harmonic case, $b_{1}$ is the only
non-zero term, and the frequencies $\omega_{j}$ within a stripe are
independent of $j$ and thus identical (or when $\mu=\nu$ so that there simply
is no frequency). In all other cases, each frequency depends on $j$ and thus
the sum in Eq. (67) will not result in cancellations of any particular
frequency, no matter what the Fourier coefficients
$R_{j^{\prime}k^{\prime}}^{\nu\mu s}$ may be, provided they are nonzero.

That they are nonzero for the closed system is trivial; the sum
$R_{j^{\prime}k^{\prime}}^{\nu\mu s}$ reduces to a single term corresponding
to the initial conditions for a density matrix element
$\rho_{j^{\prime}k^{\prime}}(0)$ multiplied by a positive combinatorial
factor.

To show that the Fourier coefficients are non-zero for the open system is the
final ingredient to complete our more rigorous argument. We begin by noting
that the effect of arbitrary $n$-photon gain and loss channels is to modify
the real part of $\tilde{\gamma}_{jk}$ as well as the coefficients
$A_{j^{\prime}k^{\prime}}^{j^{\prime}-k^{\prime}+j,j}$.

This latter point is seen by inserting the formal solution for the density
matrix elements Eq. (62) into the coupled equations Eq. (63),

$\displaystyle\frac{d}{dt}\rho_{jk}(t)$ $\displaystyle=$
$\displaystyle-\sum_{j^{\prime}k^{\prime}}\tilde{\gamma}_{j^{\prime}k^{\prime}}A_{j^{\prime}k^{\prime}}^{(jk)}e^{-\tilde{\gamma}_{j^{\prime}k^{\prime}}t}\delta_{j^{\prime}-k^{\prime},j-k}$
$\displaystyle=$
$\displaystyle-\tilde{\gamma}_{jk}\sum_{j^{\prime}k^{\prime}}A_{j^{\prime}k^{\prime}}^{(jk)}e^{-\tilde{\gamma}_{j^{\prime}k^{\prime}}t}\delta_{j^{\prime}-k^{\prime},j-k}+\sum_{n}\left(g^{(n)\,+}_{jk}A_{j^{\prime}k^{\prime}}^{j+n,k+n}+g^{(n)\,-}_{jk}A_{j^{\prime}k^{\prime}}^{j-n,k-n}\right)e^{-\tilde{\gamma}_{j^{\prime}k^{\prime}}t}\delta_{j^{\prime}-k^{\prime},j-k},$

where the first line comes from performing the time derivative and the second
line comes from direct substituion into the right hand side of Eq. (63).
Making use of the linear independence of the functions
$e^{-\tilde{\gamma}_{j^{\prime}k^{\prime}}t}$ to establish equality for each
term, we rearrange to yield an expression

$\begin{split}A_{j^{\prime}k^{\prime}}^{(jk)}&(\tilde{\gamma}_{j^{\prime}k^{\prime}}-\tilde{\gamma}_{jk})+\\\
&\sum_{n}\left(g^{(n)\,+}_{jk}A_{j^{\prime}k^{\prime}}^{j+n,k+n}+g^{(n)\,-}_{jk}A_{j^{\prime}k^{\prime}}^{j-n,k-n}\right)=0.\end{split}$
(72)

Since the $g^{(n)\,\pm}_{jk}$ are non-zero positive numbers, we immediately
see that including an additional $n$-photon loss or gain channel must modify
the coefficient $A_{j^{\prime}k^{\prime}}^{(jk)}$. A corollary of this is that
each coefficient $A_{j^{\prime}k^{\prime}}^{(jk)}$ in a stripe is a function
of _every_ $\tilde{\gamma}_{jk}$ and $g^{(n)\,\pm}_{jk}$ within that same
stripe. Now, consider if a single $A_{j^{\prime}k^{\prime}}^{(jk)}$ were zero.
This implies that
$\sum\limits_{n}g^{(n)\,+}_{jk}A_{j^{\prime}k^{\prime}}^{j+n,k+n}+g^{(n)\,-}_{jk}A_{j^{\prime}k^{\prime}}^{j-n,k-n}=0$.
This requires an incredible amount of fine-tuning, since each
$A_{j^{\prime}k^{\prime}}^{j\pm n,k\pm n}$ depends on all system parameters.
Even if this holds, for the diagonal elements of the density matrix, the
coefficients $A_{j^{\prime}k^{\prime}}^{j\pm n,k\pm n}$ are necessarily real
and positive and so this condition can never be met. Thus, every coefficient
$A_{j^{\prime}k^{\prime}}^{(jk)}$ in the sum in Eq. (62) is non-zero.

Having shown that the Fourier coefficients $R_{j^{\prime}k^{\prime}}^{\nu\mu
s}$ of the expectation value are generally non-zero, we see that any
expectation value of the form
$\left\langle\hat{a}^{\dagger\mu}\hat{a}^{\nu}\right\rangle_{s}(t)$ inherits a
set of frequencies $\omega_{j}$ present in the $(\mu-\nu)$th stripe of the
density matrix, and that this set of frequencies is unchanged by open system
effects; while each frequency may contribute differently, it is still present
in the signal for all finite time. This set is also unchanged by operator
ordering, as it must be, since different orderings correspond to the same
operator physically. Mathematically, this is because the commutation relation
$[\hat{a},\hat{a}^{\dagger}]=1$ preserves the difference in powers of the
creation and annihilation operators; powers of both are removed together.
Because all frequencies that are present in a closed system expectation value
are also present in the open system expectation values, any periodicity of the
closed system expectation values must also be inherited by the open system
expectation values (The effect of the real part of $\tilde{\gamma}_{jk}$
coming from open system effects is to rescale the Fourier coefficients). Thus,
we have re-derived our results: recurrences of the closed system are inherited
by the open system irrespective of operator ordering.

We now turn to our second task, deriving analytic solutions to the Kerr
system’s evolution at zero temperature with single-photon loss utilizing this
Fock-space representation in terms of diagonal stripes of the density matrix.
We begin by truncating the Fock space at some
$N\gg\left\langle\hat{n}\right\rangle$, such that for each terminating row and
column of the density matrix

$\displaystyle\frac{d}{dt}\rho_{Nk}(t)$ $\displaystyle=$
$\displaystyle-\tilde{\gamma}_{Nk}\rho_{Nk}(t)\quad\forall k,$
$\displaystyle\frac{d}{dt}\rho_{jN}(t)$ $\displaystyle=$
$\displaystyle-\tilde{\gamma}_{jN}\rho_{jN}(t)\quad\forall j.$ (73)

This allows the final density matrix elements of each diagonal stripe to be
determined exactly

$\displaystyle\rho_{Nk}(t)$ $\displaystyle=$
$\displaystyle\rho_{Nk}(0)\,e^{-\tilde{\gamma}_{Nk}t}\quad\forall k$
$\displaystyle\rho_{jN}(t)$ $\displaystyle=$
$\displaystyle\rho_{jN}(0)\,e^{-\tilde{\gamma}_{jN}t}\quad\forall j.$ (74)

For simplicity, we focus on the bottom row of the density matrix, whose
elements are given exactly by the second line of Eq. (B). From these, we will
be able to recursively generate solutions to other density matrix elements
within the same stripe, and generate the rest of the density matrix by complex
conjugation.

To do this, we substitute our ansatz solution in terms of the generalized
discrete Fourier transform Eq. (62) into the simplified master equation for
the Kerr system at zero temperature

$\displaystyle\frac{d}{dt}\rho_{jk}(t)=-\tilde{\gamma}_{jk}\rho_{jk}(t)+g^{(1)\,+}_{jk}\rho_{j+1,k+1}(t),$
(75)

yielding

$\displaystyle\sum_{j^{\prime},k^{\prime}}\left(\left(\tilde{\gamma}_{j^{\prime}k^{\prime}}-\tilde{\gamma}_{jk}\right)A_{j^{\prime}k^{\prime}}^{(jk)}-g^{(1)\,+}_{jk}A_{j^{\prime}k^{\prime}}^{(j+1,k+1)}\right)e^{-\tilde{\gamma}_{j^{\prime}k^{\prime}}t}\delta_{j^{\prime}-k^{\prime}-(j-k)}=0.$
(76)

Here we proceed using the same methodology as we did in Eqs. (B) and (72):
rearranging and isolating each term in the sum.

Since the frequencies within any individual stripe are distinct, each term in
the sum is linearly independent and we conclude

$\displaystyle
A_{j^{\prime}k^{\prime}}^{(j,k)}=\frac{g^{(1)\,+}_{jk}}{\tilde{\gamma}_{j^{\prime}k^{\prime}}-\tilde{\gamma}_{jk}}A_{j^{\prime}k^{\prime}}^{(j+1,k+1)}.$
(77)

We now have a way to generate the Fourier coefficients in terms of the next
coefficient in a stripe. From Eq. (B), we see that there is only a single
nonzero term in the Fourier expansion of the bottom row of the density matrix.
The coefficients are

$\displaystyle
A_{j^{\prime}k^{\prime}}^{(N,k)}=\begin{cases}\rho_{Nk}(0),&j^{\prime}=N,k^{\prime}=k\\\
0,&{\rm otherwise}\end{cases}.$ (78)

Furthermore, because of the unidirectionality of the coupling in Eq. (75) for
zero temperature (that is, there is no way for population to travel from lower
Fock states to higher ones), we find

$\displaystyle A_{j^{\prime}k^{\prime}}^{(j,k)}=0,$ $\displaystyle
j>j^{\prime},\,k>k^{\prime}.$ (79)

We now shift notation slightly and explicitly limit our attention to just the
bottom left triangle of the truncated matrix: let $j\rightarrow N-i$ and
$k\rightarrow N-m-i$, where $N$ is the Fock space truncation defined
previously, and $m$ is how far off diagonal the stripe is, ranging from $0$
(the diagonal) to $N$ (the bottom left stripe consisting of a single element
of the truncated density matrix), and $i$ is the position within the stripe
from its bottom-rightmost terminating element, ranging from $0$ to $N-m$. In
this notation and making use of the truncated Fock space, we can now write an
analytic solution for the density matrix of the open quantum system.

We begin by rewriting, in this new notation, our formal solution Eq. (62)

$\displaystyle\rho_{N-i,N-m-i}(t)=\sum\limits_{i^{\prime}=0}^{i}A_{N-i^{\prime},N-m-i^{\prime}}^{(N-i,N-m-i)}e^{-\tilde{\gamma}_{N-i^{\prime},N-m-i^{\prime}}t},$

where we have made use of the result in Eq. (79) that
$A_{N-i^{\prime},N-m-i^{\prime}}^{(N,N-m)}=0$ for $i^{\prime}>i$ to truncate
the sum at $i$ instead of $N-m$. We also rewrite in the new notation the
recursion relation Eq. (77)

$\displaystyle\\!\\!A_{N-i^{\prime},N-m-i^{\prime}}^{(N-i,N-m-i)}=\frac{g^{(1)\,+}_{N-i,N-m-i}}{\tilde{\gamma}_{N-i^{\prime},N-m-i^{\prime}}-\tilde{\gamma}_{N-i,N-m-i}}A_{N-i^{\prime},N-m-i^{\prime}}^{(N-i+1,N-m-i+1)}$
$\displaystyle=\left(\prod\limits_{u=1}^{i-1}\frac{g^{(1)\,+}_{N-i+u,N-m-i+u}}{\tilde{\gamma}_{N-i^{\prime},N-m-i^{\prime}}-\tilde{\gamma}_{N-i+u,N-m-i+u}}\right)A_{N-i^{\prime},N-m-i^{\prime}}^{(N,N-m)}\,\,,$
(81)

where in the second line we have applied the recursion relation $i-1$ times
and have assumed $i>0$ and $i^{\prime}\leq i$ (again, because of the result in
Eq. (79)).

Defining the function

$\displaystyle
f_{Nmii^{\prime}}=\begin{cases}\prod\limits_{u=1}^{i-1}\frac{g^{(1)\,+}_{N-i+u,N-m-i+u}}{\tilde{\gamma}_{N-i^{\prime},N-m-i^{\prime}}-\tilde{\gamma}_{N-i+u,N-m-i+u}}&i>0,\\\
1&i=0,\\\ \end{cases}$

and substituting Eq. (B) into our formal solution Eq. (B), we find

$\begin{split}&\rho_{N-i,N-m-i}(t)\\\
&\qquad=\sum\limits_{i^{\prime}=0}^{i}f_{Nmii^{\prime}}A_{N-i^{\prime},N-m-i^{\prime}}^{(N,N-m)}e^{-\tilde{\gamma}_{N-i^{\prime},N-m-i^{\prime}}t}.\end{split}$
(83)

We have reduced the problem of the dynamics of a single stripe of the density
matrix elements $\rho_{N-i,N-m-i}(t)$ to finding the matrix of coefficients
$A_{N-i^{\prime},N-m-i^{\prime}}^{(N,N-m)}$, which we can now derive from
boundary terms, and evaluating Eq. (83). Notably, now the coefficients
$A_{N-i^{\prime},N-m-i^{\prime}}^{(N,N-m)}$ are independent of $i$ (the
location within the stripe) and can form a matrix, which is spanned by the
coordinates $i^{\prime}$ and $m$. One can then recursively solve for the
elements $A_{N-i^{\prime},N-m-i^{\prime}}^{(N,N-m)}$ for each $i$ at $t=0$ and
making use of the boundary condition of the known initial state

$\displaystyle\rho_{N-i,N-m-i}(0)=\sum\limits_{i^{\prime}=0}^{i}f_{Nmii^{\prime}}A_{N-i^{\prime},N-m-i^{\prime}}^{(N,N-m)}.$
(84)

For $i=0$ (the terminating row of the density matrix) we already have the
solution for density matrix element evolution in Eq. (B) and have written the
values of the coefficients in Eq. (78). Since there is a single term only, we
conclude for the $i=0$ case $A_{N,N-m^{\prime}}^{(N,N-m)}=\rho_{N,N-m}(0)$. By
then considering the $i=1$ case and using the boundary condition (Eq. (84)),
we can recursively generate the coefficients for larger and larger $i$ up to
its maximum value of $N-m$. For $i=1$, we find

$\displaystyle\rho_{N-1,N-m-1}(0)$ $\displaystyle=$ $\displaystyle
f_{Nm10}A_{N,N-m}^{(N,N-m)}$ $\displaystyle+f_{Nm11}A_{N-1,N-m-1}^{(N,N-m)}$
$\displaystyle\rightarrow A_{N-1,N-m-1}^{(N,N-m)}$ $\displaystyle=$
$\displaystyle\frac{\rho_{N-1,N-m-1}(0)-f_{Nm10}A_{N,N-m}^{(N,N-m)}}{f_{Nm11}}$

where in the second line we have made use of the result from $i=0$. For the
general case we find

$\displaystyle A_{N-i,N-m-i}^{(N,N-m)}$ $\displaystyle=$
$\displaystyle\begin{cases}\rho_{N,N-m}(0)&i=0,\\\
\frac{1}{f_{Nmii}}\left(\rho_{N-i,N-m-i}(0)-\sum\limits_{i^{\prime}=0}^{i-1}f_{Nmii^{\prime}}A_{N-i^{\prime},N-m-i^{\prime}}^{(N,N-m)}\right)&i>0.\end{cases}$
(86)

The coefficients generated from Eq. (86) can then be substituted into Eq. (83)
to give the analytic solution for the full time evolution of the density
matrix elements.

This recipe scales poorly ($\sim\alpha_{0}^{6}$), is tedious to derive, fails
numerically for the closed system ($f_{Nmii^{\prime}}$ develops
singularities), and is memory-intensive. However, it has the nice property
that, since it front-loads the difficulty of the problem into calculating the
matrix of coefficients $A$, it is intrinsically stable in time. With this
recipe, we have generated a movie giving the qualitative features of the
evolution of the bosonic anharmonic oscillator for a small value of
$\alpha_{0}$ Propp (2023a).

Analytic solutions to non-linearly evolving open quantum system dynamics have
thus far been rare to find in quantum optics. It is our hope that this method,
though cumbersome, may inspire others to find more solutions.

##  Appendix C Kitten State Expectation Values

In this section, we derive the three cases in Eq. (IV) explicitly, solidifying
the relationship between the $N$-kitten states and the operators
$\hat{X}_{\theta}^{n}$ via their shared symmetry $\mathbb{Z}_{n}$ (the
$n$-fold cyclic group symmetry).

Recall from Eq. (33) that an operator $\hat{X}_{\theta}^{n}$ is decomposable
into a sum of symmetrically ordered operators
$\\{\hat{a}^{\dagger\mu}\hat{a}^{\nu}\\}_{\rm sym}$, where the difference
$|\mu-\nu|\leq n$ for each term. By the results of Appendix B, recurrence
times are independent of operator ordering and open system effects, so it will
suffice to study normally ordered operators
$\hat{a}^{\dagger\mu}\hat{a}^{\nu}$ of the closed quantum system.

The question now is, given an $N$-kitten state
$\left|\psi(t_{M,N})\right\rangle$ as defined in Eq. (II), what is
$\langle\hat{a}^{\dagger\mu}\hat{a}^{\nu}\rangle$ and when is the
semiclassical TWA a reasonable approximation to it? If, for some $\mu$ and
$\nu$, this expectation is independent of $M$ and $N$, then that expectation
value must be entirely calculable in terms of semiclassical behavior since it
does not depend on the quantum structure of the state (i.e. which kitten state
has formed). To answer this question, it is sufficient to consider the
following three cases: a) $\mu=\nu$, b) $\mu-\nu=pN$, and c) $\mu-\nu\neq pN$
where $p$ is a nonzero integer.

Case a) $\mu=\nu$: Since $[\hat{H}_{\rm Kerr},\hat{n}]=0$, all operators
proportional to powers of the number operator are conserved in the evolution
$\left|\psi(t_{M,N})\right\rangle=e^{i\hat{H}_{\rm
Kerr}t_{M,N}}\left|\psi(0)\right\rangle$. Thus we conclude
$\langle\hat{a}^{\dagger\mu}\hat{a}^{\nu}\rangle=|\alpha_{0}|^{2\mu}$ when
$\mu=\nu$. This is an entirely semiclassical quantity (and easy to calculate,
as it does not even require the semiclassical approximation of the unitary
dynamics).

Case b) $\mu-\nu=pN$: First note that $\left|\psi(t_{M,N})\right\rangle$ is an
eigenstate of $\hat{a}^{N}$ with eigenvalue $\pm\alpha_{0}^{N}$, where the
sign is negative only if $N$ is even and $M$ is odd. To see this, consider
applying $\hat{a}^{N}$ to a single term in the superposition. Each application
of $\hat{a}$ multiplies the coherent state by the complex number:
$\alpha_{0}e^{\frac{i\pi(2k+M)}{N}}$ if $N$ is even and
$\alpha_{0}e^{\frac{i\pi 2k}{N}}$ if $N$ is odd. After $N$ applications the
amplitude of this number is independent of $k$ with the sign dependence noted
above.

Thus, if the difference of $\mu$ and $\nu$ is an integer multiple of $N$, we
can utilize the eigenproperty to pull out $\mu-\nu$ powers of $\alpha_{0}^{*}$
and are left with an expectation value of a power of the number operator as in
the first case. Explicitly,

$\langle\hat{a}^{{\dagger}\mu}\hat{a}^{\nu}\rangle=(\pm
1)^{p}\alpha_{0}^{*\mu-\nu}|\alpha_{0}|^{2\nu}\qquad\mu-\nu=p\,N,$ (87)

where, again, the sign is negative if $N$ is even and $M$ is odd and positive
otherwise.

This quantity is $M$ and $N$ dependent and thus will not be captured by
semiclassical calculations that are insensitive to the particular quantum
state. Hence, the TWA will not give the correct value here.

Case c) $\mu-\nu\neq pN$: First we show the following. For
$2\pi|\alpha_{0}|\gg N$, the expectation value
$\left\langle\psi(t_{M,N})\right|\hat{a}^{q}\left|\psi(t_{M,N})\right\rangle$
tends to zero for $0<q<N$ and is exactly zero for all $N$ in the
infinite-$|\alpha_{0}|$ limit. Here the integer $q$ plays the role of $pN$
(which must be an integer since $\mu-\nu$ is an integer). For now we assume
$p<1$ so that $q<N$, but this restriction will be relaxed later.

The expectation value
$\left\langle\psi(t_{M,N})\right|\hat{a}^{q}\left|\psi(t_{M,N})\right\rangle$
is rewritten in terms of an analytic function $f(\omega)$, such that

$\displaystyle\left\langle\psi(t_{M,N})\right|\hat{a}^{q}\left|\psi(t_{M,N})\right\rangle=\tilde{\alpha}_{0}^{q}f(\omega_{N}^{q})+\mathscr{O}(\left\langle\alpha_{i}|\alpha_{j}\right\rangle),$
(88)

where $\omega_{N}$ is an $N$th root of unity $e^{2\pi i/N}$ and we have
defined $\tilde{\alpha}_{0}=e^{\frac{i\pi M}{N}}\alpha_{0}$ for $N$-even
$M$-odd and $\tilde{\alpha}_{0}=\alpha_{0}$ otherwise.
$\mathscr{O}(\left\langle\alpha_{i}|\alpha_{j}\right\rangle)$ goes to zero in
the limit of orthogonal coherent states (that is, $2\pi|\alpha_{0}|\gg N$ for
coherent states distributed uniformly around the circle with radius
$|\alpha_{0}|$). $f(\omega)$ is the sum of a finite geometric sequence

$\displaystyle
f(\omega)=1+\omega+\omega^{2}+\omega^{3}+\dots+\omega^{N-1}=\frac{1-\omega^{N}}{1-\omega}.$
(89)

Figure 7: The cyclic group $\mathbb{Z}_{n}$ for $n=6$. Initially, all states
have the same phase, marked by the position on the circle by a blue $\theta$
(as all the phases are measured with respect to this un-moving point). After a
single application of $\hat{a}$, they will have phases equally distributed
over the circle at the location of the black dots. After a second application
of $\hat{a}$, they will now be equally distributed around the three pink
markers. After a third application, they will be equally distributed around
two points: $\theta$ and its antipode. After a fourth application, they will
be equally distributed around the three green markers. And after a fifth
application, they will be equally distributed around the black dots again
before finally “recohering” after a sixth ($n$th) application at the dot with
the blue theta. At each stage the group action maintains an equal distribution
over the support: the roots of unity corresponding to each way of factoring
$n$. For a video describing this, see Propp (2023b).

Since $\omega^{qN}=1^{q}=1$, it follows that $f(\omega_{N}^{q})=0$. Hence, the
corresponding term in the expectation value in Eq. (88) is zero as well.

This result can also be understood geometrically, considering the unit circle
and the cyclic group $\mathbb{Z}_{N}$ as illustrated in Fig. 7 (as well as in
the accompanying video Propp (2023b)); at each stage of the group action
(rotating the phase of each state vector by the phase of its coherent state),
an initially equal distribution over the support remains equally distributed,
even as the structure of the support changes due to common factors within the
roots of unity.

Both this property and the eigenvalue relation above will also apply to the
associated mixed states distributed uniformly in the same manner. The only
difference is that now this property holds exactly even outside the
large-$\alpha_{0}$ limit. Since the density matrix is diagonal, there are no
contributions from terms
$\mathscr{O}(\left\langle\alpha_{i}|\alpha_{j}\right\rangle)$ as there are for
the superposition.

We now know the expectation value in Eq. (88) is approximately zero for $q<N$.
From this, we can conclude that the expectation value
$\left\langle\hat{a}^{\dagger\mu}\hat{a}^{\nu}\right\rangle$ where
$\mu-\nu=pN$ with $p\notin\mathbb{Z}$ is also approximately zero using the
eigenstate property. Taking $\mu<\nu$, we rewrite the expectation value as
$\left\langle\hat{a}^{\dagger\mu}\hat{a}^{\nu-\mu}\hat{a}^{\mu}\right\rangle$.
The action of $\hat{a}^{\dagger\mu}$ to the left and $\hat{a}^{\mu}$ to the
right result in matching rotations of the phases in the kitten state and an
overall multiplication of the state by $|\alpha_{0}|^{2\mu}$. Thus, the
expectation value reduces to
$|\alpha_{0}|^{2\mu}\left\langle\hat{a}^{\nu-\mu}\right\rangle$, and we can
make use of our results from Eq. (88) with $q=\nu-\mu$ to say that this zero,
provided $\nu-\mu<N$. If instead $\nu-\mu>N$, we make use of the eigenstate
property to pull out $N$ powers of the annihilation operator $\hat{a}$ and
replace them with $\tilde{\alpha}_{0}^{N}$, repeating the process $w$ times
until there are $q=\nu-\mu-wN<N$ powers of $\hat{a}$ remaining in the
expectation value and a pre-factor $\tilde{\alpha}_{0}^{wN}$. Since $\nu-\mu$
is not an integer multiple of $N$, $q=\nu-\mu-wN$ is not an integer multiple
of $N$, and the expectation value is approximately zero. Similarly, for
$\mu>\nu$, taking complex conjugation to give the associated expressions.

Thus, for a sufficiently large $|\alpha_{0}|$, the expectation value is
approximately zero and can be calculated semiclassically. Note that when the
number of kittens $N$ becomes large enough such that $2\pi|\alpha_{0}|\lesssim
N$, their non-orthogonality contributes to oscillations, giving rise to the
non-precision of the recurrences (due to neighboring high-$N$ kitten states)
as well as the deviations from the TWA at short times (also high-$N$ kitten
states) in Fig. 4.

To summarize, we have rederived the three cases of Eq. (IV). We know that the
recurrences for normally ordered expectation values given in Eq. (26) that we
calculated in the Heisenberg picture is a manifestation of the symmetry of
kitten states being that of the symmetric group $\mathbb{Z}_{n}$, which is
shared by operators $\hat{a}^{\dagger\mu}\hat{a}^{\mu}$ when $\mu-\nu=pN$,
with $p\in\mathbb{Z}^{\pm}$. By the various results of the previous sections,
we can also appreciate that this is true of symmetric products of powers of
$\hat{a}^{\dagger}$ and $\hat{a}$. The operators $\hat{X}_{\theta}^{n}$, when
expanded in such products, will contain terms with $\mu-\nu=\pm n$, hence the
relationships between the operators $\hat{X}_{\theta}^{n}$ and the $N$-kitten
states when $n=pN$.