# Driven-dissipative quantum Kerr resonators: new exact solutions, photon
blockade and quantum bistability

David Roberts  Department of Physics, University of Chicago, Chicago, IL
60637, USA  Pritzker School for Molecular Engineering, University of Chicago,
5640 South Ellis Avenue, Chicago, Illinois 60637, U.S.A.  Aashish A. Clerk
Pritzker School for Molecular Engineering, University of Chicago, 5640 South
Ellis Avenue, Chicago, Illinois 60637, U.S.A.

(August 27, 2024)

###### Abstract

We present a new approach for deriving exact, closed-form solutions for the
steady state of a wide class of driven-dissipative nonlinear resonator that is
distinct from more common complex-$P$ function methods. Our method generalizes
the coherent quantum absorber approach of Stannigel et al. Stannigel  _et al._
(2012) to include nonlinear driving and dissipation, and relies crucially on
exploiting the Segal-Bargmann representation of Fock space. Our solutions and
method reveal a wealth of previously unexplored observable phenomena in these
systems, including new generalized photon-blockade and anti-blockade effects,
and an infinite number of new parameter choices that yield quantum
bistability.

##  I Introduction

Exact solutions of interacting, driven-dissipative quantum problems are rare,
and thus occupy a special place in the study of open quantum systems. A
canonical example is the solution of the driven-dissipative Kerr resonator.
Here, a bosonic mode with a Kerr nonlinearity (i.e. a Hubbard $U$ interaction)
is subject to a coherent linear drive and Markovian single photon loss. As
shown by Drummond and Walls Drummond and Walls (1980), one can exactly solve
for the steady state of this system using a complex-$P$ phase space
representation. Later work showed that models including two-photon driving and
loss are also solvable using this technique Drummond and Walls (1980); Bartolo
_et al._ (2016); Elliott and Ginossar (2016). These driven nonlinear cavity
systems have renewed relevance, as they can be directly implemented in
superconducting circuit QED setups (see, e.g., Kirchmair  _et al._ (2013);
Leghtas  _et al._ (2015); Touzard  _et al._ (2018); Lescanne  _et al._ (2019);
Grimm  _et al._ (2019)). Their ability to exhibit multiple steady states has
utility in quantum information processing Mirrahimi  _et al._ (2014); Goto
(2016); Puri  _et al._ (2017).

While the existence of exact solutions here are remarkable, they are somewhat
physically opaque and unwieldy (e.g. they are typically expressed as infinite
sums of special functions). Their derivation is also somewhat intricate,
requiring a non-trivial integration to relate the solution of an effective
classical problem to the underlying quantum system. More direct methods for
obtaining and possibly extending these solutions are thus highly desirable.
For the simplest version of the Kerr-cavity problem (single-photon drive and
loss only), Stannigel et al. Stannigel  _et al._ (2012) were able to reproduce
the exact solution of Ref. Drummond and Walls (1980) using a simple, purely
algebraic approach. While extremely elegant, it was unclear whether this
approach could be extended to more complex problems.

Figure 1: (a) Generalized driven Kerr cavity problem, where a single
interacting bosonic mode is subject to linear and nonlinear coherent drives
$\Lambda_{j}$, as well as independent one and two photon loss (rates
$\kappa_{1}$, $\kappa_{2}$). (b) The coherent quantum absorber (CQA) method
represents each dissipative bath as a chiral waveguide, and couples a second
auxiliary $b$ cavity downstream. By picking its Hamiltonian judiciously, the
entire composite system can relax to a pure state, providing an efficient
means for finding the steady state of cavity-$a$.

In this paper, we show that such an extension is indeed possible: the
“coherent quantum absorber” (CQA) method of Ref. Stannigel  _et al._ (2012)
can be extended to a wide class of systems which include nonlinear coherent
driving as well as multiple dissipators (see Fig. 1). Our extension employs a
new ingredient: the Segal-Bargmann representation of a single-mode pure state
wavefunction Bargmann (1961); Segal (1962, 1963); Bargmann (1967). This
enables non-trivial transformations that are crucial for finding exact
solutions. Our approach yields several new insights. We find and describe new
parameter regimes where the steady-state exhibits a surprising generalized
photon-blockade phenomenon. In particular, we show how the use of a nonlinear
driving term allows for photon blockade even for nonlinearities much weaker
than dissipation rates; unlike so-called “unconventional photon blockade" Liew
and Savona (2010); Bamba  _et al._ (2011); Lemonde  _et al._ (2014), the
effect we describe results in non-Gaussian states and a complete suppression
of higher-$n$ photon probabilities.

We also find an infinite number of points in parameter space where our
generalized driven Kerr system exhibits quantum bistability (i.e. a two-
dimensional decoherence-free subsystem), despite the lack of photon-number
parity conservation. The required parameters can be achieved asymptotically in
the limit of weak single-photon loss. Our solution also provides a simple
intuitive picture when there is a unique steady state (see Fig. 2): the steady
state is formed by mixing a pure state with vacuum at a 50-50 beamsplitter,
and then discarding one of the outputs. At a technical level, our work also
provides new, simple closed-form expressions for the steady-state Wigner
function and normally-ordered moments.

The remainder of this paper is organized as follows. In Sec. II we introduce
the basic system. In Sec. III, we review the CQA method, and in Sec. IV
present our extension to nonlinear driving and multiple dissipators. Sec. V
summarizes the new physical phenomena uncovered by our exact solution, while
Sec. VI and Sec. VII discuss regimes of classical and quantum bistability. We
conclude in Sec. VIII.

##  II System

We consider a driven Kerr resonator whose coherent dynamics is described by
the Hamiltonian

$\displaystyle\hat{H}_{a}$
$\displaystyle=\frac{K}{2}\hat{a}^{\dagger}\hat{a}^{\dagger}\hat{a}\hat{a}-\Delta\hat{a}^{\dagger}\hat{a}$
$\displaystyle+\bigg{[}\bigg{(}\Lambda_{1}\hat{a}^{\dagger}+\frac{\Lambda_{2}}{2}\hat{a}^{\dagger}\hat{a}^{\dagger}+\Lambda_{3}\hat{a}^{\dagger}\hat{a}^{\dagger}\hat{a}\bigg{)}+h.c.\bigg{]}.$
(1)

We work in a rotating frame, and have assumed that all drives have an equal
detuning $\Delta$ from the cavity resonance frequency (which allows us to have
a time-independent rotating-frame Hamiltonian). Here $K$ is the Kerr
nonlinearity, and $\Lambda_{1},\Lambda_{2}$ are the complex amplitudes of
standard coherent one and two photon driving terms. $\Lambda_{3}$ represents
an unusual kind of nonlinear single-photon driving term; as we will see, it
enables a striking new kind of photon blockade effect that does not require
strong nonlinearity. We show in Appendix A how this $\Lambda_{3}$ drive can be
implemented using the superconducting circuit architecture of Ref. Sivak  _et
al._ (2019).

The full dissipative dynamics includes one and two photon loss processes, and
is described by the Lindblad master equation

$\displaystyle\frac{d}{dt}\hat{\rho}$
$\displaystyle=-i[\hat{H}_{a},\hat{\rho}]+\kappa_{1}\mathcal{D}[\hat{a}]\hat{\rho}+\kappa_{2}\mathcal{D}[\hat{a}^{2}]\hat{\rho}\equiv\mathcal{L}_{0}\hat{\rho},$
(2)

where
$\mathcal{D}[\hat{X}]\hat{\rho}\equiv\hat{X}\hat{\rho}\hat{X}^{\dagger}-(1/2)\left\\{\hat{X}^{\dagger}\hat{X},\hat{\rho}\right\\}$
is the usual Lindblad dissipative superoperator, and $\kappa_{1}$
($\kappa_{2}$) are the one (two) photon decay rates. Note that the dissipative
evolution corresponds to coupling the system to two distinct zero-temperature
baths.

We will focus exclusively on finding the steady states of this kind of system,
i.e. density matrices $\hat{\rho}_{\rm ss}$ satisfying

$\mathcal{L}_{0}\hat{\rho}_{\rm ss}=0.$ (3)

We briefly summarize prior work on this model. For $\Lambda_{3}\equiv 0$,
exact solutions for $\hat{\rho}_{\rm ss}$ have been found using the complex
$P$-function approach Drummond and Walls (1980); Drummond  _et al._ (1981);
Elliott and Ginossar (2016); Bartolo  _et al._ (2016). The solutions express
matrix elements of $\hat{\rho}_{\rm ss}$ in the Fock basis as sums of special
functions. In the semiclassical limit, solutions for the steady state can be
found using an alternate approach developed by Dykman and co-workers Marthaler
and Dykman (2006); Dykman (2012); unlike the complex-$P$ approach, these can
also be used to describe dissipation at a non-zero temperature. Systems with
higher-order coherent driving terms (like our $\Lambda_{3}$) have been studied
previously, (see, e.g., Guo  _et al._ (2013); Svensson  _et al._ (2018); Lörch
_et al._ (2019)), but were not previously known to be solvable.

While the prior work on driven Kerr resonators is a remarkable achievement, it
leaves several mysteries unanswered. First, in the presence of single photon
loss, the unique steady state that one finds always yields a positive-definite
Wigner function. Given the nonlinearity in the system, it is not a priori
obvious that this should necessarily be the case. Second, in the absence of
single photon processes (i.e. $\Lambda_{1}=\kappa_{1}=0$), this system
exhibits multiple steady states Krippner  _et al._ (1994); Wolinsky and
Carmichael (1988); Hach III and Gerry (1994); Gilles  _et al._ (1994);
Mirrahimi  _et al._ (2014). We are not aware of any discussion of this using
the complex-$P$ approach. For $\Delta=0$, the system is simple enough that the
multiple steady states can be found via elementary means, in terms of
superpositions of coherent states Goto (2016); Puri  _et al._ (2017).
Conditions needed for Wigner function negativity were recently discussed in
Ref. Braasch  _et al._ (2019), though these are not directly applicable to our
system. In the sections that follow, we discuss an alternate, physically-
transparent method for solving this class of problems that addresses the open
issues mentioned above.

##  III Exact solutions using the quantum absorber method

Our approach to solving driven-dissipative Kerr problems is to adapt and
extend the so-called “coherent quantum absorber” (CQA) approach first
introduced by Stannigel et al. Stannigel  _et al._ (2012) to solve the
simplest driven Kerr problem where there are no two photon drive or loss
processes. We quickly recap the philosophy of this approach, and then show how
it can be extended to deal with more complex problems involving two and even
three photon processes.

###  III.1 Recap of the basic approach

Consider first the case where our system in Eq. (2) has no two photon loss
($\kappa_{2}=0$). The starting point of the CQA method is to represent the one
photon loss as arising from a coupling to a chiral (i.e. unidirectional)
waveguide. Further, one imagines coupling a second auxiliary bosonic mode
(annihilation operator $\hat{b}$, system Hamiltonian $\hat{H}_{b}$) to the
waveguide, downstream from the physical $a$ cavity (see Fig. 1). Given the
chirality of the waveguide, the dynamics of this auxiliary cavity can have no
impact on the physical cavity $a$. The entire composite system can be
described using standard cascaded quantum systems theory Carmichael (1993);
Gardiner (1993); Gardiner and Zoller (2000). The dynamics of the reduced
density matrix $\hat{\rho}_{ab}$ describing both cavities is described by a
Lindblad master equation of the form:

$\displaystyle\frac{d}{dt}\hat{\rho}_{ab}$
$\displaystyle=-i[\hat{H}_{ab},\hat{\rho}_{ab}]+\kappa_{1}\mathcal{D}[\hat{a}-\hat{b}]\hat{\rho}_{ab},$
(4) $\displaystyle\hat{H}_{ab}$
$\displaystyle=\hat{H}_{a}+\hat{H}_{b}-\frac{i\kappa_{1}}{2}(\hat{a}^{\dagger}\hat{b}-h.c.).$
(5)

Note that one can rigorously trace out cavity $b$ from this equation,
recovering Eq. (2) for cavity $a$ alone.

While the introduction of the auxiliary cavity $b$ has no impact on cavity
$a$, it provides a useful tool for finding its steady state. As shown in Ref.
Stannigel  _et al._ (2012), for a general cavity $a$ Lindblad master equation
having only single-photon loss (i.e. Eq. (2) with $\kappa_{2}=0$ and arbitrary
$\hat{H}_{a}$), one can always construct a Hamiltonian $\hat{H}_{b}$ for the
auxiliary cavity $b$ such that the composite system has a pure steady state.
This steady state state necessarily has vanishing emission to the waveguide–
it is a “dark” state. Letting $\hat{\rho}_{ab,ss}$ denote the steady-state
density matrix of the two-cavity problem, this means:

$\displaystyle\hat{\rho}_{ab,ss}$
$\displaystyle=\left|\psi\right>\left<\psi\right|,\hskip
17.07182pt\left(\hat{a}-\hat{b}\right)\left|\psi\right>=0.$ (6)

Note that the dark state condition implies that $\left|\psi\right>$ is
essentially a single mode state. Introducing new composite mode operators

$\hat{c}_{\pm}\equiv\frac{\hat{a}\pm\hat{b}}{\sqrt{2}},$ (7)

one notes that the dark state condition forces the composite mode
$\hat{c}_{-}$ to be in vacuum. Hence, one just needs to solve for the (pure)
state of the composite $\hat{c}_{+}$ mode.

In physical terms, the CQA approach seeks to construct $\hat{H}_{b}$ such that
the auxiliary cavity $b$ acts as a “perfect absorber” for all photons emitted
into the waveguide by cavity $a$. By tracing out cavity $b$, one obtains the
desired steady state for the physical cavity-$a$ problem. One generically
obtains an impure state, as the two cavities will be entangled in the state
$\left|\psi\right>$.

While such a construction is always possible, in practice it would seem to be
of no utility, as one can only construct the required $\hat{H}_{b}$ by first
independently solving for the cavity-$a$ steady $\hat{\rho}_{a,ss}$, Despite
this seeming obstacle, Ref. Stannigel  _et al._ (2012) demonstrated that for a
range of problems, one could essentially guess the form of $\hat{H}_{b}$
without first knowing $\hat{\rho}_{a,ss}$. This educated guess is extremely
simple: $\hat{H}_{b}$ is taken to be identical to $\hat{H}_{a}$ up to an
overall minus sign. Ref. Stannigel  _et al._ (2012) applied this to the
simplest driven Kerr problem ($\Lambda_{2}=\Lambda_{3}=\kappa_{2}=0$ in Eq.
(138)), in which case

$\hat{H}_{b}=-\frac{K}{2}\hat{b}^{\dagger}\hat{b}^{\dagger}\hat{b}\hat{b}+\Delta\hat{b}^{\dagger}\hat{b}-\left[\Lambda_{1}\hat{b}^{\dagger}+h.c.\right].$
(8)

With this choice, Stannigel et al. were able to find a pure-state solution of
the cascaded master equation in Eq. (4) by solving a simple one-term recursion
relation. By then tracing out cavity $b$, they recovered (in a much simpler
manner) the classic solution of Drummond et al. Drummond and Walls (1980) for
the linear-drive Kerr problem.

###  III.2 Extension to nonlinear driving and two-photon loss

It is natural to ask whether the absorber method approach can be extended to
solve problems with nonlinear driving and two-photon loss. An immediate issue
is the presence of two independent dissipators in the master equation Eq. (2).
We find that the CQA approach is easily modified to deal with this situation.
As shown in Fig. 1(b), one can represent the two-photon loss process as a
nonlinear coupling to a second chiral waveguide.

One again needs to add something downstream along this waveguide to absorb the
emitted excitations. While there are many possible options, we find the
simplest approach is sufficient: we assume that there is still a single
auxiliary cavity $b$ that now couples to both these independent chiral
waveguides. The cascaded master equation now takes the form:

$\displaystyle\frac{d}{dt}\hat{\rho}_{ab}$
$\displaystyle=-i[\hat{H}_{ab},\hat{\rho}_{ab}]+\kappa_{1}\mathcal{D}[\hat{a}-\hat{b}]\hat{\rho}_{ab}+\kappa_{2}\mathcal{D}[\hat{a}^{2}-\hat{b}^{2}]\hat{\rho}_{ab},$
(9)

with

$\displaystyle\hat{H}_{ab}\equiv\hat{H}_{a}-\hat{H}_{b}-\frac{i\kappa_{1}}{2}(\hat{a}^{\dagger}\hat{b}-h.c.)-\frac{i\kappa_{2}}{2}(\hat{a}^{\dagger}\hat{a}^{\dagger}\hat{b}^{2}-h.c.).$
(10)

Again, tracing out cavity $b$ from the above equation recovers the cavity $a$
master equation given in Eq. (2), independent of the choice of $\hat{H}_{b}$.

The next step is the same as before: we want to pick $\hat{H}_{b}$ so that
cavity $b$ absorbs all photons emitted by cavity $a$ into either of the two
chiral waveguides. We thus want a pure steady state $\left|\psi\right>$ of the
two cavity system that is a dark state of both collective loss operators
appearing in Eq. (9). Fortunately, these dark state conditions are not
independent: having $(\hat{a}-\hat{b})\left|\psi\right>=0$ as before ensures
that the state is dark with respect to emission to either waveguide.

Finally, there remains the question of how exactly to find the desired
$\hat{H}_{b}$. As we show in Sec. IV, the simple educated guess of taking
$\hat{H}_{b}$ to be the negative of $\hat{H}_{a}$ still works in the presence
of two photon driving and loss, and even for a wider class of problems.

###  III.3 Connection to Segal-Bargmann representations

A second crucial element in our extension of the CQA method is to combine it
with the Segal-Bargmann (SB) representation of single-mode pure-state
wavefunctions in terms of holomorphic functions Bargmann (1961); Segal (1962,
1963); Bargmann (1967). This provides an extremely efficient way of solving
the complex recursion relations that determine the desired dark state
wavefunction $\left|\psi\right>$. More importantly, it is an extremely useful
tool for developing physical intuition. It renders the operation of tracing
out the auxiliary cavity $b$ trivial, and allows one to directly obtain the
Wigner function of the cavity-$a$ steady state.

####  III.3.1 Basics of the representation

Consider a single bosonic mode in a pure state $\left|\psi\right>$ that is
written in terms of Fock states $\left|m\right>$ as:

$\left|\psi\right>=\sum_{m=0}^{\infty}\alpha_{m}|m\rangle.$ (11)

In the SB representation, this state is associated with a holomorphic function
$\psi_{\rm SB}(z)$ defined on the complex plane:

$\psi_{\rm SB}(z)=\sum_{m=0}^{\infty}\frac{\alpha_{m}}{\sqrt{m!}}z^{m}.$ (12)

The space of these functions forms a Hilbert space that is unitarily
equivalent to the original Fock space, with an induced inner product:

$\displaystyle\langle\psi_{\rm SB},\phi_{\rm SB}\rangle_{{\rm
SB}}\equiv\frac{1}{\pi}\int_{\mathbb{C}}d^{2}z\,\psi_{\rm SB}^{*}(z)\phi_{\rm
SB}(z)e^{-|z|^{2}}.$ (13)

The SB wavefunction has a direct physical interpretation: its modulus
determines the Husimi $Q$-function of the state $\left|\psi\right>$. Letting
$\left|z\right>$ denote a coherent state with amplitude $z$, we have

$Q(z)\equiv\frac{1}{\pi}\left|\left<z\vphantom{\psi}\right|\left.\psi\vphantom{z}\right>\right|^{2}=\frac{1}{\pi}\left|\psi_{\rm
SB}(z^{*})\right|^{2}e^{-|z|^{2}}.$ (14)

Finally, the canonical $\hat{c}$ and creation $\hat{c}^{\dagger}$ operators
become linear differential operators in the Bargmann space:
$\hat{c}\mapsto\partial/\partial z$, $\hat{c}^{\dagger}\mapsto z$.

Figure 2: Simple picture of the unique steady state of the generalized driven
Kerr resonator with non-zero single-photon loss. One starts with a pure state,
single-mode wavefunction $\left|\psi_{+}\right>$. This is mixed with vacuum
noise at a 50-50 beamsplitter; the output ports represent the final steady
state of the physical $a$ cavity and the auxiliary $b$ cavity. This operation
implies that the cavity-$a$ steady state is $\left|\psi_{+}\right>$ convolved
with vacuum noise. As a result, cavity-$a$’s steady-state Wigner function
$W_{a}(z)$ is equal (up to scaling) to the $Q$ function $Q_{\psi_{+}}(z)$ of
the pure state $\left|\psi_{+}\right>$.

####  III.3.2 Tracing out the auxiliary cavity

As we will see, the CQA method reduces to finding a single-mode, pure-state
wavefunction $\left|\psi_{+}\right>$ for the collective mode
$\hat{c}_{+}=(\hat{a}+\hat{b})/\sqrt{2}$; the orthogonal mode $\hat{c}_{-}$
must be in vacuum to have a dark state. To find the corresponding state of the
physical cavity $a$, one transforms from the $\hat{c}_{\pm}$ basis to the
$\hat{a}/\hat{b}$ basis, and then traces out the state of the auxiliary cavity
$b$. This operation has a very simple physical interpretation (see Fig. 2): it
corresponds to mixing the state $\left|\psi_{+}\right>$ with vacuum noise at a
50-50 beamsplitter, and then discarding one of the output modes.

At a heuristic level, this operation implies that in phase space, the
cavity-$a$ steady state will be equivalent to that of the state
$\left|\psi_{+}\right>$ convolved with an extra half-quantum of vacuum noise.
Recall that this is the same transformation that converts a Wigner function
into a $Q$-function. As a result, we find a very simple expression for the
cavity-$a$ steady state Wigner function. Letting $\psi_{+,{\rm SB}}(z)$ denote
the SB representation of the pure state $\left|\psi_{+}\right>$, we have:

$W_{a}(z)=2Q_{+}(\sqrt{2}z)=\frac{2}{\pi}\left|\psi_{+,{\rm
SB}}\left(\sqrt{2}z^{*}\right)\right|^{2}e^{-2|z|^{2}}$ (15)

We see that the SB “wavefunction” $\psi_{+,{\rm SB}}(z)$ has a direct physical
interpretation: its modulus determines the cavity-$a$ steady-state Wigner
function. We also see that this Wigner function must necessarily be positive
(as it is equivalent to the $Q$-function of the state $\left|\psi_{+}\right>$,
and the $Q$ function is always positive). The above relation follows from the
fact that Wigner functions transform in the expected way (i.e. like classical
probability distributions) under a beamsplitter transformation.

A more direct relation results from examination of the cavity-$a$ steady-state
$P$-function: since the $P$-function is insensitive to vacuum noise, the
output of the beamsplitter in the $P$-representation is simply equal to the
$P$-function $P_{+}(z)$ of the state $|\psi_{+}\rangle$, rescaled by a factor
of $\sqrt{2}$:

$\displaystyle P_{a}(z)=2P_{+}(\sqrt{2}z).$ (16)

This immediately implies that the cavity-$a$ $P$ function is generically
singular and thus non-positive Cahill (1969). Thus, in the pathological
behavior of the cavity-$a$ $P$-function, we have a tell-tale signature of
nonclassicality.

##  IV CQA solution of the general driven Kerr cavity

We now use the results of Sec. III to solve Eq. (2) for a driven Kerr
resonator subject to both one and two photon driving, and one and two-photon
loss. This will allow us to reproduce results previously derived using
complex-$P$ methods Drummond  _et al._ (1981); Bartolo  _et al._ (2016);
Minganti  _et al._ (2016); Elliott and Ginossar (2016), but in a manner that
allows greater physical intuition. We are also able to solve an extended model
which includes a nonlinear single-photon driving term; this model has not been
previously solved. Our approach yields several new physical insights: the
possibility of photon blockade and “anti-blockade” phenomena, and the
possibility of near quantum-bistability without parity conservation. We focus
in this section on the case where there is non-zero single-photon loss
($\kappa_{1}\neq 0$), implying the existence of a unique steady state. In Sec.
VII, we turn to the case where there is no single photon driving or loss; we
are able to use the CQA method to provide insights into the bistability in
this system, and how this changes from quantum to classical bistability with
the addition of a drive detuning.

###  IV.1 Solution without nonlinear single-photon driving

We are interested in the driven-Kerr system described by Eqs. (138) and (2)
with $\Lambda_{3}=0$ and $\kappa_{1}>0$. The CQA approach represents this
system using the equivalent two-cavity cascaded system in Eqs. (9) and (10).
We seek a pure-state steady-state $\left|\psi\right>$ that is necessarily dark
with respect to dissipation, meaning that Eq. (6) is satisfied:
$\sqrt{2}\hat{c}_{-}\left|\psi\right>\equiv(\hat{a}-\hat{b})\left|\psi\right>=0$.
Our steady state can thus be written as a tensor product of a non-trivial
state of the $\hat{c}_{+}$ collective mode, and a vacuum state for the
$\hat{c}_{-}$ mode:

$\left|\psi\right>=\left|\psi_{+}\right>\left|0_{-}\right>.$ (17)

In order for $\left|\psi\right>$ to be a steady state, it also needs to be an
eigenstate of the cascaded Hamiltonian $\hat{H}_{ab}$ with energy $E$. Writing
$\hat{H}_{ab}$ in terms of $\hat{c}_{\pm}$, and using the fact that
$\hat{c}_{-}\left|\psi\right>$ vanishes, the eigenvalue equation becomes

$\frac{1}{2}\hat{c}^{\dagger}_{-}\hat{\mathcal{H}}_{+}\left|\psi\right>=E\left|\psi\right>,$
(18)

with

$\displaystyle\hat{\mathcal{H}}_{+}$
$\displaystyle\equiv(K-i\kappa_{2})\hat{c}_{+}^{\dagger}\hat{c}_{+}^{2}-(2\Delta+i\kappa_{1})\hat{c}_{+}+2\Lambda_{2}\hat{c}_{+}^{\dagger}+2\sqrt{2}\Lambda_{1}.$
(19)

Our choice of the auxiliary-cavity Hamiltonian $\hat{H}_{b}$ thus leads to a
cascaded Hamiltonian that necessarily creates an excitation in the
$\hat{c}_{-}$ mode. It follows that we must have $E=0$. Having
$\left|\psi\right>$ be a stationary state then reduces to a single mode
problem:

$\hat{\mathcal{H}}_{+}\left|\psi_{+}\right>=0,$ (20)

i.e. we need to find a pure state $\left|\psi_{+}\right>$ that is annihilated
by the non-Hermitian operator $\hat{\mathcal{H}}_{+}$.

The seemingly obvious next step is to follow the approach used in Ref.
Stannigel  _et al._ (2012): express $\left|\psi_{+}\right>$ in the Fock state
basis, and turn Eq. (20) into a recursion relation for the expansion
coefficients $\alpha_{j}$:

$\displaystyle\Big{[}(K-i\kappa_{2})m-(2\Delta+i\kappa_{1})\Big{]}\sqrt{m+1}\alpha_{m+1}$
(21)
$\displaystyle+2\Lambda_{2}\sqrt{m}\alpha_{m-1}+2\sqrt{2}\Lambda_{1}\alpha_{m}$
$\displaystyle=0.$

In special cases, this reduces to an easily solvable single-term recursion
relation: either the case of no two-photon driving $\Lambda_{2}=0$ Stannigel
_et al._ (2012), or the case of no one-photon driving $\Lambda_{1}=0$ Mamaev
(2018). In the more general case, the resulting two-term recursion relation is
more unwieldy.

A more direct way of getting the desired solution is to use the SB
representation $\psi_{+,{\rm SB}}(z)$ of the state $\left|\psi_{+}\right>$.
Eq. (20) is then transformed into a second-order ordinary differential
equation:

$\displaystyle\Big{[}z\frac{\partial^{2}}{\partial z^{2}}$
$\displaystyle-D\frac{\partial}{\partial
z}+(\lambda_{2}z+\lambda_{1})\Big{]}\psi_{+,{\rm SB}}(z)=0,$ (22)

where

$\displaystyle D$
$\displaystyle=\frac{2\widetilde{\Delta}}{\widetilde{K}},\leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\
\lambda_{2}=\frac{2\Lambda_{2}}{\widetilde{K}},\leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\
\lambda_{1}=\frac{2\sqrt{2}\Lambda_{1}}{\widetilde{K}}.$ (23)

Here, $\widetilde{\Delta}\equiv\Delta+i\kappa_{1}/2$ and $\widetilde{K}\equiv
K-i\kappa_{2}$ are, respectively, effective complex detuning and Kerr
nonlinearity parameters.

Without two-photon driving (i.e. $\lambda_{2}=0$), Eq. (22) is a standard
hypergeometric equation. It has a unique analytic solution:

$\displaystyle\psi_{+,{\rm SB}}(z)$
$\displaystyle=Nz^{(D+1)/2}J_{-(D+1)}\left(2\sqrt{\lambda_{1}z}\right),$ (24)

where $J_{n}(x)$ is a Bessel function and $N$ is a normalization constant.
Using the correspondence between the SB wavefunction and Fock state amplitudes
(c.f. Eq. (12)), we recover the infinite series result given in Ref. Stannigel
_et al._ (2012), which in turn corresponds to the classic solution of Ref.
Drummond and Walls (1980). The closed form we have here has additional
virtues. Via Eq. (15), it directly yields a closed form expression for the
steady-state Wigner function of the physical cavity $a$; this is in contrast
to expressions involving infinite sums that are the usual result of
complex-$P$ solutions. Our expression for this case agrees with that derived
earlier (via an alternate method) Kheruntsyan (1999).

We turn now to the more interesting case where $\lambda_{2}\neq 0$. Eq. (21)
is now a more nontrivial second-order recursion relation. The SB
representation allows us, however, to simplify the system via non-standard
transformations. An example is a “non-unitary gauge transformation”

$\displaystyle\psi_{+,\rm SB}(z)\equiv e^{-\theta(z)}\phi(z),$ (25)

where $\theta(z)$ is the “gauge potential”. This transformation shifts the
differentiation operator by the gradient of $\theta(z)$,
$\partial_{z}\mapsto\partial_{z}-\partial_{z}\theta(z)$. Here, we try the
simplest potential $\theta(z)\equiv\epsilon z$, with $\epsilon$ some constant.
Note that as $\theta$ is not purely imaginary, the resulting transformation on
the Hilbert space is non-unitary. In the Fock representation, it is equivalent
to acting on the state by the exponential of a raising operator:

$\left|\psi_{+}\right>\propto
e^{-\epsilon\hat{c}_{+}^{\dagger}}\left|\phi\right>.$ (26)

After our transformation, the problematic two-photon driving term is
effectively shifted by an amount $\epsilon^{2}$:

$\displaystyle\Big{[}z\frac{\partial^{2}}{\partial z^{2}}-$
$\displaystyle(2\epsilon z+D)\frac{\partial}{\partial z}$
$\displaystyle+(\lambda_{2}+\epsilon^{2})z+(\lambda_{1}+\epsilon
D)\Big{]}\phi(z)=0.$ (27)

It can thus be eliminated by choosing $\epsilon$ such that

$\epsilon_{\pm}=\pm i\sqrt{\lambda_{2}}.$ (28)

We will call these non-unitary gauges plus-gauge and minus-gauge. Choosing,
e.g. the plus gauge $\epsilon\equiv\epsilon_{+}$, we see that the gauge-
transformed state $\phi(z)$ satisfies Kummer’s differential equation (see
Brychkov (2008)), so that:

$\displaystyle\phi(z)=N_{0}\bigg{[}\,_{1}F_{1}\bigg{(}-\frac{\lambda_{1}+\epsilon
D}{2\epsilon};-D;2\epsilon z\bigg{)}\bigg{]},$ (29)

where $N_{0}$ is a normalization factor, and $\,{}_{1}F_{1}(r_{1};r_{2};z)$ is
Kummer’s hypergeometric function, the same special function which appears in
the hydrogen atom problem (see, e.g. Sakurai and Napolitano (2017)). We stress
that that the special case where $D$ is a positive integer must be treated
specially; this is discussed in Sec. V. Note also that in the
$\epsilon\rightarrow 0$ limit, the solution above tends smoothly to the
Bessel-function solution in Eq. (24).

The above result combined with Eq. (15) immediately yields a closed-form
expression for the steady-state Wigner function of the physical $a$ cavity of
interest:

$\displaystyle
W_{a,ss}(z)=N|\phi(\sqrt{2}z^{*})|^{2}\,e^{-2|z+\epsilon/\sqrt{2}|^{2}},$ (30)

where $N$ is a normalization constant. Note that if $\phi(z)=1$, then
$W_{a,ss}(z)$ corresponds to a coherent state with amplitude
$\alpha=\sqrt{-\lambda_{2}/2}$. Thus, a non-unity $\phi(z)$ describes
corrections to the dark state being just a simple coherent state. Note also
that if one had chosen the minus gauge in Eq. (25), one obtains an identical
solution (see Appendix D).

###  IV.2 Including nonlinear single-photon driving

We now allow $\Lambda_{3}\neq 0$ in Eq. (138). We are still able to exactly
solve for the steady state in this case; unless $\kappa_{2}=0$, it has a
qualitatively different form from the $\Lambda_{3}=0$ case. The CQA method
proceeds as in Sec. IV.1. We again write the two-mode dark state as
$|\psi\rangle=|\psi_{+}\rangle|0_{-}\rangle$, and the eigenvalue equation
again reduces to finding the kernel of a non-Hermitian operator
$\hat{\mathcal{H}}_{+}$:

$\displaystyle\hat{\mathcal{H}}_{+}=\Big{(}\widetilde{K}\hat{c}_{+}^{\dagger}+\sqrt{2}\Lambda_{3}^{*}\Big{)}\hat{c}_{+}^{2}$
$\displaystyle+\Big{(}2\sqrt{2}\Lambda_{3}\hat{c}_{+}^{\dagger}-2\widetilde{\Delta}\Big{)}\hat{c}_{+}$
$\displaystyle+\left(2\Lambda_{2}\hat{c}_{+}^{\dagger}+2\sqrt{2}\Lambda_{1}\right).$
(31)

Comparing against Eq. (19), we see that the presence of $\Lambda_{3}$ creates
a term proportional to $\hat{c}_{+}^{2}$. Attempting to solve directly for
$\left|\psi_{+}\right>$ in the Fock basis leads a complicated recursion
relation, as now we have terms that add a photon
$(\propto\hat{c}_{+}^{\dagger})$, as well as those that subtract two photons
$(\propto\hat{c}_{+}^{2})$. One obtains a third-order recursion, in place of
the second-order recursion that we had before.

One can nonetheless still solve for the dark state in closed-form. We first
perform a displacement,

$\displaystyle|\xi_{+}\rangle=\hat{D}(\alpha_{+})|\psi_{+}\rangle,$ (32)

where $\alpha_{+}=\sqrt{2}\Lambda_{3}/\widetilde{K}^{*}$, and
$\hat{D}(\alpha)\equiv e^{\alpha\hat{c}^{\dagger}_{+}-h.c.}$ is the standard
displacement operator. We can then remove the two-photon drive by applying a
non-unitary gauge transformation (as before), yielding a differential equation
which again has a simple solution in terms of Kummer’s confluent
hypergeometric function:

$\displaystyle\phi(z)$
$\displaystyle=N_{0}\bigg{[}\,_{1}F_{1}\bigg{(}-\frac{\lambda_{1}+\epsilon_{+}D}{\epsilon_{+}-\epsilon_{-}};-D;(\epsilon_{+}-\epsilon_{-})z\bigg{)}\bigg{]}.$
(33)

Here, $\epsilon_{\pm}$ correspond to the non-unitary gauge choices in which
the displaced two-photon drive vanishes (c.f. Eq. (28)):

$\displaystyle\lambda_{2}-\lambda_{3}\epsilon+\epsilon^{2}$ $\displaystyle=0$
(34)

To be manifestly consistent with the solution of the driven Kerr cavity
without nonlinear coherent driving, we have again written the solution in the
plus gauge. Finally, $\lambda_{3},\lambda_{2},\lambda_{1},D$ are the following
general complex constants:

$\displaystyle D$
$\displaystyle=\frac{2}{\widetilde{K}}\bigg{(}\widetilde{\Delta}+\frac{2|\Lambda_{3}|^{2}}{\widetilde{K}}\bigg{)},$
(35) $\displaystyle\lambda_{1}$
$\displaystyle=\frac{\sqrt{2}\Lambda_{3}}{|\widetilde{K}|^{2}}\bigg{(}\frac{4|\Lambda_{3}|^{2}}{\widetilde{K}}+2\widetilde{\Delta}\bigg{)}+\frac{2\sqrt{2}}{\widetilde{K}}\bigg{(}\Lambda_{1}-\frac{\Lambda_{2}\Lambda_{3}^{*}}{\widetilde{K}}\bigg{)},$
(36) $\displaystyle\lambda_{3}$
$\displaystyle=\frac{2\sqrt{2}\Lambda_{3}}{\widetilde{K}}\bigg{(}1-\frac{\widetilde{K}}{\widetilde{K}^{*}}\bigg{)},\leavevmode\nobreak\
\leavevmode\nobreak\
\lambda_{2}=\frac{2\Lambda_{3}^{2}}{|\widetilde{K}|^{2}}\bigg{(}\frac{\widetilde{K}}{\widetilde{K}^{*}}-2\bigg{)}+\frac{2\Lambda_{2}}{\widetilde{K}}.$
(37)

We have again defined $\tilde{\Delta}=\Delta+i\kappa_{1}/2$,
$\tilde{K}=K-i\kappa_{2}$. For the case where $\Lambda_{3}\rightarrow 0$,
these parameters revert to those given before Eq. (23). Note that for
vanishing two-photon loss, $\widetilde{K}$ is real, and hence Eq. (37) implies
that $\lambda_{3}=0$. In this case, the cubic drive does not give us anything
qualitatively new, as it can be completely eliminated by our displacement
transformation. In contrast, for non-zero $\kappa_{2}$, cubic driving gives
rise to genuinely new phenomena.

As before, the solution above directly determines the steady-state Wigner
function of the physical cavity:

$\displaystyle W_{a,ss}(z-\alpha)$
$\displaystyle=N|\phi(\sqrt{2}z^{*})|^{2}\,e^{-2|z+\epsilon_{+}/\sqrt{2}|^{2}},$
(38)

where $\alpha\equiv\alpha_{+}/\sqrt{2}$, and $N$ is a normalization constant.
Note that, if $\lambda_{3}\equiv 0$, then the non-unitary gauge choices in Eq.
(34) satisfy $\epsilon_{+}=-\epsilon_{-}$, and so
$\epsilon_{+}-\epsilon_{-}\to 2\epsilon_{+}$, and we recover the standard
solution Eq. (30).

##  V Steady-state phase diagram of the generalized driven Kerr resonator

We now use our exact solutions in Eqs. (29) and (33) to explore the parameter
dependence of the steady state of our generalized driven-dissipative Kerr
resonator. The steady-state is largely controlled by just two dimensionless
parameters $r_{1},r_{2}$. For the usual case $\Lambda_{3}=0$ (no three photon
drive), these are:

$\displaystyle r_{1}$ $\displaystyle\equiv\frac{\lambda_{1}+\epsilon
D}{2\epsilon}=\frac{\Delta+i\frac{\kappa_{1}}{2}}{K-i\kappa_{2}}-\frac{i\Lambda_{1}}{\sqrt{\Lambda_{2}(K-i\kappa_{2})}}$
(39) $\displaystyle r_{2}$ $\displaystyle\equiv
D=\frac{2\Delta+i\kappa_{1}}{K-i\kappa_{2}}.$ (40)

The various drive amplitudes $\Lambda_{j}$ enter only through $r_{1}$; in
contrast, $r_{2}$ is a generalized detuning parameter which is independent of
drive amplitudes. With a non-zero $\Lambda_{3}$, one has
$r_{1}=(\lambda_{1}+\epsilon_{+}D)/(\epsilon_{+}-\epsilon_{-})$, $r_{2}=D$,
where $\lambda_{1},D,\epsilon_{\pm}$ are defined in Eqs. (34)-(37).

As we now show, the steady state exhibits remarkable properties whenever
system parameters are tuned to make one or both of $r_{1},r_{2}$ be non-
negative integers (see Fig. 3). At these points in parameter space, the
solution can exhibit generalized forms of photon blockade and anti-blockade,
as well as new kinds of bistability. This latter result generalizes the
previously studied cat-state bistability that occurs when
$\Lambda_{1}=\Delta=\kappa_{1}=0$ (i.e. $r_{1}=r_{2}=0$) Mirrahimi  _et al._
(2014). We stress that all of these features have clear observable signatures,
and are quantum in nature. In what follows, we focus primarily on the standard
case $\Lambda_{3}=0$. We also highlight the fact that with the addition of a
nonlinear coherent drive, the observable consequences of the photon blockade
and anti-blockade phenomena can be made even more dramatic.

Figure 3: Steady-state phase diagram for the generalized, driven-dissipative
Kerr resonator. $r_{2}$ is a dimensionless detuning parameter, whereas $r_{1}$
is a drive-dependent dimensionless parameter; both are defined in Eqs.
(39)-(40). The phase diagram indicates parameter choices that lead to unusual
steady states (as discussed in the main text).

###  V.1 Basic intuition

Recall that the steady state is determined by a single-mode pure-state
$\left|\psi_{+}\right>$ (c.f. Fig. (2)), and that further, this state is
related to a simpler state $\left|\phi\right>$ via a “non-unitary gauge
transformation” (c.f. Eq. (26)). We could always expand the transformed state
$\left|\phi\right>$ in the Fock basis as:

$\left|\phi\right>=\sum_{m=0}^{\infty}\beta_{m}\left|m\right>_{+}.$ (41)

Defining the scaled Fock state amplitudes (c.f. Eq. (41))

$c_{m}=\beta_{m}\frac{\sqrt{m!}}{\left(2\epsilon\right)^{m}}$ (42)

the ODE defining the gauge-transformed state $\left|\phi\right>$ in Eq. (29)
is equivalent to the simple recursion relation ($m\geq 0$)

$(m-r_{2})c_{m+1}=(m-r_{1})c_{m}$ (43)

with $r_{1},r_{2}$ defined in Eqs. (39),(40). The significance of
$r_{1},r_{2}$ being positive integers is now clear: in this case, there is the
possibility of the recursion relation terminating (i.e. vanishing for certain
values of $m$). This termination corresponds to a kind of quantum interference
effect, and will be at the heart of the new blockade, anti-blockade and
bistability phenomena we describe.

Note that we can directly go from the Fock state structure of
$\left|\phi\right>$ to the SB wavefunction of the desired, untransformed state
$\left|\psi_{+}\right>$. For $\Lambda_{3}=0$ the SB wavefunction of
$\left|\psi_{+}\right>$ is

$\displaystyle\psi_{+,{\rm SB}}(z)\propto e^{-\epsilon
z}\sum_{m=0}^{\infty}c_{m}\frac{\left(2\epsilon z\right)^{m}}{m!},$ (44)

with $\epsilon=i\sqrt{\lambda_{2}}$. For the more general case with non-zero
$\Lambda_{3}$, up to a displacement, we have:

$\displaystyle\psi_{+,\rm SB}(z)\propto
e^{-\epsilon_{+}z}\sum_{m=0}^{\infty}c_{m}\frac{\\{(\epsilon_{+}-\epsilon_{-})z\\}^{m}}{m!}$
(45)

where $\epsilon_{\pm}$ are defined in Eq. (34). Recall that these SB
wavefunctions directly determine the steady-state Wigner function of the
physical cavity via Eq. (15).

###  V.2 Pure unique steady states: $r_{1}=r_{2}$

The first surprising phenomena we describe is the emergence of unique pure
steady states even with nonlinearity. In general, the combination of
dissipation and nonlinearity leads us to anticipate impure cavity-$a$ steady
states. Surprisingly, there are a range of parameters where the unique steady
state of cavity $a$ is a pure coherent state (as would be expected from a
damped, linearly-driven, linear cavity). This occurs when parameters are
chosen such that $r_{1}=r_{2}$ (without either being a positive integer). In
terms of physical parameters, and for $\Lambda_{3}=0$, this requires tuning
the one and two photon drives $\Lambda_{1},\Lambda_{2}$ so that:

$-\frac{\Lambda_{1}}{\sqrt{-\Lambda_{2}(K-i\kappa_{2})}}=\frac{\Delta+i\frac{\kappa_{1}}{2}}{K-i\kappa_{2}}$
(46)

For this parameter tuning, Eq. (43) implies that all the scaled Fock state
amplitudes $c_{m}$ are identical. This in turn implies from Eq. (44) and (12)
that the state $\left|\psi_{+}\right>$ is a coherent state with amplitude
$\gamma=i\sqrt{2\Lambda_{2}/(K-i\kappa_{2})}=i\sqrt{\lambda_{2}}$. As sending
coherent states through a beamsplitter also generates coherent states at the
output, this also implies that the cavity-$a$ steady state is a simple, pure
coherent state of amplitude $\gamma/\sqrt{2}$. This follows directly from Eq.
(44) and the general expression in Eq. (15) for the steady-state cavity-$a$
Wigner function. Note that this steady-state coherent state amplitude is
consistent with the semiclassical cavity-$a$ equations of motion.

###  V.3 Higher-order photon blockade: $r_{1}=n_{0}$

Surprising effects also occur when drives and detuning are chosen so that
$r_{1}=n_{0}$, where $n_{0}$ is a non-negative integer. The recursion relation
in Eq. (43) now terminates at $m=n_{0}$: Fock state amplitudes $c_{m}$ vanish
for all $m\geq n_{0}+1$. This is an example of a generalized strong photon-
blockade phenomena: the gauge-transformed steady-state $\left|\phi\right>$ has
strictly zero probability to have more than $n_{0}$ photons. Unlike standard
photon blockade Imamoglu  _et al._ (1997), the mechanism here does not require
infinitely strong nonlinearity. Also, unlike the so-called “unconventional”
photon blockade Liew and Savona (2010); Bamba  _et al._ (2011); Lemonde  _et
al._ (2014), the blockade here is complete: there is strictly no probability
to have more than $n_{0}$ photons in the state.

Figure 4:  Generalized Photon Blockade. (a) Mean steady-state cavity-$a$
photon number as a function of single photon drive
$\delta\Lambda_{1}\equiv\Lambda_{1}-\Lambda_{1}^{(0)}$, where the offset
$\Lambda_{1}^{(0)}=(0.01-10i)\cdot K$ is determined from Eq. (47) and our
choice of system parameters. The periodic, sharp drop in photon number
corresponds to a generalized photon blockade phenomena, which occurs whenever
the parameter $r_{1}$ (c.f. Eq. (39))is a non-negative integer. Solid lines:
analytic exact solution, diamonds: master equation numerics. Photon numbers
associated with the semiclassical stationary stable amplitudes are also
plotted (dashed red lines). (b) Zoom-in of one of the blockade anti-
resonances. Loss values are $\kappa_{1}=K/100,K/20,$ and $K/10$, with more
faded green corresponding to greater loss. (c) Steady-state Wigner function,
for two choices of $\delta\Lambda_{1}$ corresponding to being either at (off)
a blockaded parameter value; black dots indicate the three semiclassical
amplitudes that exist for these parameters. For all results, $\Delta=5K$,
$\Lambda_{2}=4K$ , $\kappa_{1}=10^{-2}K$, and $\Lambda_{3}=\kappa_{2}=0$. By
using nonlinear coherent driving $\Lambda_{3}$, this blockade phenomenon can
be made sharp (i.e. there is a sharp cutoff in the photon number
distribution).

While the “gauge-transformed” state $\left|\phi\right>$ exhibits blockade,
physical phenomena is controlled by the untransformed state
$\left|\psi_{+}\right>$. Eq. (26) shows that this state is a “smeared” version
of the blockaded state. Despite this, the physical cavity $a$ steady state
still shows a pronounced suppressed photon population whenever the parameter
$r_{1}$ is tuned to an integer. This blockade-induced suppression can be
observed by considering how the steady state changes as a function of the
single photon drive amplitude $\Lambda_{1}$ (as this tunes $r_{1}$ but not
$r_{2}$). From Eq. (39), one sees that blockade occurs periodically as a
function of $\Lambda_{1}$, with the $n$th-order blockade occuring when

$\displaystyle\Lambda_{1}\equiv\Lambda_{1}^{(0)}-in\sqrt{[K-i\kappa_{2}]\Lambda_{2}},$
(47)

where
$\Lambda_{1}^{(0)}=-i(\Delta+i\kappa_{1}/2)\sqrt{\Lambda_{2}}(K-i\kappa_{2})^{-3/2}$
is a constant offset. Note that achieving a blockade requires tuning both the
phase and magnitude of the single photon drive amplitude $\Lambda_{1}$.

Fig. 4 shows representative results for $\kappa_{2}=\Lambda_{3}=0$: the
average cavity photon number shows a sharp suppression whenever $\Lambda_{1}$
is tuned to make $r_{1}$ a positive integer. Note the remarkable fact that the
width of these blockade suppressions (as a function of $\Lambda_{1}$) are much
smaller than $\kappa_{1}$. We stress that in the main plot Fig. 4(a), it is
only the single photon drive that is being tuned; all other parameters are
held fixed.

Fig. 4a also plots the photon number associated with each stable, stationary
semi-classical amplitude (obtained by solving the classical, noise-free
equation of motion). These semiclassical solutions do not exhibit any sharp
behaviour as a function of $\Lambda_{1}$. The sharp behaviour of the quantum
steady state that occurs when $r_{1}$ is tuned to a positive integer
corresponds to the quantum steady state solution suddenly switching (as a
function of $\Lambda_{1}$) from being localized near the high amplitude
classical solution to being localized near a low amplitude classical solution
(see Fig. 4(c)). The physics here is thus intimately related to physics of
quantum activiation and quantum tunneling Marthaler and Dykman (2006); Dykman
(2007), e.g. the dynamical switching between different semiclassical
solutions. For more detailed discussion of semiclassical switching behavior in
steady-states of Kerr resonators, see Bartolo  _et al._ (2016); Minganti  _et
al._ (2018). We stress that the behaviour here cannot be understood in terms
of the metapotential $M(x,y)$ often used in studies of nonlinear cavities
(see, e.g., Dykman (2012); Puri  _et al._ (2017)). The metapotential is simply
the classical Hamiltonian viewed as a function of the canonical quadratures
$x,y$. It is a completely smooth function of parameters. For $\Lambda_{1}=0$,
it has two degenerate extrema, correspondig to the two stable classical steady
states. Adding an approximately purely imaginary $\Lambda_{1}$ (as we do in
Fig. 4) tilts this metapotential, but does not break the degneracy between the
classical solutions. Hence, this does not provide any insight into why the
quantum steady state localizes around one classical amplitude versus another.

Figure 5:  Exact photon blockade using a nonlinear single-photon drive. (a)
Mean steady-state cavity-$a$ photon number as a function of single photon
drive $\Lambda_{1}$. The periodic, sharp drop in photon number corresponds to
a generalized photon blockade phenomena, which occurs whenever the parameter
$r_{1}$ (c.f. Eq. (39))is a non-negative integer. Solid lines: analytic exact
solution, diamonds: master equation numerics. Photon numbers associated with
the semiclassical stable amplitudes are also plotted (red lines). (b) Zoom-in
of one of the blockade anti-resonances. Loss values are
$\kappa_{1}=K/100,K/20,$ and $K/10$, with more faded green corresponding to
greater loss. (c) Steady-state photon statistics, for two choices of
$\Lambda_{1}$ corresponding to being either at (off) a blockaded parameter
value. For all results, $\Delta=K$, $\Lambda_{2}=0$ , $\kappa_{1}=K/100$,
$\kappa_{2}=K/1000$, and $\Lambda_{3}=K$.

Finally, we note that when $r_{1}=n_{0}$, the SB wavefunction for the dark
state $\left|\psi_{+}\right>$ (which directly determines the cavity-$a$ Wigner
function) reduces to an associated Laguerre polynomial $L^{(\alpha)}_{m}(z)$:

$\displaystyle\psi_{+,{\rm SB}}(z)\underset{r_{1}\to
n_{0}}{\propto}e^{-\epsilon z}L_{n_{0}}^{(1-D)}\left(2\epsilon z\right),$ (48)

where for $\Lambda_{3}=0$, we have $\epsilon=i\sqrt{\lambda_{2}}$.

###  V.4 Sharp photon blockade with weak nonlinearities

We now show that the generalized photon blockade phenomena is most striking
and intuitive in the case where there is no two photon driving, but only
linear and nonlinear one photon driving $\Lambda_{1},\Lambda_{3}\neq 0$.
Recall first our most general solution which includes a non-zero
$\Lambda_{2}$. In terms of the displacement parameter
$\alpha_{+}=\sqrt{2}\Lambda_{3}/\widetilde{K}^{*}$, the desired steady-state,
pure state of the $+$ mode is given by:

$\displaystyle|\psi_{+}\rangle=\hat{D}(\alpha_{+})e^{-\epsilon_{+}\hat{c}_{+}^{\dagger}}|\phi\rangle$
(49)

where $|\phi\rangle$ is the core state defined in Eq. (41). The photon
blockade phenomenon is best intuitively understood when there is nonlinear
one-photon driving, but no two-photon driving. In the limit that $\Lambda_{2}$
vanishes, $\alpha_{+}\to-\epsilon_{+}$, meaning that, after direct application
of the Baker-Campbell Hausdorff identity, the displacement transformation
partially cancels the exponential factor:

$\displaystyle\psi_{+,\rm SB}(z)\underset{\Lambda_{2}\to 0}{\sim}\phi_{\rm
SB}(z+\alpha_{+})$ (50)

where $\phi_{\rm SB}$ is the Segal-Bargmann representation of the core state
defined in Eq. (41). Therefore, we see that the limit $\Lambda_{2}\to 0$ is
physically important because it removes the smearing factors that spoil the
bare physics contained in the state $|\phi\rangle$.

Further, in the limit of vanishing $\Lambda_{2}$

$\displaystyle r_{1}\underset{\Lambda_{2}\to
0}{\sim}-\Lambda_{1}/\Lambda_{3}.$ (51)

It thus follows that a sharp photon blockade occurs in the cavity $a$ steady
state each time the ratio of nonlinear to linear one-photon driving is a
negative integer (representative results are shown in Figure 5).

In this limit, it is easy to understand the origin of the photon blockade
phenomena as the result of destructive interference between linear and
nonlinear one-photon driving. The Hamiltonian in this case is

$\displaystyle\hat{H}$
$\displaystyle=\frac{K}{2}\hat{a}^{\dagger}\hat{a}^{\dagger}\hat{a}\hat{a}-\Delta\hat{a}^{\dagger}\hat{a}$
$\displaystyle+\Lambda_{3}(\hat{n}-r_{1})\hat{a}^{\dagger}+\Lambda_{3}^{*}\hat{a}(\hat{n}-r_{1}^{*})$
(52)

When $r_{1}$ is tuned to a non-negative integer $n_{0}$, the Hamiltonian has
strictly no matrix elements connecting Fock states with photon number $n_{0}$
or less to states with photon number $n_{0}+1$ or greater. The result is that
the system becomes “trapped” in the subspace of states having $n_{0}$ or less
photons.

As this mechanism for photon blockade depends on matrix elements and not
energy detunings, it is effective even in regimes where dissipation is much
stronger than nonlinearity: while dissipation can smear out energies, it does
not smear out matrix elements, meaning that the interference preventing
excitation of the $n_{0}+1$ Fock state is robust. To see this explicitly,
consider the simplest case $n_{0}=1$, where the system gets stuck in a
subspace with at most one photon. We also consider for simplicity a system
where the only nonlinearity is the nonlinear drive $\Lambda_{3}$ (i.e.
$K=\kappa_{2}=0$), and where there is no drive detuning $\Delta$. In this
case, the steady state depends only on a single dimenionless parameter
$\Lambda\equiv\Lambda_{3}/\kappa_{1}$, and can be found using elementary
means. This state only involves the vacuum state $\left|0\right>$ and $n=1$
Fock state $\left|1\right>$ and is given by:

$\displaystyle\hat{\rho}_{\text{ss}}$
$\displaystyle=\frac{(4\Lambda^{2}+1)\left|0\right>\left<0\right|+4\Lambda^{2}\left|1\right>\left<1\right|+2i\Lambda\left(\left|1\right>\left<0\right|-h.c.\right)}{8\Lambda^{2}+1}$
(53)

In the limit of weak nonlinearity $\Lambda\rightarrow 0$, the blockade of
course still remains sharp: there is still zero probability for the state to
have $2$ or more photons, even though the one-photon probability is
$\Lambda^{2}/2\ll 1$. We stress that this mechanism is completely distinct
from the so-called unconvential photon blockade Liew and Savona (2010);
Lemonde  _et al._ (2014), which also only requires weak nonlinearities, but
which is restricted to Gaussian states, and which does not produce a sharp
blockade (e.g. there is non-zero probability to have more than one photon).

While at first glance the nonlinear one-photon drive term may seem quite
exotic, it is within reach of experiment. In Appendix A, we show how this
driving term could be realized in circuit QED using the recently developed
SNAIL architecture Sivak  _et al._ (2019). Generalized photon blockade may
have applications in quantum information science settings where nonlinearity
is a limited resource.

###  V.5 Photon anti-blockade: $r_{2}=m_{0}$

Tuning the parameter $r_{2}$ to be an integer $m_{0}$ in the recurrence
relation Eq. (43) also results in unusual behaviour of our dark states. For
zero-dissipation, $r_{2}=m_{0}$ is simply the condition for the Fock states
$n=0$ and $n=m_{0}$ of our physical $a$ cavity to be degenerate in the absence
of any driving (i.e. the detuning and Kerr terms cancel out) Bartolo  _et al._
(2016). Such resonances are analogous to multi-photon resonances that are used
to directly drive transmon qubits from the ground state to the $n$th excited
state (as a transmon can also be approximately modelled as a Kerr resonator).
Our exact solution shows that this resonance condition has strong consequences
even with dissipation and drive. When $r_{2}=m_{0}$, the only solution to the
recurrence relation has the coefficients $c_{1}$ through $c_{m_{0}}$ be
exactly zero. This implies that the gauge-transformed dark state in Eq. (41)
will have strictly zero probability to have a photon number equal to $m_{0}$
or smaller (while higher Fock states will be occupied). We call this
phenomenon a photon “anti-blockade”. As with the photon-blockade phenomena,
this will also have implications for our physical cavity, via Eq. (26).

Figure 6:  Photon anti-blockade. Average cavity-$a$ steady-state photon number
as a function of drive detuning $\Delta$, with drive amplitudes fixed at
$\Lambda_{1}=\Lambda_{2}=K/2$, $\Lambda_{3}=0$. Resonances here correspond to
having tuned the parameter $r_{2}$ (c.f. Eq. (40)) to be near a non-negative
integer. Other parameters are $\kappa_{1}=0.01K$, $\kappa_{2}=0$.

For any non-zero amount of dissipation, it is clear from Eq. (40) that we can
never have $r_{2}$ exactly be a positive integer (as
$\kappa_{1},\kappa_{2}\geq 0$). This remains true even in the presence of a
nonlinear coherent drive, where $r_{2}=D$, with $D$ given by Eq. (35).
Nonetheless, for weak dissipation (namely $\kappa_{1}\ll\Delta$,
$\kappa_{2}\ll K$), one can still tune $r_{2}$ to be extremely close to an
integer. In this regime, one still has strong signatures of the anti-blockade
behaviour. For the physical cavity $a$, this translates into a kind of
resonant enhancement of photon number and skewed photon number statistics.
Representative behaviour is shown in Fig. 6. Note that this resonance
phenomenon was observed in Ref. Bartolo  _et al._ (2016), though connections
to photon statistics and the properties of the analytic steady-state solution
were not discussed.

###  V.6 Generalized bistability: $(r_{1},r_{2})=(n_{1},n_{2})$

Having understood photon blockade and anti-blockade phenomena, the natural
remaining case is when both these phenomena coexist. This occurs when
parameters are chosen so that $(r_{1},r_{2})=(n_{1},n_{2})$, where $n_{2}\geq
n_{1}$ are both non-negative integers. Eq. (43) then yields both a photon-
blockaded solution, and a distinct anti-blockaded solution. These correspond
to two distinct dark states of the $+$ mode, described respectively by SB
wavefunctions:

$\displaystyle\psi_{1,{\rm SB}}(z)$ $\displaystyle=\frac{e^{-\epsilon
z}}{N_{1}^{1/2}}\sum_{m=0}^{r_{1}}\leavevmode\nobreak\ \leavevmode\nobreak\
\frac{c_{m}(2\epsilon z)^{m}}{m!},$ (54) $\displaystyle\psi_{2,{\rm SB}}(z)$
$\displaystyle=\frac{e^{-\epsilon
z}}{N_{2}^{1/2}}\sum_{m=r_{2}+1}^{\infty}\frac{c_{m}(2\epsilon z)^{m}}{m!}.$
(55)

Any linear combination of these solutions is also a dark steady-state. We
refer to this situation as “quantum bistability”: the extended, two-cavity
cascaded system in Fig. 1(b) has an infinite number of steady states,
corresponding to any superposition state of the form:

$\left|\tilde{\psi}[a_{1},a_{2}]\right>=a_{1}\left|\psi_{1}\right>_{+}\left|0\right>_{-}+a_{2}\left|\psi_{2}\right>_{+}\left|0\right>_{-}$
(56)

This steady-state structure is conventionally referred to as a (two-
dimensional) decoherence-free subspace. This also implies multi-stability for
the physical $a$ cavity, which exhibits a two-parameter continuous family of
steady states,

$\hat{\rho}_{a,{\rm
ss}}=\textrm{tr}_{b}\left[\left|\tilde{\psi}[a_{1},a_{2}]\right>\left<\tilde{\psi}[a_{1},a_{2}]\right|\right].$
(57)

The upshot is that the generalized driven-dissipative Kerr cavity has a
multitude of distinct parameter points that yield multi-stability, despite any
obvious symmetry.

Unfortunately, we have the same issue as with the anti-blockade phenomena:
non-zero dissipation makes it impossible to exactly tune to bistable parameter
values except for the case $n_{1}=n_{2}=0$. This is because the constraint of
having one or both of $\kappa_{1},\kappa_{2}$ be positive implies $r_{2}$
cannot be exactly equal to a positive integer (c.f. Eq. (40)). The only
exactly-achievable bistable point is the case $n_{1}=n_{2}=0$, which can be
reached if $\kappa_{1}=0,\kappa_{2}>0$. This parameter point corresponds to
the well-studied cat-state bistability in a two-photon driven Kerr resonator
Wolinsky and Carmichael (1988).

Despite these caveats, the new bistable points are physically relevant: for
weak dissipation, one can come arbitrarily close to them in parameter space,
with striking observable consequences for the steady state. We explore this
further in the next section. We also discuss in Appendix B how one can exactly
achieve the physics of these new bistable points using a non-cascaded version
of the two cavity setup depicted in Fig. 1(b).

###  V.7 Simultaneous/coexisting blockade and anti-blockade

What if $r_{1},r_{2}$ are both non-negative integers, and $r_{2}<r_{1}$? In
this case, neither the photon-blockaded nor the photon-resonant solution is
permitted. Instead, a medium-photon number solution exists, and serves as the
unique dark state. For $\Lambda_{3}=0$, we have

$\displaystyle\psi_{+,\rm SB}(z)=\frac{e^{-\epsilon
z}}{N^{1/2}}\sum_{r_{2}+1}^{r_{1}}c_{m}\frac{(2\epsilon z)^{m}}{m!}$ (58)

with $\epsilon\equiv i\sqrt{\lambda_{2}}$. Without the exponential prefactor,
this state would exhibit both photon blockade and anti-blockade (i.e. its
photon number distribution would be cut-off at small and large photon
numbers).

##  VI Consequences of new quantum bistable points

As discussed in the previous section, there are an infinite number of points
in parameter space where our generalized driven-dissipative Kerr resonator is
almost quantum bistable (c.f. Fig. 3). With non-zero one-photon loss, one
cannot exactly achieve the required parameter tuning for bistability, but one
can come arbitrarily close to a given bistable parameter point. In this
section, we explore the physical consequences of this near-bistability. We
show that there is an extremely strong sensitivity to small parameter changes
when one is in this near-bistable regime, and that the unique steady state can
be understood as “picking-out” a unique state from the bistable manifold in
Eq. (57).

Suppose we chose parameters that result in $(r_{1},r_{2})$ being close to
integers $(n_{1},n_{2})$:

$r_{1}=n_{1}+\delta r_{1},\hskip 28.45274ptr_{2}=n_{2}+\delta r_{2},$ (59)

These small deviations kill the bistability. However, for small $\delta r_{j}$
the resulting pure steady state of the $+$ mode is a particular linear
combination of the states $\phi_{j}(z)$ that span the bistable manifold at
$\delta r_{j}=0$. Moreover, the precise form of this combination is extremely
sensitive to parameter variations.

For example, consider the simple case where the unperturbed recursion
parameters are $(r_{1},r_{2})=(n,n)$. In this case, the recursion relation Eq.
(43) simplifies to

$\displaystyle c_{m+1}=\frac{(m-n)-\delta r_{1}}{(m-n)-\delta r_{2}}c_{m}.$
(60)

In the regime that $\delta r_{1},\delta r_{2}\ll 1$, we can see that the ratio
$c_{m+1}/c_{m}$ is essentially $1$, except for the ratio $c_{n+1}/c_{n}=\delta
r_{1}/\delta r_{2}$. Therefore, as $\delta r_{1},\delta r_{2}\to 0$, the
unique steady state solution (i.e. solution to the recursion relation) has the
limiting form

$\displaystyle\psi_{\text{+,SB}}(z)\underset{\delta r_{1},\delta r_{2}\to
0}{\sim}\psi_{\text{1,SB}}(z)+\frac{\delta r_{1}}{\delta
r_{2}}\psi_{\text{2,SB}}(z)$ (61)

as a superposition of the bistable solutions given in Eqs. (54),(55). Note
that in writing this equation, we must pick the overall phase of $\psi_{2,{\rm
SB}}$ such that the ratio between $c_{n+1}$ (appearing in
$\psi_{\text{2,SB}}$) and $c_{n}$ (appearing in $\psi_{\text{1,SB}}$) is
precisely $\delta r_{1}/\delta r_{2}$.

As a result, the unique steady state Wigner function of the physical $a$
cavity will be:

$\displaystyle
W_{a,ss}(z)\simeq\frac{e^{-2\left|z\right|^{2}}}{N}\bigg{|}\psi_{1,{\rm
SB}}(\sqrt{2}z^{*})+\frac{\delta r_{1}}{\delta r_{2}}\,\cdot\psi_{2,{\rm
SB}}(\sqrt{2}z^{*})\bigg{|}^{2}.$ (62)

This equation is the crucial result of this subsection: for parameters that
bring us close to a quantum bistable point, it provides a simple way to
understand the system’s steady state and its extreme sensitivity to small
parameter changes.

Figure 7: Extreme parameter sensitivity near a quantum-bistable point. Bottom
row: Wigner function for the purification $\psi_{+,\rm SB}$ of the Kerr-cavity
steady state, for various parameter choices that are close to the
$(r_{1},r_{2})=(0,0)$ quantum bistable point. For all plots
$\Lambda_{3}=\Delta=\kappa_{2}=0$, $\kappa_{1}=10^{-2}K$, and
$\Lambda_{2}=4K$, and $\Lambda_{1}$ increases from left to right: (a)
$\Lambda_{1}=0$ ($Q=0$), (b), $\Lambda_{1}=0.01K$ ($Q=1$), (c)
$\Lambda_{1}=0.02K$ ($Q=2$), and (d) $\Lambda_{1}=0.1K$ ($Q=10$). The small
value of $\kappa_{1}$ and $\Lambda_{1}$ imply that one is not exactly at the
bistable point; the $Q$ parameter then controls the form of the unique steady
state, c.f. Eq. (68). By tuning the single-photon drive amplitude
$\Lambda_{1}$, one can pick out a particular superposition in the "bistable"
manifold by varying $Q$. Top row: corresponding metapotential $M(z)$ for the
same parameter choices. The metapotential is essentially unchanged for this
range of $\Lambda_{1}$, showing that it cannot be used to understand the large
changes in the quantum steady state.

###  VI.1 Cat-state bistability: $(r_{1},r_{2})=(0,0)$

The simplest bistable point is where $r_{1}=r_{2}=0$. From Eqs. (39),(40), we
see that this requires there to be no single photon drive or loss, nor any
detuning: $\Lambda_{1}=\Delta=\kappa_{1}=0$. This corresponds to the well-
known quantum bistability that occurs in a two-photon driven Kerr resonator
Goto (2016); Puri  _et al._ (2017); Grimm  _et al._ (2019), a system where
photon number parity is conserved. The two distinct solutions to the
recurrence relation in Eq. (43) are $c_{j}=\delta_{j,0}$ and
$c_{j}=1-\delta_{j,0}$ (c.f. Eqs. (54)-(55)). This corresponds to two distinct
dark states for the $+$ mode, with SB wavefunctions

$\displaystyle\psi_{1,\textrm{SB}}(z)$ $\displaystyle=e^{-\epsilon z}$ (63)
$\displaystyle\psi_{2,\textrm{SB}}(z)$ $\displaystyle=e^{\epsilon
z}-e^{-\epsilon z}$ (64)

$\psi_{1,\textrm{SB}}(z)$ corresponds to a coherent state with amplitude
$\epsilon\equiv i\sqrt{\lambda_{2}}$, whereas $\psi_{2,\textrm{SB}}(z)$
corresponds to an odd cat state (odd superposition of coherent states with
amplitude $\epsilon$). Note that we have picked the global phase of
$\psi_{2,\rm SB}$ to be compatible with Eq. (62).

We thus have a direct connection between this parity-based bistability and the
photon blockade and anti-blockade discussed above: bistability corresponds to
both these phenomena occurring simultaneously. As always, any amount of
single-photon loss will kill the bistability and yield a unique steady state
(though relaxation to this state could be extremely slow). Our approach gives
a simple way to understand the unique steady state when there is weak single
photon loss, and possibly other weak perturbations (such as single photon
driving and/or a detuning). These imperfections cause a shift in the recursion
parameters away from $\delta r_{1}=\delta r_{2}=0$:

$\displaystyle\delta r_{2}$
$\displaystyle=\frac{2\Delta+i\kappa_{1}}{K-i\kappa_{2}}$ (65)
$\displaystyle\delta r_{1}$ $\displaystyle=\frac{\delta
r_{2}}{2}-\frac{i\Lambda_{1}}{\sqrt{\Lambda_{2}(K-i\kappa_{2})}}$ (66)

For small imperfections, we can then use Eq. (62) to give us the steady-state
SB wavefunction:

$\psi_{+,\textrm{SB}}(z)=N\left[(1+Q)e^{-\epsilon z}+(1-Q)e^{\epsilon
z})\right],$ (67)

where $N$ is a normalization constant, and

$Q=\sqrt{\frac{K-i\kappa_{2}}{\Lambda_{2}/4}}\frac{i\Lambda_{1}}{2\Delta+i\kappa_{1}}.$
(68)

Eq. (67) directly gives us the Wigner function of the unique steady state via
Eq. (15). Each term in Eq. (67) on its own corresponds to a simple coherent
state (amplitudes $\pm\epsilon\equiv\pm i\sqrt{\lambda_{2}}$). This equation
also reveals something surprising: the localization of the steady state in
phase space is a non-monotonic function of $\Lambda_{1}$. The state is
delocalized both for $\Lambda_{1}=0$, and for $\Lambda_{1}$ large enough to
make $Q\gg 1$. Representative results are shown in Figure 7; We plot the
semiclassical metapotential for in this figure for each parameter choice; it
shows almost no changes, indicating that it cannot be used to understand the
strong parameter-sensitivity of the quantum steady state.

###  VI.2 Quantum bistability with a single photon drive:
$(r_{1},r_{2})=(n,n)$

A more surprising regime of near bistability is when the recursion parameters
are both tuned to be close to the same positive integer, i.e.
$(r_{1},r_{2})\simeq(n,n)$. As discussed, for an exact tuning to this point,
the expanded system exhibits quantum bistability. There are two orthogonal
solutions to the recurrence relations, given by
$c_{j}=\sum_{k=1}^{n}\delta_{j,k}$ and $c_{j}=\sum_{k=1}^{n}(1-\delta_{j,k})$
(c.f. Eq. (54-55)). These in turn correspond to two distinct $+$-mode states

$\displaystyle\psi_{1,\textrm{SB}}(z)$ $\displaystyle=N_{1}e^{-\epsilon
z}\,\Gamma(n+1,2\epsilon z)$ (69) $\displaystyle\psi_{2,\textrm{SB}}(z)$
$\displaystyle=N_{2}e^{-\epsilon z}\Bigg{(}1-\frac{\Gamma(n+1,2\epsilon
z)}{\Gamma(n+1)}\Bigg{)}$ (70)

where $\Gamma(r,z)\equiv\int_{z}^{\infty}t^{r-1}e^{-t}dt$ is the incomplete
Gamma function.

In the absence of any loss, tuning $r_{1}=r_{2}=n$ requires a detuning
$\Delta=n/2K$ and a single photon drive
$\Lambda_{1}=-i(n/2)\sqrt{\Lambda_{2}K}$. If we now include single photon loss
(but keep $\kappa_{2}=0$), and also shift $\Lambda_{1}$ slightly from the
above value, the recurrence parameters are slightly shifted as well:

$\displaystyle r_{1}$ $\displaystyle=$ $\displaystyle
n+\frac{i\kappa_{1}}{2K}\equiv n+\delta r_{1}$ (71) $\displaystyle r_{2}$
$\displaystyle=$ $\displaystyle
n+\frac{i\kappa_{1}}{4K}-i\frac{\delta\Lambda_{1}}{\sqrt{K\Lambda_{2}}}\equiv
n+\delta r_{2}$ (72)

Hence, via Eq. (62), by slightly varying the one photon drive amplitude, one
can pick out completely different linear combinations of the two different
bistable states as the single unique steady state. This leads to an extreme
sensitivity of the final state to small changes in $\Lambda_{1}$. Note that by
picking parameters so that $\delta r_{1}=\delta r_{2}$, the steady state
becomes a coherent state with amplitude $\gamma=\sqrt{-\Lambda_{2}/K}$,
whereas if $\delta r_{1}=0$, it has a bimodal form.

###  VI.3 Metastability due to proximal quantum bistability

Tuning parameters to be close to a quantum bistable point also has
consequences for dynamics. The characteristic decay rates of the system
correspond to the non-zero eigenvalues of the Liouvillian $\mathcal{L}_{0}$
(c.f. Eq. (2)). We find that tuning to a regime of near-bistability gives rise
to an extremely slow population-decay mode, and also a clear dissipative gap
separating the rate of this slow-mode from other decay modes. Formally, if we
let $\gamma_{j}$ denote the decay modes of the Liouvillian (i.e. negative real
parts of the eigenvalues of $\mathcal{L}_{0}$), and order rates such that
$\gamma_{1}\leq\gamma_{2}\leq....$, then in near-bistable regimes:

$\gamma_{1}\ll\kappa_{1},\,\,\,\,\,\gamma_{2}\gg\gamma_{1}$ (73)

Note that this hierarchy of dissipative rates has already been described for
the more familiar $(r_{1},r_{2})=(0,0)$ “cat-state” bistable point Puri  _et
al._ (2017); we show that this is also true for our new bistable points. An
exact description of this dynamical behaviour is outside the scope of the CQA
method. It can however be studied numerically. Representative behavior of a
driven Kerr cavity whose parameters are close to either the
$(r_{1},r_{2})=(2,2)$ or $(4,4)$ bistable points are shown in Fig. 8(a).

For near-bistable parameters, the CQA approach provides insight into the
nature of the slow decay mode of $\mathcal{L}_{0}$. As one might expect, this
mode corresponds to slow relaxation within the bistable manifold of states.
For more general works on metastability in open quantum systems, see
Macieszczak  _et al._ (2016); Kessler  _et al._ (2012). To make this precise,
recall that if one tuned exactly to a bistable parameter point, cavity-$a$ has
a continuous three-parameter family of possible steady states corresponding to
Eq. (57) (and incoherent mixtures of these states). Density matrices in this
bistable manifold lie in the span of the four operators ($i,j=1,2$):

$\displaystyle\hat{M}_{ij}$
$\displaystyle=\text{tr}_{b}\left[\left(\left|\psi_{i}\right>\left<\psi_{j}\right|\right)_{+}\left(\left|0\right>\left<0\right|\right)_{-}\right]$
(74)

By Appendix C, these operators have Wigner transforms

$\displaystyle W_{ij}(z)$
$\displaystyle=N\psi_{i,\text{SB}}(\sqrt{2}z^{*})\psi_{j,\text{SB}}^{*}(\sqrt{2}z^{*})e^{-2|z|^{2}}$
(75)

with $N$ a normalization constant.

Figure 8:  Slow dynamics near generalized bistable regimes. (a) Solid line:
ratio of the two smallest relaxation rates (i.e. dissipative rates of the
system Liouvillian $\mathcal{L}_{0}$), as a function of $\kappa_{1}$. Dashed
line: $\kappa_{1}/\gamma_{1}$. $\kappa_{1}\rightarrow 0$ corresponds to being
at a bistable parameter point $(r_{1},r_{2})=(n,n)$, either $n=2$ (light
green) or $n=4$ (dark green). Parameters are $\Lambda_{2}=6K$,
$\Lambda_{3}=\kappa_{2}=0$, $\Delta=nK/2$, and
$\Lambda_{1}=-in\sqrt{\Lambda_{2}}/2$. One sees that the slow rate
$\gamma_{1}$ is much slower than $\kappa_{1}$, and that there is a pronounced
dissipative gap. (b) Solid line: The measure $1-P$ (c.f. Eq. (76)) of how
closely the slowest system decay mode (with rate $\gamma_{1}$) corresponds to
dynamics in the bistable manifold. Dashed lines: same, but measuring how
closely this mode is described by coherent states centered at the
semiclassical stable amplitudes. One clearly sees that the bistable manifold
gives a far better description. Same parameters as in (a).

The slow mode (rate $\gamma_{1}$) has an associated right eigenvector
$\hat{M}_{\rm slow}$, i.e. $\mathcal{L}_{0}\hat{M}_{\rm
slow}=-\gamma_{1}\hat{M}_{\rm slow}$. If the slow dynamics is entirely in the
bistable manifold, then $\hat{M}_{\rm slow}$ should lie completely within the
span of the $\hat{M}_{ij}$. To see whether this is the case, we pick
parameters for near-bistability, and numerically calculate the Hilbert-Schmidt
norm $P$ of the projection of the slow mode onto the bistable subspace:

$\displaystyle
P\equiv\sum_{i,j=1,2}|\text{Tr}[\widetilde{M}_{ij}^{{\dagger}}\hat{M}_{\text{slow}}]|^{2}$
(76)

Here $\widetilde{M}_{ij}$ is an orthonormal basis for the span of
$\hat{M}_{ij}$ (obtained via the standard Gram-Schmidt process). As $0\leq
P\leq 1$, the quantity $1-P$ measures how much of the slow mode’s dynamics
lies outside the bistable manifold.

Representative results for $1-P$ are shown in Fig. 8(b). One sees that for
small $\kappa_{1}$ (i.e. when one is close to the bistable point), the slow
mode is almost entirely described by the bistable state manifold. For
comparison, we have also tried to describe the slow mode in terms of simple
coherent states centered at the expected classical bistable steady-state
amplitudes. This involves taking

$\displaystyle\hat{M}^{ij}_{\text{cat}}$
$\displaystyle=|\alpha_{i}\rangle\langle\alpha_{j}|$ (77)

with $\alpha_{j}$ the classical amplitudes, determined by (with $K\equiv 1$):

$\displaystyle\Delta(\alpha_{j}-i\sqrt{\Lambda_{2}})-\Lambda_{2}\alpha^{*}_{j}-\alpha_{j}|\alpha_{j}|^{2}\equiv
0,\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\ j=1,2$ (78)

One sees from Fig. 8(b) that this coherent-state description does a far poorer
job of describing the dynamical slow mode, compared to the states from the
bistable manifold. While metastability in two-photon driven Kerr resonators
was discussed in Minganti  _et al._ (2018), its connection to the existence of
nearby, novel quantum bistable points (i.e. generalized cat-state regimes) in
the resonator’s phase diagram has not been previously investigated.

##  VII Parity-conserving dynamics: true quantum bistability

We now focus on a special case that has received considerable recent attention
Mirrahimi  _et al._ (2014); Goto (2016); Puri  _et al._ (2017); Grimm  _et
al._ (2019): a system where $\kappa_{1}=\Lambda_{1}=0$ in Eq. (2), implying
that the full dynamics conserves photon number parity. This in turn implies
that there are at least two distinct steady states, and opens the possibility
of true quantum bistability. Note that a comprehensive discussion of generic
Lindblad master equations with multiple steady states is provided in Ref.
Albert and Jiang (2014). Our exact-solution CQA method provides several
insights into this regime. Among other things, it allows one to understand why
adding a drive-detuning destroys quantum bistability despite parity still
being conserved, something that is not possible with $P$-function methods,
which give a unique solution Hach III and Gerry (1994); Gilles  _et al._
(1994). In addition, the CQA method also gives a succinct analytical
expression that controls which unique steady state is selected from the
bistable manifold when quantum bistability is broken.

We start by revisiting the CQA method of Sec. III for systems described by Eq.
(2) with $\kappa_{1}=\Lambda_{1}=0$. The corresponding cascaded two-cavity
system is described by Eqs. (9) and (10). The first step as always is to
insist that we have a state that is dark with respect to the cascaded
dissipators. For $\kappa_{1}=0$, we only have the a two-photon loss
dissipator, given by

$\mathcal{D}[\hat{a}^{2}-\hat{b}^{2}]=\mathcal{D}[2\hat{c}_{+}\hat{c}_{-}]$
(79)

where again the collective $\hat{c}_{\pm}$ modes are defined in Eq. (7). There
are now two distinct possibilities for a non-trivial dark state: either the
$\hat{c}_{-}$ mode is forced to be in vacuum (with the $+$ mode occupied), or
the $\hat{c}_{+}$ mode is forced to be in vacuum (with the $-$ mode occupied).
The first option is the same as what we did for $\kappa_{1}=0$; the second
option is a new possibility enabled by the lack of one photon loss.

It follows that the most general 2-cavity dark state has the form:

$\left|\psi_{\rm
dk}\right>=\alpha_{+}\left|\psi\right>_{+}\left|0\right>_{-}+\alpha_{-}\left|0\right>_{+}\left|\theta\right>_{-}$
(80)

This structure is a direct consequence of parity conservation, which
guarantees the existence of at least two orthogonal steady states (one even
parity, one odd parity). This structure also implies that the general argument
in Sec. III ensuring a positive cavity-$a$ steady-state Wigner function no
longer holds, as $\left|\psi_{\rm dk}\right>$ can have both $+$ and $-$ modes
occupied.

###  VII.1 Zero detuning: quantum bistability

Consider first the case $\Delta=0$, meaning that we have a resonantly-driven
Kerr parametric oscillator subject to two photon loss. Without dissipation,
this system has degenerate coherent state eigenstates Goto (2016); Puri  _et
al._ (2017). Including two-photon loss, the dissipative system exhibits true
quantum bistability: the steady-state manifold corresponds to a two-
dimensional decoherence-free subspace, in the language of Albert (2018). We
show how this structure emerges via the CQA method.

The first step of the CQA method is to identify possible pure dark states of
the collective cascaded-systems dissipators; in our case, this is Eq. (80). To
have this state be a steady state, it must also be an eigenstate of the
cascaded Hamiltonian (c.f. Eq. (10)). For $\Delta=0$, this leads to the
equation:

$\displaystyle\hat{c}_{-}^{\dagger}\hat{c}_{+}^{\dagger}\Bigg{(}\alpha_{+}\left[\hat{c}_{+}^{2}+\lambda_{2}\right]\left|\psi\right>_{+}\left|0\right>_{-}$
$\displaystyle+\,\alpha_{-}\left[\hat{c}_{-}^{2}+\lambda_{2}\right]\left|0\right>_{+}\left|\theta\right>_{-}\Bigg{)}=E\left|\psi_{\rm
dk}\right>$ (81)

with $\lambda_{2}=2\Lambda_{2}/(K-i\kappa_{2})$. Since $\hat{H}_{\text{casc}}$
always adds an excitation to both modes $\hat{c}_{-}$ and $\hat{c}_{+}$, the
only possible energy eigenvalue is $E=0$. The equation then decouples into
separate equations for $\left|\psi\right>_{+},\left|\theta\right>_{-}$ which
are easily solved. Crucially, each of these equations admits two possible
solutions.

As a result, one finds that the most general dark state solution can be
written in terms of coherent states as:

$\displaystyle\left|\psi_{\rm dk}\right>$
$\displaystyle=\sum_{\pm}\left(\mu_{\pm}\left|\pm\epsilon\right>_{+}\left|0\right>_{-}+\nu_{\pm}\left|0\right>_{+}\left|\pm\epsilon\right>_{-}\right)$
(82)

where the coherent state amplitude $\epsilon=i\sqrt{\lambda_{2}}$, as in the
previous section. We see that the cascaded two-cavity system has a four
dimensional subspace of possible steady-state, dark states.

The last step is to determine the corresponding steady state structure of the
physical $a$ cavity. As discussed in Sec. III, this effectively corresponds to
taking a given two-cavity state, sending it through a 50-50 beamsplitter, and
then discarding one of the outputs. This procedure is easy to carry out on the
general state in Eq. (82), as coherent states transform in a simple manner
under a beamsplitter operation. In the basis of the physical $a$ cavity and
auxiliary $b$ cavity, our general dark state has the form:

$\displaystyle\left|\psi_{\rm dk}\right>$
$\displaystyle=\sum_{\pm}\Big{(}\mu_{\pm}\left|\pm\tilde{\epsilon}\right>_{a}\left|\pm\tilde{\epsilon}\right>_{b}$
$\displaystyle\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\
+\nu_{\pm}\left|\pm\tilde{\epsilon}\right>_{a}\left|\mp\tilde{\epsilon}\right>_{b}\Big{)}$
(83)
$\displaystyle=\Big{(}\mu_{+}\left|\tilde{\epsilon}\right>_{a}+\nu_{-}\left|\tilde{\epsilon}\right>_{a}\Big{)}\left|\tilde{\epsilon}\right>_{b}$
$\displaystyle\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\
+\Big{(}\mu_{-}\left|-\tilde{\epsilon}\right>_{a}+\nu_{+}\left|+\tilde{\epsilon}\right>_{a}\Big{)}\left|-\tilde{\epsilon}\right>_{b},$
(84)

with $\tilde{\epsilon}\equiv\epsilon/\sqrt{2}$. As there is in general
entanglement between the physical $a$ cavity and the auxiliary $b$ cavity, one
in general is left with an impure state for cavity $a$. However, pure
cavity-$a$ steady states are indeed possible; consider for example the case
where $\mu_{-}=\nu_{+}=0$.

The upshot is that we have a steady state manifold for cavity $a$ that is two
dimensional, and spanned by the states $\left|\pm\tilde{\epsilon}\right>_{a}$
(in agreement with previous work Mirrahimi  _et al._ (2014)). In simple terms,
the steady-state manifold corresponds to a quantum bit, i.e. a full single-
qubit Bloch sphere Albert and Jiang (2014). This is what we mean by the system
exhibiting quantum bistability.

###  VII.2 Non-zero detuning: classical bistability

####  VII.2.1 Loss of quantum bistability

We next consider the case of a adding a non-zero detuning $\Delta$. As has
been discussed previously Mirrahimi  _et al._ (2014), this causes the steady-
state manifold to transition from being a two-dimensional decoherence free
subspace (i.e. quantum bistability) to having the structure of a simple
classical bit. Formally, it corrsponds to an orthogonal direct sum of two one-
dimensional noiseless subsystems (one for each parity sector). We will find
that, in contrast to $P$-function methods Minganti  _et al._ (2016), CQA is
able to analytically detect this transition, and gives closed-form expressions
for each of the direct-summands in the steady-state manifold. Solving the
system again using the CQA method, the requirement of having our general dark
state in Eq. (80) be a energy eigenstate of the cascaded Hamiltonian leads to
the equations

$\displaystyle\hat{c}_{-}^{\dagger}\left(\hat{c}_{+}^{\dagger}\hat{c}_{+}^{2}-D\hat{c}_{+}+\lambda_{2}\hat{c}_{+}^{\dagger}\right)|\psi\rangle_{+}$
$\displaystyle=0,$ (85)
$\displaystyle\hat{c}_{+}^{\dagger}\left(\hat{c}_{-}^{\dagger}\hat{c}_{-}^{2}-D\hat{c}_{-}+\lambda_{2}\hat{c}_{-}^{\dagger}\right)|\theta\rangle_{-}$
$\displaystyle=0.$ (86)

where $D=2\Delta/(K-i\kappa_{2})$. As before, the equations determining
$\left|\psi\right>_{+}$ and $\left|\theta\right>_{-}$ are identical
(reflecting parity conservation). The equation in each case can be solved by
using a SB representation for the state, and turning the operator equations
into differential equations. We get the same ODE in each case:

$\displaystyle\Big{(}z\frac{\partial^{2}}{\partial
z^{2}}-D\frac{\partial}{\partial z}+\lambda_{2}z\Big{)}\psi_{\rm SB}(z)=0$
(87)

With the same equation for $\theta_{\rm SB}$, and with $\lambda_{2},D$ having
the same definitions as earlier in the main text.

At the qualitative level, one can see how true quantum bistability is lost in
the presence of detuning: for zero detuning $D\equiv 0$, the ODE above has no
singular points, and thus the standard existence theorem (§12.22 in Ince
(1956)) guarantees two independent, analytic solutions. As discussed earlier,
this leads to quantum bistability for the physical mode $a$. However, the term
$\propto D\partial_{z}$ introduces a singular point into the ODE at $z=0$, and
the existence of two dark steady states is no longer guaranteed. Indeed, the
singular point at $z=0$ produces a branch-cut discontinuity in one of the
solutions. Generically, only one analytic solution survives:

$\displaystyle\psi_{\text{SB}}(z)$
$\displaystyle=\frac{1}{N^{1/2}}\,_{0}F_{1}(1/2-D/2;-\lambda_{2}z^{2}/4),$
(88)

where $N$ is a normalization constant.

As we will see, this two-fold reduction in the number of dark steady-states
has dramatic consequences for the bistability of the physical mode $a$. As
there is a unique choice for both $\left|\psi\right>_{+}$ and
$\left|\theta\right>_{-}$, the most general dark state has the form of Eq.
(80), and corresponds to a two-dimensional subspace. In what follows, it will
be useful to write this general dark state as

$\left|\psi_{\rm dk}\right>=\mu_{\rm e}\left|\Phi_{\rm e}\right>+\mu_{\rm
o}\left|\Phi_{\rm o}\right>$ (89)

with

$\displaystyle\left|\Phi_{\rm e/o}\right>$ $\displaystyle=\frac{1}{\sqrt{2\pm
2N^{-1}}}\left(\left|\psi\right>_{+}\left|0\right>_{-}\pm\left|0\right>_{+}\left|\psi\right>_{-}\right).$
(90)

We now trace-out the auxiliary $b$ cavity. Note that the pure dark states
above span a subspace of dimension 2. Incoherent mixtures in this subspace are
also stationary states; hence the cascaded 2 cavity problem has a steady-state
manifold corresponding to a Bloch sphere. We imagine starting with an abitrary
mixed state in this subspace (described by a 2 cavity density matrix), and
then tracing out cavity $a$ to determine the corresponding cavity $a$ state.
Understanding the full range of cavity $a$ states produced here determines the
steady-state manifold of cavity $a$.

This procedure leads us to consider four linearly-independent cavity-$a$
operators (that determine the cavity $a$ density matrix after tracing out
cavity $b$):

$\displaystyle\hat{M}^{++}_{a,ss}$
$\displaystyle\equiv\text{tr}_{b}[|\psi\rangle_{+}|0\rangle_{-}\langle\psi|_{+}\langle
0|_{-}]$ $\displaystyle\hat{M}^{--}_{a,ss}$
$\displaystyle\equiv\text{tr}_{b}[|0\rangle_{+}|\psi\rangle_{-}\langle
0|_{+}\langle\psi|_{-}]$ $\displaystyle\hat{M}^{+-}_{a,ss}$
$\displaystyle\equiv\text{tr}_{b}[|\psi\rangle_{+}|0\rangle_{-}\langle
0|_{+}\langle\psi|_{-}]$ $\displaystyle\hat{M}^{-+}_{a,ss}$
$\displaystyle\equiv(\hat{M}^{+-}_{a,ss})^{\dagger}$ (91)

To understand the structure of these operators, we consider their
corresponding $Q$-functions (easily obtainable using the SB representation):

$\displaystyle Q^{\pm\pm}_{a,ss}(z)$ $\displaystyle=\int d^{2}u\,\psi_{\rm
SB}^{*}\Big{(}\frac{z\pm u}{\sqrt{2}}\Big{)}\psi_{\rm SB}\Big{(}\frac{z\pm
u}{\sqrt{2}}\Big{)}$ (92) $\displaystyle Q^{\pm\mp}_{a,ss}(z)$
$\displaystyle=\int d^{2}u\,\psi_{\rm SB}^{*}\Big{(}\frac{z\pm
u}{\sqrt{2}}\Big{)}\psi_{\rm SB}\Big{(}\frac{z\mp u}{\sqrt{2}}\Big{)}$ (93)

We obtain an important result: these four operators are not all independent.
Because of the symmetry of each integral under the mapping $u\to-u$, we have

$\displaystyle\hat{M}^{++}_{a,ss}=\hat{M}^{--}_{a,ss},\leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\
\hat{M}^{+-}_{a,ss}=\hat{M}^{-+}_{a,ss}=(\hat{M}^{+-}_{a,ss})^{\dagger}$ (94)

These equalities imply a loss of information in tracing out cavity-$b$, and
result in the cavity-$a$ steady state manifold being simply two-dimensional.
It is spanned by the quantities

$\displaystyle\hat{\rho}_{a,ss}^{+}$ $\displaystyle\equiv\hat{M}_{a,ss}^{++}$
$\displaystyle\hat{\rho}_{a,ss}^{-}$ $\displaystyle\equiv\hat{M}_{a,ss}^{+-}.$
(95)

with $\hat{\rho}_{a,ss}^{\pm}$ both Hermitian. We now have enough information
to calculate each steady state exactly: since the steady-state manifold is
two-dimensional, and since parity is a conserved quantity, every density
matrix in the manifold must then be an impure mixture of the form

$\displaystyle\hat{\rho}_{a,ss}$
$\displaystyle=p\hat{\rho}_{e}+(1-p)\hat{\rho}_{o}$ (96)

where the extremal states $\hat{\rho}_{e/o}$ are uniquely characterized by the
property of having definite photon number parity (even and odd respectively).

Thus, in summary, in this case there is a distinct steady state in both the
even and odd photon number sectors; any mixture of these is also a possible
steady state. The steady state manifold is indexed by just a single number
$0\leq p\leq 1$, which simply corresponds to the dynamically-conserved
probability of having an even photon number parity. In simpler terms, the
cavity-$a$ steady-state manifold corresponds to a classical bit Albert and
Jiang (2014). To conclude our discussion of $\Delta\neq 0$, we use the CQA
method to compute exactly each steady state in the bistable manifold. We begin
by noting that the states $|\Phi_{e/o}\rangle$ in Eq. (90) have definite
photon-number parity, and thus so do the corresponding states of the physical
cavity-$a$ (obtained by tracing over cavity-$b$). Therefore, by uniqueness of
the extremal states, these states must be precisely $\hat{\rho}_{e/o}$:

$\displaystyle\hat{\rho}_{e}$
$\displaystyle=\text{tr}_{b}[|\Phi_{e}\rangle\langle\Phi_{e}|]$
$\displaystyle\hat{\rho}_{o}$
$\displaystyle=\text{tr}_{b}[|\Phi_{o}\rangle\langle\Phi_{o}|].$ (97)

To compute these steady-states, we note that by substituting Eq. (90) into Eq.
(97) we can expand, e.g.

$\displaystyle\hat{\rho}_{e/o}$ $\displaystyle=\frac{N}{N\pm
1}(\hat{\rho}_{a,ss}^{+}\pm\hat{\rho}_{a,ss}^{-})$ (98)

where $N$ is just the normalization constant $N$ for the dark state
$|\psi\rangle_{+}$, which has the exact expression

$\displaystyle N$
$\displaystyle=\,_{1}F_{2}(1/2;1/2-D/2,(1/2-D/2)^{*};|\lambda_{2}/2|^{2}).$
(99)

Inverting the above linear relation, we get

$\displaystyle\hat{\rho}^{+}_{a,ss}$
$\displaystyle=\frac{1}{2}\left[\frac{N+1}{N}\hat{\rho}_{e}+\frac{N-1}{N}\hat{\rho}_{o}\right]$
(100)

This equation immediately leads to exact expressions for $\hat{\rho}_{e/o}$,
which are given in Appendix H. Furthermore, by comparison with Eq. (91),
$\hat{\rho}^{+}_{a,ss}$ also happens to be the steady-state in the presence of
an infinitessimal amount of bistability-breaking single-photon loss.
Therefore, in Eqs. (99-100), CQA is able to smoothly describe the transition
from a Kerr oscillator having two quantum steady states to having only one. In
the weak-driving limit $\lambda_{2}\to 0$, the hypergeometric series defining
$N$ collapses to just the first term (c.f. Eq. (99)), and we get $N\to 1$, so

$\displaystyle\hat{\rho}_{a,ss}^{+}\underset{\lambda_{2}\to
0}{\sim}\hat{\rho}_{e}.$

In contrast, in the strong-driving limit $N$ diverges, and thus

$\displaystyle\hat{\rho}_{a,ss}^{+}\underset{\lambda_{2}\to\infty}{\sim}\frac{\hat{\rho}_{e}+\hat{\rho}_{o}}{2}.$
(101)

A final piece of physical intuition: since $N$ is a function only of the
modulus $|\lambda_{2}|$, the relative bias (towards either $\hat{\rho}_{o/e}$)
is independent of the phase $\phi$ of the drive $\lambda_{2}\equiv
e^{i\phi}|\lambda_{2}|$.

##  VIII Conclusions

In this, work, we have presented a generalization of the coherent quantum
absorber method developed by Stannigel et. al. Stannigel  _et al._ (2012) for
solving the simplest driven Kerr resonator problem. Our generalization
exploited the Segal-Bargmann representation, and allows one to analytically
solve for the steady state of driven-dissipative Kerr cavity models with
nonlinear driving and nonlinear loss. We used these analytic solutions to
describe a host of new physical phenomena, including generalized photon-
blockade phenomena, and new regimes of near quantum bistability. These
phenomena should be experimentally accessible in a number of different
platforms, including superconducting circuit experiments.

Our work naturally suggests many new open questions and directions for future
study. For example, can the new bistable parameter points we have identified
be utilized for quantum-information applications? Are there other forms of
nonlinear dissipation and driving that could also be included in our system
that still leave it amenable to solution via the CQA method? Can this approach
be extended to nonlinear driven-dissipative systems with more than one cavity?

At a fundamental level, there is also the basic question of why the CQA method
is able to yield exact solutions to systems that are on the surface highly
non-trivial (because of strong nonlinearities and driving). Is there some
general physical principle here, or perhaps a dissipative version of
integrability that underlies this method? These are all questions we hope to
explore in future works.

###### Acknowledgements.

This work was supported by the Air Force Office of Scientific Research MURI
program, under grant number FA9550-19-1-0399.

##  Appendix A Circuit QED realization of the model

We now show that it is possible to realize the generalized driven-Kerr
oscillator using cQED devices that already currently exist, in particular, a
superconducting nonlinear asymmetric inductive element (a.k.a. SNAIL device)
Sivak  _et al._ (2019). The Hamiltonian for a SNAIL, as a function of applied
magnetic flux, can be written as (following the notation of Sivak  _et al._
(2019)):

$\displaystyle\hat{H}=\omega\hat{a}^{\dagger}\hat{a}+g_{3}(\Phi)(\hat{a}+\hat{a}^{\dagger})^{3}+g_{4}(\Phi)(\hat{a}+\hat{a}^{\dagger})^{4}.$
(102)

We now introduce time-dependence into the flux parameter $\Phi$ in such a way
that $g_{3}(t)\equiv g_{3}(\Phi(t))$ is oscillating at the cavity frequency
with amplitude $g_{3}^{(0)}$, whereas $g_{4}(t)\equiv g_{4}^{(0)}$ is
essentially constant (c.f. Fig. 1 in Sivak  _et al._ (2019)). In the frame
rotating at the cavity frequency, the time-dependent Hamiltonian then has the
form

$\displaystyle\hat{U}\hat{H}\hat{U}^{\dagger}$
$\displaystyle=3g_{3}^{(0)}\\{(\hat{a}^{\dagger}\hat{a}^{\dagger}\hat{a}+\hat{a}^{\dagger}\hat{a}\hat{a})+(\hat{a}^{\dagger}+\hat{a})\\}$
$\displaystyle+3g_{4}^{(0)}\\{2\hat{a}^{\dagger}\hat{a}^{\dagger}\hat{a}\hat{a}+4\hat{a}^{\dagger}\hat{a}+1\\}$
$\displaystyle\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ +\text{counter-rotating terms.}$ (103)

Under the rotating-wave approximation, if the drives are weak we can neglect
all counter-rotating terms, which yields the effective Kerr Hamiltonian

$\displaystyle\hat{H}_{\rm RWA}$
$\displaystyle=\frac{K}{2}\hat{a}^{\dagger}\hat{a}^{\dagger}\hat{a}\hat{a}-\Delta\hat{a}^{\dagger}\hat{a}$
$\displaystyle+(\Lambda_{1}\hat{a}^{\dagger}+\Lambda_{3}\hat{a}^{\dagger}\hat{a}^{\dagger}\hat{a}+h.c.),$
(104)

where $\Lambda_{1}=\Lambda_{3}=3g_{3}^{(0)}$, and $K=12g_{4}^{(0)}=-\Delta$.
In conclusion, realization of the nonlinear coherent driving effect, for weak
driving strengths, is possible using a superconducting nonlinear asymmetric
inductive element, by modulating its flux parameter at the cavity frequency.
We can also see from this analysis how it would be even harder to realize the
three-photon additional / removal terms $(\hat{a}^{\dagger})^{3},\hat{a}^{3}$
within this scheme, as this would require modulating the external flux $\Phi$
three-times more rapidly (specifically: 18 GHz, for the device considered in
Sivak  _et al._ (2019)).

##  Appendix B Exact realization of new quantum bistable regimes using a two-
cavity non-cascaded setup

In Sec. V.6, we discussed how the generalized driven-dissipative Kerr problem
could be tuned to be arbitrarily close to points in parameter space where we
have true quantum bistability. Exact tuning to a bistable point was not
possible due to the constraint that neither $\kappa_{1}$ nor $\kappa_{2}$
could be made negative.

An exact realization of these quantum bistable points is nonetheless possible
if one works with the two cavity system in Fig. 1. Making one of $\kappa_{1}$
or $\kappa_{2}$ negative now has a simple physical interpretation: we simply
reverse the chirality of one of the waveguides in the absorber setup (see Fig.
9). Reversing the chirality of the (e.g. linearly-coupled) waveguide leads to
the dynamics of the master equation Eq. (4) with the same dissipators but with
the Hamiltonian (c.f. Eq. (5)) changed to

$\displaystyle\hat{H}_{ab}\to\hat{H}_{a}-\hat{H}_{b}+\frac{i\kappa_{1}}{2}(\hat{a}^{\dagger}\hat{b}-h.c.)-\frac{i\kappa_{2}}{2}(\hat{a}^{\dagger}\hat{a}^{\dagger}\hat{b}^{2}-h.c.).$
(105)

Again, using the absorber method, we can solve this master equation in a
manner identical to before, i.e. with
$|\psi\rangle=|\psi_{+}\rangle|0_{-}\rangle$, except now we have
$\kappa_{1}\to-\kappa_{1}$. So the master equation specified by Eq. (105)
constitutes an analytic extension of the steady state to negative values of
$\kappa_{1,2}$ (and thus arbitrary values of $D$), and thus can exhibit
quantum bistability. For a depiction of the setup, see Fig. (9).

In this case Eq. (61) actually becomes a relation for selecting a pure state
in the bistable manifold:

$\displaystyle|\psi_{+}\rangle$ $\displaystyle=\delta
r_{2}|\psi_{+,1}\rangle+\delta r_{1}|\psi_{+,2}\rangle.$ (106)

In this case, $|\psi_{+,j}\rangle$ are the photon-added coherent states of the
symmetric mode defined in Sec. III. The states are perhaps best understood in
the Fock basis. For $\Lambda_{3}=0$,

$\displaystyle|\psi_{+,1}\rangle$ $\displaystyle=\sum_{\leavevmode\nobreak\
m=0\leavevmode\nobreak\ \leavevmode\nobreak\
}^{n}\frac{(2\epsilon\hat{c}_{+}^{\dagger})^{m}}{m!}|-\epsilon\rangle$ (107)
$\displaystyle|\psi_{+,2}\rangle$
$\displaystyle=\sum_{m=n+1}^{\infty}\frac{(2\epsilon\hat{c}_{+}^{\dagger})^{m}}{m!}|-\epsilon\rangle$
(108)

where $|z\rangle$ as usual denotes a coherent state with amplitude $z$. Note
that their sum is Gaussian, i.e. a coherent state, as is expected from
properties of Kummer’s hypergeometric function. In the more general case of
the ”off-diagonal" bistable points (i.e. the $(n,m)$ points with $n\neq m$),
one stabilizes even more exotic states, whose sum may no longer be Gaussian.

Figure 9: Realizing quantum bistability by breaking chirality. The generalized
bistable points in our phase diagram (c.f. Fig. 3) are exactly realizable by
using a two-cavity setup (b) which is not cascaded, i.e. where the chirality
of one of the waveguides is reversed.

##  Appendix C Steady-state Wigner function

We will now rigorously prove the connection between the Wigner function of the
steady-state $\hat{\rho}_{a,ss}$ of the driven Kerr cavity, and the modulus
squared of the SB representation of its purification $|\psi_{+}\rangle$. We
will also show how the calculation generalizes to the case where there are
multiple dark states. The result here relies on a deep fact relating operator
ordering conventions for a quantum-mechanical mode, and the heat semigroup on
the corresponding classical phase space. This was originally pointed out by
Glauber and Cahill in Cahill and Glauber (1969), and we will review the
salient results here. Specifically, given a (possibly non-Hermitian operator)
$\hat{A}$ of a quantum-mechanical mode, define its normally-ordered symbol
$\sigma_{N}$ to be

$\displaystyle\sigma^{N}_{A}(z):=:\hat{A}:\big{|}_{\hat{a}^{\dagger},\hat{a}\mapsto
z,z^{*}},$ (109)

where $:\hat{A}:$ is the operator $\hat{A}$ but re-expressed in normal order,
i.e. with all of the creation operators to the left of the annihilation
operators. Analogously, we can define the symmetrically-ordered symbol:

$\displaystyle\sigma_{A}(z):=:\hat{A}:_{S}\big{|}_{\hat{a}^{\dagger},\hat{a}\mapsto
z,z^{*}},$ (110)

where $:\hat{A}:_{S}$ is the operator $\hat{A}$ but re-expressed according to
the symmetric ordering convention (as defined in Cahill and Glauber (1969)).
The symmetrically-ordered symbol is proportional to the standard Wigner
transform, which can be formally computed via an integral:

$\displaystyle\sigma_{A}(z)$ $\displaystyle\propto\int
d^{2}\xi\,\text{Tr}[e^{\xi^{*}\hat{a}-\xi\hat{a}^{\dagger}}\hat{A}]e^{\xi
z^{*}-\xi^{*}z}$ (111)

For positive semi-definite operators, e.g. a density matrix
$\hat{\rho}\equiv\hat{\rho}^{\dagger}$, the symmetrically-ordered and
normally-ordered symbols coincide with the Wigner- and Q-functions
respectively:

$\displaystyle Q(z)=\frac{1}{\pi}\sigma_{\rho}^{N}(z),\leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
W(z)=\frac{1}{\pi}\sigma_{\rho}(z).$ (112)

What Glauber and Cahill showed in Cahill and Glauber (1969) is that operator
symbols corresponding to different ordering conventions are related by the
heat semigroup. In particular, we have the following theorem:

Theorem (Cahill and Glauber (1969)). Let $\hat{A}$ be a Hilbert-Schmidt
operator (i.e. $\text{Tr}[\hat{A}^{\dagger}\hat{A}]<\infty$). Then its
normally-ordered symbol can be obtained by ”cooling" (i.e. running the heat
equation on) the symmetrically-ordered symbol for a time $t=1/8$, i.e.

$\displaystyle\sigma_{A}^{N}(z)=\int_{\mathbb{C}}d^{2}z^{\prime}\mathcal{K}(t,z,z^{\prime})\sigma_{A}(z^{\prime})\bigg{|}_{t=1/8}$
(113)

where $K(t,z,z^{\prime})$ is the heat kernel on $\mathbb{C}$, which can be
exactly computed and comes out to

$\displaystyle\mathcal{K}(t,z,z^{\prime})=\frac{e^{\frac{-|z-z^{\prime}|^{2}}{4t}}}{4\pi
t}.$ (114)

We can use the theorem above to directly trace-out the ancilla cavity used in
the absorber method in the main text. Following the notation of the main text,
suppose we have an orthonormal basis
$|\psi_{1}\rangle,\cdots,|\psi_{k}\rangle$ of the space of dark states, i.e.

$\displaystyle\hat{c}_{-}|\psi_{j}\rangle=0.$ (115)

Clearly, these states form a $k$-dimensional Bloch sphere, spanned by their
outer-products:

$\displaystyle|\psi_{i}\rangle\langle\psi_{j}|.$ (116)

According to CQA, to obtain the corresponding stationary modes $\hat{M}_{ij}$
of the physical cavity-$a$ Lindbladian, we must trace-out the ancilla mode:

$\displaystyle\hat{M}_{ij}$
$\displaystyle=\text{Tr}_{b}[|\psi_{i}\rangle\langle\psi_{j}|].$ (117)

We now compute the Wigner transform of the above stationary modes. First, we
take advantage of the fact that a dark state factorizes across the two-modes
$\hat{c}_{\pm}$ as $|\psi_{j}\rangle=|\psi_{j,+}\rangle|0_{-}\rangle$:

$\displaystyle\hat{M}_{ij}$
$\displaystyle=\text{Tr}_{b}[(|\psi_{i,+}\rangle\langle\psi_{j,+}|)\,(|0_{-}\rangle\langle
0_{-}|)].$ (118)

Letting $\sigma_{+,ij}(z)$ denote the symmetrically-ordered symbol (i.e.
Wigner transform) of the outer-product $|\psi_{i,+}\rangle\langle\psi_{j,+}|$,
and $\sigma_{-}(z)$ denote the symmetrically-ordered symbol of the vacuum
state $|0_{-}\rangle\langle 0_{-}|$, let

$\displaystyle\sigma_{ij}(z)$ (119)

denote the symmetrically-ordered symbol (i.e. the Wigner transform) of the
stationary mode $\hat{M}_{ij}$. We can then rewrite the expression for the
partial trace completely in terms of symmetrically-ordered symbols: in this
case, the partial trace becomes an integral, and the symmetrized-
antisymmetrized nature of the input states means that the integral convolves
the operator symbols. The symbol $\sigma_{-}(z)$ then acts as a Gaussian
filter for the symbol $\sigma_{+,ij}(z)$:

$\displaystyle\sigma_{ij}(z)$
$\displaystyle=\int_{\mathbb{C}}d^{2}z^{\prime}\leavevmode\nobreak\
\sigma_{-}\bigg{(}\frac{z-z^{\prime}}{\sqrt{2}}\bigg{)}\sigma_{+,ij}\bigg{(}\frac{z+z^{\prime}}{\sqrt{2}}\bigg{)}$
(120) $\displaystyle=\int_{\mathbb{C}}d^{2}z_{+}\leavevmode\nobreak\
\frac{2e^{-2|z_{+}-\sqrt{2}z|^{2}}}{\pi}\sigma_{+,ij}(z_{+}),$ (121)

where we have defined symmetrized and anti-symmetrized phase-space variables
$z_{\pm}\equiv(z\pm z^{\prime})/\sqrt{2}$. The above filtering operation is
the same exact operation which ”reorders" a normally-ordered symbol into a
symmetrically-ordered symbol, up to a rescaling of the phase space
$z\mapsto\sqrt{2}z$. Indeed, we can rewrite it in terms of the heat kernel:

$\displaystyle\sigma_{ij}(z)$
$\displaystyle=2\int_{\mathbb{C}}d^{2}u\,\mathcal{K}(t,\sqrt{2}z,u)\sigma_{+,ij}(u)\bigg{|}_{t=1/8}$
(122) $\displaystyle=2\sigma^{N}_{+,ij}(\sqrt{2}z),$ (123)

where $\sigma^{N}_{+,ij}$ is the normally-ordered symbol of the mode
$|\psi_{i,+}\rangle\langle\psi_{j,+}|$. This form highly constrains the Wigner
function of the steady state of any single-mode system with single-photon loss
that is solvable via CQA.

In particular, now we can compute the symbol exactly in terms of the Segal-
Bargmann representation. It is easy to show that the normally-ordered symbol
of an operator has the simple form:

$\displaystyle\sigma_{A}^{N}(z)=\langle z|\hat{A}|z\rangle$ (124)

where $|z\rangle$ denotes a coherent state with amplitude $z$. By expanding
$\hat{A}$ in terms of outer-products as

$\displaystyle\hat{A}=\sum_{ij}\alpha_{ij}|\psi_{i}\rangle\langle\psi_{j}|,$
(125)

and utilizing the property $\psi_{\rm SB}(z)=\langle z^{*}|\psi\rangle
e^{-|z|^{2}/2}$, Bargmann in Bargmann (1967) was able to show that this
implies

$\displaystyle\sigma_{A}^{N}(z)=\sum_{ij}\alpha_{ij}\psi_{i,\rm
SB}(z^{*})\psi_{j,\rm SB}(z^{*})e^{-|z|^{2}}.$ (126)

By substituting the exact expression for the normally-ordered symbol into Eq.
(123), we can finally state the main result utilized in the main text:

$\displaystyle\sigma_{ij}(z)=2\psi_{i,\rm SB}(\sqrt{2}z^{*})\psi_{j,\rm
SB}(\sqrt{2}z^{*})e^{-2|z|^{2}}.$ (127)

By taking linear combinations of the above stationary modes, the Wigner
function of any stationary density matrix of the physical cavity $a$ has the
closed-form

$\displaystyle W_{a,ss}(z)=\frac{2}{\pi}\sum_{ij}\alpha_{ij}\psi_{i,\rm
SB}(\sqrt{2}z^{*})\psi_{j,\rm SB}^{*}(\sqrt{2}z^{*})e^{-2|z|^{2}},$ (128)

where $\alpha_{ij}=\alpha_{ji}^{*}$ is a positive semi-definite matrix with
unit trace.

We can also write the CQA ansatz in a manifestly positive-form. Letting
$\\{p_{j}\\}$ denote the eigenvalues of the positive semi-definite matrix
$\\{\alpha_{ij}\\}$, and letting $\\{\phi_{j}(z)\\}$ denote the Segal-Bargmann
representations of the corresponding eigenvectors, the Wigner function can be
equivalently written as

$\displaystyle
W_{a,ss}(z)=\frac{2}{\pi}\sum_{j}p_{j}|\phi_{j}(\sqrt{2}z^{*})|^{2}e^{-2|z|^{2}}\geq
0,$ (129)

where normalization forces $\sum_{j}p_{j}=1$.

As a simple example, when the dark-state subspace is one-dimensional, the
Wigner function is the squared-modulus of the SB representation of the unique,
normalized dark state in that subspace:

$\displaystyle
W_{a,ss}(z)=\frac{2}{\pi}|\phi(\sqrt{2}z^{*})|^{2}e^{-2|z|^{2}}.$ (130)

In summary, we have derived an exact, closed-form expression for the steady-
state Wigner function of a cavity that is solvable via CQA. What is most
striking from this analysis is the absence of any consideration of the
Hamiltonian of the cavity: the CQA method, if it works, will predict that the
Wigner function will be positive-definite, simply as a consequence of the
presence of single-photon loss in the system.

####  C.0.1 Other phase-space representations of the steady state

One can also obtain the steady-state $P$ function using the results of the
above analysis. Accordingly, we define the generalized Weierstrass transform
$\mathcal{W}_{t}\equiv e^{-t\Delta}$ as the following integral transform:

$\displaystyle(\mathcal{W}_{t}f)(z)\equiv\int_{\mathbb{C}}d^{2}u\mathcal{K}(t,z,u)f(u).$
(131)

Let $(S_{\alpha}f)(z)\equiv f(\alpha z)$ also denote the linear operator which
rescales the function argument. We then have the following commutation
relation:

$\displaystyle(S_{\sqrt{s}}\mathcal{W}_{t}f)(z)$
$\displaystyle=\int_{\mathbb{C}}d^{2}u\frac{e^{-\frac{|u-\sqrt{s}z|}{4t}}}{4\pi
t}f(u)$
$\displaystyle=\int_{\mathbb{C}}d^{2}v\frac{e^{-\frac{|v-z|^{2}}{4t/s}}}{4\pi
t/s}f(\sqrt{s}z)=(\mathcal{W}_{t/s}S_{\sqrt{s}}f)(z)$ (132)

From Eq. (132), it follows immediately that
$\mathcal{W}_{t}S_{\sqrt{s}}=S_{\sqrt{s}}\mathcal{W}_{s\cdot t}$. With this
identity in mind, we can solve for the steady-state $P$-function. Letting
$Q_{+},P_{+}$ denote the $P$\- and $Q$-functions of the single-mode pure state
$|\psi_{+}\rangle$, we have

$\displaystyle P_{a,ss}(z)$
$\displaystyle=(\mathcal{W}_{-1/8}S_{\sqrt{2}}\,Q_{+})(z)=(S_{\sqrt{2}}\mathcal{W}_{-1/4}Q_{+})(z)$
$\displaystyle=2P_{+}(\sqrt{2}z).$ (133)

Therefore, the $P$-function of the Kerr resonator is simply the $P$-function
of the single-mode pure state $|\psi_{+}\rangle$, up to a rescaling of phase
space $z\mapsto\sqrt{2}z$.

##  Appendix D “Gauge-invariance” of the dark state

Eq. (34) in the main text makes it manifest that there exist two distinct
(non-unitary) gauge choices in which the troublesome two-photon term vanishes
in the dark state equation Eq. (18). In the main text, we solved the dark
state equations in the plus-gauge. This leads to the question of what would
happen if we solved the dark state conditions in the minus-gauge. If we had
solved for the dark state in the minus gauge, we would have ended up with the
solution

$\displaystyle\widetilde{\psi}_{+,\text{SB}}(z)$ $\displaystyle\equiv
e^{-\epsilon_{-}z}\,_{1}F_{1}\bigg{(}-\frac{\lambda_{1}+\epsilon_{-}D}{\epsilon_{-}-\epsilon_{+}};-D;(\epsilon_{-}-\epsilon_{+})z\bigg{)}.$
(134)

However, the results are gauge-invariant, as we can write

$\displaystyle\widetilde{\psi}_{+,\text{SB}}(z)=e^{-\epsilon_{-}z+(\epsilon_{-}-\epsilon_{+})z}$
$\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\
\,{}_{1}F_{1}\bigg{(}-D-\frac{\lambda_{1}+\epsilon_{-}D}{\epsilon_{+}-\epsilon_{-}};-D;(\epsilon_{+}-\epsilon_{-})z\bigg{)}$
(135) $\displaystyle=e^{-\epsilon_{+}z}$ $\leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\,{}_{1}F_{1}\bigg{(}-\frac{(\epsilon_{+}-\epsilon_{-})D}{\epsilon_{+}-\epsilon_{-}}-\frac{\lambda_{1}+\epsilon_{-}D}{\epsilon_{+}-\epsilon_{-}};-D;(\epsilon_{+}-\epsilon_{-})z\bigg{)}$
$\displaystyle=e^{-\epsilon_{+}z}\,_{1}F_{1}\bigg{(}-\frac{\lambda_{1}+\epsilon_{+}D}{\epsilon_{+}-\epsilon_{-}};-D;(\epsilon_{+}-\epsilon_{-})z\bigg{)}=\psi_{+,\rm
SB},$ (136)

where $\psi_{+,\rm SB}$ is the dark state Eq. (33) in the main text, and in
the first line (Eq. (135)) we utilized Kummer’s transformation (see e.g. Ref.
Brychkov (2008)), which is a fundamental symmetry of the confluent
hypergeometric differential equation:

$\,{}_{1}F_{1}(r_{1};r_{2};z)$
$\displaystyle=e^{z}\,_{1}F_{1}(r_{2}-r_{1};r_{2};-z).$ (137)

##  Appendix E Beyond nonlinear single-photon driving: breakdown of CQA

The most general Kerr Hamiltonian, that is, containing all possible terms of
lower order than the Kerr nonlinearity, is

$\displaystyle\hat{H}_{a}$
$\displaystyle=\frac{K}{2}\hat{a}^{\dagger}\hat{a}^{\dagger}\hat{a}\hat{a}-\Delta\hat{a}^{\dagger}\hat{a}$
$\displaystyle+\bigg{[}\bigg{(}\Lambda_{1}\hat{a}^{\dagger}+\frac{\Lambda_{2}}{2}\hat{a}^{\dagger}\hat{a}^{\dagger}+\Lambda_{3}\hat{a}^{\dagger}\hat{a}^{\dagger}\hat{a}+\Lambda_{4}\,\hat{a}^{\dagger}\hat{a}^{\dagger}\hat{a}^{\dagger}\bigg{)}+h.c.\bigg{]}.$
(138)

This begs the question of why everything of degree three or lower is exactly
solvable by CQA, except the $(\hat{a}^{\dagger})^{3}$ term. The explanation
for this is rather simple: when acting the cascaded Hamiltonian on the dark
state ansatz, one gets the equation

$\displaystyle(\hat{c}_{-}^{\dagger}\mathcal{\hat{H}}_{+}+\frac{\Lambda_{4}}{\sqrt{8}}(\hat{c}_{-}^{\dagger})^{3})|\psi_{+}\rangle=0$
(139)

where

$\displaystyle\mathcal{\hat{H}}_{+}$
$\displaystyle\to\mathcal{\hat{H}}_{+}+\frac{3\Lambda_{4}}{\sqrt{8}}(\hat{c}_{+}^{\dagger})^{2}.$
(140)

The shift in $\hat{\mathcal{H}}_{+}$ is innocuous. However, the term cubic in
$\hat{c}_{-}^{\dagger}$ is lethal: both terms in Eq. (139) have respectively
one and three photons in the minus mode, and are thus generically orthogonal.
Thus the only solution to the dark state condition Eq. (139) is
$|\psi_{+}\rangle=0$.

##  Appendix F Stationary density matrix and moments of a driven Kerr cavity

We show here how to compute exact analytic expressions from our steady-state
solution for density matrix elements in the Fock basis, as well as normal-
ordered cavity moments. Although the expressions obtained here are
considerably more complex/physically opaque, this will allow us to make
contact with older results obtained via $P$-function methods Bartolo  _et al._
(2016); Elliott and Ginossar (2016). Expanding the purification of the density
matrix $|\psi\rangle=|\psi_{+}\rangle\otimes|0_{-}\rangle$, and writing the
symmetric component in the SB representation yields

$\displaystyle\psi_{+,\rm
SB}(z)\equiv\sum_{l=0}^{\infty}\psi_{l}\frac{z^{l}}{l!},$ (141)

which implicitly defines coefficients $\psi_{l}\equiv\psi_{+,\rm SB}^{(l)}(0)$
which are the derivatives of the Bargmann state evaluated at the origin $z=0$
in phase space. In the special cases $\lambda_{1}\equiv\lambda_{3}\equiv 0$,
reproducing results in Bartolo  _et al._ (2016)), or the more generic regime
$\lambda_{2}\equiv 0$, which represents new results, we can actually evaluate
the sums, resulting in compact, closed-form expressions.

####  F.0.1 Steady-state density matrix

In terms of these Taylor coefficients, the steady state density matrix can be
computed in the Fock basis:

$\displaystyle\langle m|\hat{\rho}_{a,ss}|n\rangle$ (142)
$\displaystyle=\sum_{l=0}^{\infty}\langle
m,l|\bigg{(}\sum_{j,k=0}^{\infty}\frac{\psi_{j}\psi_{k}^{*}}{j!k!}(\hat{c}_{+}^{\dagger})^{j}|0\rangle\langle
0|\hat{c}_{+}^{k}\bigg{)}|n,l\rangle$
$\displaystyle=\frac{1}{(2^{m+n}n!m!)^{1/2}}$
$\displaystyle\sum_{l,j,k=0}^{\infty}\frac{\psi_{j}\psi_{k}^{*}}{j!k!}\langle
0|\frac{(\hat{c}_{+}-\hat{c}_{-})^{l}}{\sqrt{2^{l}l!}}(\hat{c}_{+}+\hat{c}_{-})^{m}(\hat{c}_{+}^{\dagger})^{j}|0\rangle$
$\displaystyle\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \cdot\langle
0|\hat{c}_{+}^{k}(\hat{c}_{+}^{\dagger}+\hat{c}_{-}^{\dagger})^{n}\frac{(\hat{c}_{+}^{\dagger}-\hat{c}_{-}^{\dagger})^{l}}{\sqrt{2^{l}l!}}|0\rangle$
$\displaystyle=\frac{1}{(2^{m+n}n!m!)^{1/2}}$
$\displaystyle\sum_{j,k,l=0}^{\infty}\frac{\psi_{j}\psi_{k}^{*}}{j!k!}\frac{1}{2^{l}l!}\langle
0|\hat{c}_{+}^{m+l}(\hat{c}_{+}^{\dagger})^{j}|0\rangle\langle
0|\hat{c}_{+}^{k}(\hat{c}_{+}^{\dagger})^{n+l}|0\rangle.$ (143)

Using identities of the form $\langle
0|\hat{c}_{+}^{m+l}(\hat{c}_{+}^{\dagger})^{j}|0\rangle=\delta_{m+l,j}j!$,
etc., we get the remarkably simple result:

$\displaystyle\langle
m|\hat{\rho}_{a,ss}|n\rangle=\frac{1}{\sqrt{2^{m+n}n!m!}}\sum_{l=0}^{\infty}\frac{\psi_{m+l}\psi_{n+l}^{*}}{2^{l}l!}.$
(144)

This expression matches similar expressions obtained using complex-P
solutions, as we will see later in this section.

####  F.0.2 Cavity moments

We can also express the normally-ordered moments of a driven Kerr cavity
exactly in terms of the scaled Fock-state amplitudes $\psi_{l}$. The
calculation is slightly more straightforward:

$\displaystyle\text{Tr}[\hat{\rho}_{a,ss}(a^{\dagger})^{n}a^{m}]$
$\displaystyle=\langle\psi|(a^{\dagger})^{n}a^{m}|\psi\rangle$
$\displaystyle=\frac{1}{\sqrt{2^{m+n}}}\langle\psi|(\hat{c}_{+}^{\dagger}+\hat{c}_{-}^{\dagger})^{n}(\hat{c}_{+}+\hat{c}_{-})^{m}|\psi\rangle,$
(145)

where $|\psi\rangle$, as before, is the purification of the density matrix
obtained from the absorber method. Expanding the dark state yields

$\displaystyle\text{Tr}[\hat{\rho}_{a,ss}(a^{\dagger})^{n}a^{m}]$
$\displaystyle=\frac{1}{\sqrt{2^{m+n}}}\sum_{j,k=0}^{\infty}\frac{\psi_{j}^{*}\psi_{k}}{j!k!}\langle
0|\hat{c}_{+}^{j}(\hat{c}_{+}^{\dagger})^{n}\hat{c}_{+}^{m}(\hat{c}_{+}^{\dagger})^{k}|0\rangle$
$\displaystyle=\frac{1}{\sqrt{2^{m+n}}}\sum_{j,k=0}^{\infty}\frac{\psi_{j}^{*}\psi_{k}}{\sqrt{j!k!}}\langle
j_{+}|(\hat{c}_{+}^{\dagger})^{n}\hat{c}_{+}^{m}|k_{+}\rangle.$ (146)

Defining a new variable $l$ such that $j\equiv n+l$, we find that $k=m+l$, and
that furthermore $l\geq 0$. So our sum simplifies to

$\displaystyle\text{Tr}[\hat{\rho}_{a,ss}(a^{\dagger})^{n}a^{m}]$
$\displaystyle\frac{1}{\sqrt{2^{m+n}}}\sum_{l=0}^{\infty}\frac{\psi_{n+l}^{*}\psi_{m+l}}{\sqrt{(m+l)!(n+l)!}}\frac{\sqrt{(m+l)!}}{\sqrt{l!}}\frac{\sqrt{(n+l)!}}{\sqrt{l!}}.$
(147)

We thus obtain the simple result

$\displaystyle\text{Tr}[\hat{\rho}_{a,ss}(a^{\dagger})^{n}a^{m}]$
$\displaystyle=\frac{1}{\sqrt{2^{m+n}}}\sum_{l=0}^{\infty}\frac{\psi_{m+l}\psi_{n+l}^{*}}{l!}.$
(148)

This is the formula used to produce exact-solution plots of average photon
number in Fig. 4; a similar-looking expression was derived independently in
Bartolo  _et al._ (2016), using $P$-function methods.

####  F.0.3 Normalization

Throughout this section, we have assumed that the normalization of
$|\psi_{+}\rangle$ is known. Supposing that this is not the case, and
$|\psi_{+}\rangle$ is written instead in the form

$\displaystyle\psi_{+,\rm
SB}(z)=\frac{1}{\sqrt{N}}\sum_{l=0}^{\infty}\widetilde{\psi_{l}}\frac{z^{l}}{l!},$
(149)

we can write an exact expression for $N$:

$\displaystyle
N=\sum_{l=0}^{\infty}\frac{\widetilde{\psi_{l}}}{\sqrt{l!}}\frac{\widetilde{\psi_{l}}^{*}}{\sqrt{l!}}=\sum_{l=0}^{\infty}\frac{|\widetilde{\psi_{l}}|^{2}}{l!}.$
(150)

####  F.0.4 Expression for $\psi_{l}$ in general regime

The scaled Fock-state amplitudes $\psi_{l}$ can be computed in closed-form in
terms of the Gauss hypergeometric function. We can then utilize this closed
form to show that our exact expressions derived here agree with earlier
solutions Bartolo  _et al._ (2016); Elliott and Ginossar (2016) in the limit
of $\Lambda_{3}\to 0$:

$\displaystyle\psi_{l}$ $\displaystyle\equiv\partial^{l}\psi_{+,\rm SB}(0).$
(151)

The above quantity is particularly difficult to evaluate in the general case,
so we evaluate instead

$\displaystyle\xi_{l}$ $\displaystyle\equiv\partial^{l}\xi_{+,\rm SB}(0).$
(152)

where $|\xi\rangle$ is the displaced dark state in the main text. Eq. (152)
then represents the Fock-state amplitudes of the purification of the displaced
steady-state
$\hat{\rho}^{\prime}\equiv\hat{D}_{\alpha}\hat{\rho}_{a,ss}\hat{D}_{\alpha}^{\dagger}$,
where $\alpha$ is defined in the main text and vanishes when $\Lambda_{3}\to
0$.

Expanding $\xi_{+,\rm SB}(z)\equiv\Theta(z)\phi(z)$, where
$\Theta(z)\equiv\exp(-\theta(z))$ is the non-unitary gauge transformation in
the main text. Expanding via the Leibniz rule, we get

$\displaystyle\xi_{l}$
$\displaystyle=\sum_{n=0}^{l}\binom{l}{n}\partial^{l-n}\Theta(0)\partial^{n}\phi(0).$
(153)

Plugging in $\Theta(z)\equiv e^{-\epsilon_{+}z}$ and
$\phi(z)=\,_{1}F_{1}(-r_{1};-r_{2};(\epsilon_{+}-\epsilon_{-})z)$, we get

$\displaystyle\partial^{k}\Theta(0)$ $\displaystyle=(-\epsilon_{+})^{k}$ (154)
$\displaystyle\partial^{k}\phi(0)$
$\displaystyle=\frac{(-r_{1})_{k}}{(-r_{2})_{k}}(\epsilon_{+}-\epsilon_{-})^{k}$
(155)

So, in total, we get

$\displaystyle\xi_{l}$
$\displaystyle=\sum_{n=0}^{l}\binom{l}{n}\frac{(-r_{1})_{n}}{(-r_{2})_{n}}(-\epsilon_{+})^{l-n}(\epsilon_{+}-\epsilon_{-})^{n}$
$\displaystyle=(-\epsilon_{+})^{l}\sum_{n=0}^{l}(-1)^{n}\binom{l}{n}\frac{(-r_{1})_{n}}{(-r_{2})_{n}}\bigg{(}1-\frac{\epsilon_{-}}{\epsilon_{+}}\bigg{)}^{n}$
$\displaystyle=(-\epsilon_{+})^{l}\sum_{n=0}^{l}\frac{(-l)_{n}(-r_{1})_{n}}{(-r_{2})_{n}}\frac{(1-\frac{\epsilon_{-}}{\epsilon_{+}})^{n}}{n!}.$
$\displaystyle=(-\epsilon_{+})^{l}\sum_{n=0}^{\infty}\frac{(-l)_{n}(-r_{1})_{n}}{(-r_{2})_{n}}\frac{(1-\frac{\epsilon_{-}}{\epsilon_{+}})^{n}}{n!}.$
(156)

Therefore, we have a closed-form expression for the scaled Fock-state
amplitudes of the displaced dark state:

$\displaystyle\xi_{l}$
$\displaystyle=(-\epsilon_{+})^{l}\,_{2}F_{1}(-l,-r_{1};-r_{2};1-\tfrac{\epsilon_{-}}{\epsilon_{+}}).$
(157)

In the limit $\Lambda_{3}\to 0$, $\epsilon_{+}\to-\epsilon_{-}$, and so, as in
the main text, defining $\epsilon\equiv\epsilon_{+}$, we get

$\displaystyle\psi_{l}$ $\displaystyle\underset{\Lambda_{3}\to
0}{\sim}(-\epsilon)^{l}\,_{2}F_{1}(-l,-r_{1};-r_{2};2).$ (158)

where we are implicitly utilizing the fact that $\xi_{l}\to\psi_{l}$ in this
limit, as the displacement parameter $\alpha$ vanishes in the limit
$\Lambda_{3}\to 0$. From Eq. (158), it is straightforward to recover the
previous solutions Bartolo  _et al._ (2016); Elliott and Ginossar (2016) of
the Kerr resonator in the limit $\Lambda_{3}\equiv 0$.

##  Appendix G Exact results when the non-unitary gauge transformation is
trivial

The series expressions derived in Appendix F have simple closed forms when we
have $\lambda_{2}\equiv 0$, which, for $\Lambda_{2}\neq 0$, represents
previously unexplored physics. We emphasize the generic nature of this regime,
in that there are 8 real parameters to play with: $\Lambda_{1},\Lambda_{3}$,
and $K,\Delta,\kappa_{1},\kappa_{2}$. In this limit, the displaced SB
wavefunction Eq. (32) is purely hypergeometric:

$\displaystyle\xi_{+,\rm SB}(z)$
$\displaystyle=\frac{1}{N^{1/2}}\,_{1}F_{1}(-r_{1};-r_{2};-\lambda_{3}z),$
(159)

where $r_{1}\equiv-\lambda_{1}/\lambda_{3}$, and $r_{2}\equiv D$. In this
case, the coefficients $\xi_{l}$ of the displaced steady state simplify to
ratios of Pochhammer symbols:

$\displaystyle\xi_{l}=\frac{1}{N^{1/2}}\frac{(-r_{1})_{l}}{(-r_{2})_{l}}(-\lambda_{3})^{l},$
(160)

where the Pochhammer symbol is defined as
$(z)_{l}\equiv\Gamma(z+l)/\Gamma(z)$. Therefore, the normalization is
computable in closed form:

$\displaystyle N$
$\displaystyle=\sum_{l=0}^{\infty}\frac{(-r_{1})_{l}(-r_{1}^{*})_{l}}{(-r_{2})_{l}(-r_{2}^{*})_{l}}\frac{|\lambda_{3}|^{2l}}{l!}$
$\displaystyle\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\
=\,_{2}F_{2}(-r_{1},-r_{1}^{*};-r_{2},-r_{2}^{*};|\lambda_{3}|^{2}).$ (161)

Here, $\,{}_{p}F_{q}(a_{1}\cdots a_{p};b_{1},\cdots b_{q};z)$ denotes the
generalized hypergeometric function (see, e.g. Brychkov (2008)). The
normalization of the steady-state Wigner function is thus exactly computable:

$\displaystyle W_{a,ss}(z-\alpha)$
$\displaystyle=\frac{2|_{1}F_{1}(-r_{1};-r_{2};-\sqrt{2}\lambda_{3}z^{*})|^{2}e^{-2|z|^{2}}}{\pi\,_{2}F_{2}(-r_{1},-r_{1}^{*};-r_{2},-r_{2}^{*};|\lambda_{3}|^{2})^{1/2}},$
(162)

where here $\alpha\equiv\alpha_{+}/\sqrt{2}$ is the appropriately normalized
displacement factor given in Eq. (32) in the main text. We now move on to
compute the matrix elements of the density matrix in the displaced frame
(here, $\hat{D}_{\alpha}\equiv e^{-\alpha\hat{a}^{\dagger}-h.c.}$ is the
standard displacement operator of the physical cavity).

$\displaystyle\langle
m|\hat{D}_{\alpha}\hat{\rho}_{a,ss}\hat{D}_{\alpha}^{\dagger}|n\rangle=\frac{(-\lambda_{3})^{m}(-\lambda_{3}^{*})^{n}}{N\sqrt{2^{m+n}n!m!}}$
$\displaystyle\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\
\cdot\sum_{l=0}^{\infty}\frac{(-r_{1})_{m+l}(-r_{1}^{*})_{n+l}}{(-r_{2})_{m+l}(-r_{2}^{*})_{n+l}}\frac{(|\lambda_{3}|^{2}/2)^{l}}{l!}.$

Utilizing the identity $(z)_{m+l}=(z)_{m}(z+m)_{l}$, the sum closes, and we
get

$\displaystyle\langle
m|\hat{D}_{\alpha}\hat{\rho}_{a,ss}\hat{D}_{\alpha}^{\dagger}|n\rangle=\frac{\xi_{m}\xi_{n}^{*}}{\sqrt{2^{m+n}n!m!}}$
$\displaystyle\cdot\,_{2}F_{2}(n-r_{1},m-r_{1}^{*};m-r_{2};n-r_{2}^{*};|\lambda_{3}|^{2}/2)$
(163)

As for the normally-ordered cavity moments in the displaced frame, in a
similar fashion, we get an analogous closed-form in terms of a generalized
hypergeometric function:

$\displaystyle\text{Tr}[\hat{D}_{\alpha}\hat{\rho}_{a,ss}\hat{D}_{\alpha}^{\dagger}(\hat{a}^{\dagger})^{n}\hat{a}^{m}]=\frac{\xi_{m}\xi_{n}^{*}}{\sqrt{2^{m+n}}}$
$\displaystyle\cdot\,_{2}F_{2}(m-r_{1},n-r_{1}^{*};m-r_{2};n-r_{2}^{*};|\lambda_{3}|^{2})$
(164)

##  Appendix H Exact results in the parity-conserving regime

We now will complete the process started in Section VII, namely that of
tracing-out the ancilla resonator for each of the dark steady states obtained
by the CQA method. We begin with the formula in Appendix F on unique steady
states:

$\displaystyle\langle m|\hat{\rho}_{a,ss}|n\rangle$
$\displaystyle=\frac{1}{\sqrt{2^{m+n}m!n!}}\sum_{l=0}^{\infty}\frac{\psi_{m+l}\psi_{n+l}^{*}}{2^{l}l!}$
(165)

note that, as a direct consequence of $\psi_{2l-1}\equiv 0$, we have

$\displaystyle\langle 2j+1|\hat{\rho}_{a,ss}|2k\rangle=\langle
2j|\hat{\rho}_{a,ss}|2k+1\rangle=0,$ (166)

as each term in the sum over $l$ would identically vanish in these cases. In
summary,

$\displaystyle\hat{\Pi}_{e}\hat{\rho}_{a,ss}\hat{\Pi}_{o}=\hat{\Pi}_{o}\hat{\rho}_{a,ss}\hat{\Pi}_{e}=0,$
(167)

where $\hat{\Pi}_{e/o}$ are the projections onto the subspaces of the
resonator Hilbert space spanned by even/odd photon number states.

Therefore, by taking matrix elements on both sides of Eq. (100) in the main
text, one obtains

$\displaystyle\langle m|\hat{\rho}_{e}|n\rangle$
$\displaystyle=\frac{2N}{N+1}\langle
m|\hat{\rho}_{a,ss}|n\rangle,\leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ m,n\leavevmode\nobreak\ \leavevmode\nobreak\ \text{even}$
(168) $\displaystyle\langle m|\hat{\rho}_{o}|n\rangle$
$\displaystyle=\frac{2N}{N-1}\langle
m|\hat{\rho}_{a,ss}|n\rangle,\leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ m,n\leavevmode\nobreak\ \leavevmode\nobreak\ \text{odd}.$
(169)

Therefore, to compute the steady states $\hat{\rho}_{e/o}$, it suffices to
compute matrix elements of $\hat{\rho}_{a,ss}$. We note that this was done in
Bartolo  _et al._ (2016) (as this represents the unique steady-state regime
$\kappa_{1}\neq 0$), and so we’re technically done, as we could simply cite
the result here.

Figure 10: Limiting behavior of bistable states. (a). We plot the fidelity of
$\hat{\rho}_{e}$ (c.f. Eq. (168)) with an even cat state with amplitude
$\alpha=i\sqrt{\lambda_{2}/2}$ (solid line) and the vacuum state (dashed
line). Corresponding results from exact diagonalization are also given (black
dots) (b) We plot the fidelity of $\hat{\rho}_{o}$ (c.f. Eq. (169)) with an
amplitude-$\alpha$ odd cat state (solid line) and a 1-photon Fock state
(dashed line). Corresponding results from exact diagonalization are also given
(black dots). Parameter choices: In both plots, $\Lambda_{2}=5K$,
$\kappa_{2}=K$, and $\Lambda_{1},\kappa_{1}\equiv 0$.

For completeness, however, we show that the calculation of the expressions on
the RHS’s of Eqs. (168-169) can be reproduced in a straightforward manner
within the quantum absorber formalism. Assuming $m\equiv 2j$, $n\equiv 2k$ are
both even, we have

$\displaystyle\sum_{l=0}^{\infty}\frac{\psi_{m+l}\psi_{n+l}^{*}}{2^{l}l!}=\sum_{l=0}^{\infty}\frac{\psi_{2(j+l)}\psi_{2(k+l)}^{*}}{2^{2l}(2l)!}$
$\displaystyle=\psi_{2j}\psi_{2k}^{*}\sum_{l=0}^{\infty}\frac{(j+\frac{1}{2})_{l}(k+\frac{1}{2})_{l}}{(j-r_{2})_{l}(k-r_{2}^{*})_{l}}\frac{|\lambda_{2}|^{2l}}{2^{2l}(2l)!}$
$\displaystyle=\psi_{2j}\psi_{2k}^{*}\sum_{l=0}^{\infty}\frac{(j+\frac{1}{2})_{l}(k+\frac{1}{2})_{l}}{(j-r_{2})_{l}(k-r_{2}^{*})_{l}(\frac{1}{2})_{l}}\frac{|\lambda_{2}/4|^{2l}}{l!}$
(170)

Therefore, in total we have

$\displaystyle\sum_{l=0}^{\infty}\frac{\psi_{m+l}\psi_{n+l}^{*}}{2^{l}l!}$
(171) $\displaystyle\leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\
=\psi_{2j}\psi_{2k}^{*}\cdot{}_{2}F_{3}{\left[\genfrac{.}{.}{0.0pt}{}{j+\frac{1}{2}\mskip
20.0mu\leavevmode\nobreak\ k+\frac{1}{2}\leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\ }{j-r_{2}\mskip
20.0muk-r_{2}^{*}\mskip 20.0mu\frac{1}{2}};|\lambda_{2}/4|^{2}\right]}.$ (172)

Assuming $m\equiv 2j+1$, $n\equiv 2k+1$ are both odd, we have

$\displaystyle\sum_{l=0}^{\infty}\frac{\psi_{m+l}\psi_{n+l}^{*}}{2^{l}l!}=\sum_{l=0}^{\infty}\frac{\psi_{2(j+l+1)}\psi_{2(k+l+1)}^{*}}{2^{2l+1}(2l+1)!}$
$\displaystyle=4\sum_{l=1}^{\infty}\frac{l\psi_{2(j+l)}\psi_{2(k+l)}^{*}}{2^{2l}(2l)!}$
$\displaystyle=4\psi_{2j}\psi_{2k}^{*}\sum_{l=1}^{\infty}\frac{l(\frac{1}{2}+j)_{l}(\frac{1}{2}+k)_{l}}{(j-r_{2})_{l}(k-r_{2}^{*})_{l}(\frac{1}{2})_{l}}\frac{|\lambda_{2}/4|^{2l}}{l!}$