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

We present a new approach for deriving exact, closed-form solutions for the steady state of a wide class of driven-dissipative nonlinear resonators that is distinct from more common positive-$P$ function methods. Our method generalizes the coherent quantum absorber approach of Stannigel et al. 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.

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 positive-P function methods. Our method generalizes the coherent quantum absorber approach of Stannigel et al. [1] 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 [2], one can exactly solve for the steady state of this system using a positive-P phase space representation. Later work showed that models including two-photon driving and loss are also solvable using this technique [2][3][4]. These driven nonlinear cavity systems have renewed relevance, as they can be directly implemented in superconducting circuit QED setups (see, e.g., [5][6][7][8][9]). Their ability to exhibit multiple steady states has utility in quantum information processing [10][11][12].
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. [1] were able to reproduce the exact solution of Ref. [2] using a simple, purely algebraic approach. While extremely elegant, it was unclear whether this approach could be extended to more complex problems.
In this paper, we show that such an extension is indeed possible: the "coherent quantum absorber" (CQA) method of Ref. [1] 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 [13][14][15][16]. This enables non-trivial transformations that are crucial for finding exact solutions. Our approach yields several new Nonlinear mode (Generalized) Kerr resonator Equivalent Cascaded network (a) (b) Figure 1. (a) Generalized driven Kerr cavity problem, where a single interacting bosonic mode is subject to linear and nonlinear coherent drives Λj, as well as independent one and two photon loss (rates κ1, κ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.
insights. We find and describe new parameter regimes where the steady-state exhibits a surprising generalized photon-blockade phenomenon. We also find an infinite number of points in parameter space where our generalized driven Kerr system exhibits quantum bistability, despite the lack of photon-number parity conservation; these 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 closedform 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 We work in a rotating frame in which the drives are stationary. Here ∆ is the drive detuning, K is the Kerr nonlinearity, and Λ 1 , Λ 2 , Λ 3 are the complex amplitudes of the coherent one-, two-, and three-photon drives. The full dissipative dynamics includes one and two photon loss processes, and is described by the Lindblad master equation where D[X]ρ ≡XρX † − (1/2) X †X ,ρ is the usual Lindblad dissipative superoperator, and κ 1 (κ 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ρ ss satisfying L 0ρss = 0. ( We briefly summarize prior work on this model. For Λ 3 ≡ 0, exact solutions forρ ss have been found using the positive-P approach [2][3][4]17]. The solutions express matrix elements ofρ 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 [18,19]; unlike the positive-P approach, these can also be used to describe dissipation at a non-zero temperature. Systems with a non-zero three-photon drive Λ 3 = 0 have not to our knowledge been studied before, nor were they 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. Λ 1 = κ 1 = 0), this system exhibits multiple steady states [10,[20][21][22][23]. We are not aware of any discussion of this using the positive-P approach. For ∆ = 0, the system is simple enough that the multiple steady states can be found via elementary means, in terms of superpositions of coherent states [11,12]. 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. [1] 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.

A. Recap of the basic approach
Consider first the case where our system in Eq. (2) has no two photon loss (κ 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 operatorb, system Hamilto-nianĤ 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 [24][25][26]. The dynamics of the reduced density matrixρ ab describing both cavities is described by a Lindblad master equation of the form: 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. [1], for a general cavity a Lindblad master equation having only single-photon loss (i.e. Eq. (2) with κ 2 = 0 and arbitraryĤ a ), one can always construct a HamiltonianĤ 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. Lettinĝ ρ ab,ss denote the steady-state density matrix of the twocavity problem, this means: Note that the dark state condition implies that |ψ is essentially a single mode state. Introducing new composite mode operatorsĉ one notes that the dark state condition forces the composite modeĉ − to be in vacuum. Hence, one just needs to solve for the (pure) state of the compositeĉ + mode. In physical terms, the CQA approach seeks to con-structĤ 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 |ψ .
While such a construction is always possible, in practice it would seem to be of no utility, as one can only construct the requiredĤ b by first independently solving for the cavity-a steadyρ a,ss , Despite this seeming obstacle, Ref. [1] demonstrated that for a range of problems, one could essentially guess the form ofĤ b without first knowingρ a,ss . This educated guess is extremely simple: H b is taken to be identical toĤ a up to an overall minus sign. Ref. [1] applied this to the simplest driven Kerr problem (Λ 2 = Λ 3 = κ 2 = 0 in Eq. (1)), in which casê 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. [2] for the linear-drive Kerr problem.

B. Extension to nonlinear driving and two-photon loss
It is natural to ask whether the 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 exctitations. 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: witĥ Again, tracing out cavity b from the above equation recovers the cavity a master equation given in Eq.
(2), independent of the choice ofĤ b . The next step is the same as before: we want to pick 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 |ψ 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 (â −b) |ψ = 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Ĥ b . As we show in Sec. IV, the simple educated guess of takingĤ b to be the negative ofĤ a still works in the presence of two photon driving and loss, and even for a wider class of problems.

C. 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 [13][14][15][16]. This provides an extremely efficient way of solving the complex recursion relations that determine the desired dark state wavefunction |ψ . 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.

Basics of the representation
Consider a single bosonic mode in a pure state |ψ that in is written in terms of Fock states |m as: In the SB representation, this state is associated with a holomorphic function ψ SB (z) defined on the complex plane: The space of these functions forms a Hilbert space that is unitarily equivalent to the original Fock space, with an 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 |ψ+ . 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 |ψ+ convolved with vacuum noise. As a result, cavity-a's steady-state Wigner function Wa(z) is equal (up to scaling) to the Q function Q ψ + (z) of the pure state |ψ+ .
induced inner product: The SB wavefunction has a direct physical interpretation: its modulus determines the Husimi Q-function of the state |ψ . Letting |z denote a coherent state with amplitude z, we have Finally, the canonicalĉ and creationĉ † operators become linear differential operators in the Bargmann space:ĉ → ∂/∂z,ĉ † → z.

Tracing out the auxiliary cavity
As we will see, the CQA method reduces to finding a single-mode, pure-state wavefunction |ψ + for the collective modeĉ + = (â +b)/ √ 2; the orthogonal modeĉ − must be in vacuum to have a dark state. To find the corresponding state of the physical cavity a, one transforms from theĉ ± basis to theâ/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 |ψ + 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 |ψ + 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 ψ +,SB (z) denote the SB representation of the pure state |ψ + , we have: We see that the SB "wavefunction" ψ +,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 |ψ + , 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.

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 positive-P methods [3,4,17,27], but in a manner that allows greater physical intuition. We are also able to solve an extended model which includes a three-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 (κ 1 = 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.

A. Solution without three-photon terms
We are interested in the driven-Kerr system described by Eqs. (1) and (2) with Λ 3 = 0 and κ 1 > 0. The CQA approach represents this system using the equivalent twocavity cascaded system in Eqs. (9) and (10). We seek a pure-state steady-state |ψ that is necessarily dark with respect to dissipation, meaning that Eq. (6) is satisfied: Our steady state can thus be written as a tensor product of a non-trivial state of thê c + collective mode, and a vacuum state for theĉ − mode: In order for |ψ to be a steady state, it also needs to be an eigenstate of the cascaded HamiltonianĤ ab with energy E. WritingĤ ab in terms ofĉ ± , and using the fact thatĉ − |ψ vanishes, the eigenvalue equation becomes witĥ Our choice of the auxiliary-cavity HamiltonianĤ b thus leads to a cascaded Hamiltonian that necessarily creates an excitation in theĉ − mode. It follows that we must have E = 0. Having |ψ be a stationary state then reduces to a single mode problem: i.e. we need to find a pure state |ψ + that is annihilated by the non-Hermitian operatorĤ + . The seemingly obvious next step is to follow the approach used in Ref. [1]: express |ψ + in the Fock state basis, and turn Eq. (19) into a recursion relation for the expansion coefficients α j : In special cases, this reduces to an easily solvable singleterm recursion relation: either the case of no two-photon driving Λ 2 = 0 [1], or the case of no one-photon driving Λ 1 = 0 [28]. In the more general case, the resulting twoterm recursion relation is more unwieldy. A more direct way of getting the desired solution is to use the SB representation ψ +,SB (z) of the state |ψ + . Eq. (19) is then transformed into a second-order ordinary differential equation: where Here, ∆ ≡ ∆ + iκ 1 /2 and K ≡ K − iκ 2 are, respectively, effective complex detuning and Kerr nonlinearity parameters. Without two-photon driving (i.e. λ 2 = 0), Eq. (21) is a standard hypergeometric equation. It has a unique analytic solution: 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. [1], which in turn corresponds to the classic solution of Ref. [2]. 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 positive-P solutions. Our expression for this case agrees with that derived earlier (via an alternate method) [29].
We turn now to the more interesting case where λ 2 = 0. Eq. (20) 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" where θ(z) is the "gauge potential". This transformation shifts the differentiation operator by the gradient of Here, we try the simplest potential θ(z) ≡ z, with some constant. Note that as θ 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: After our transformation, the problematic two-photon driving term is effectively shifted by an amount 2 : It can thus be eliminated by choosing such that We will call these non-unitary gauges plus-gauge and minus-gauge. Choosing, e.g. the plus gauge ≡ + , we see that the gauge-transformed state φ(z) satisfies Kummer's differential equation (see [30]), so that: 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. [31]). 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 → 0 limit, the solution above tends smoothly to the Bessel-function solution in Eq. (23). 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: where N is a normalization constant. Note that if φ(z) = 1, then W a,ss (z) corresponds to a coherent state with amplitude α = −λ 2 /2. Thus, a non-unity φ(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. (24), one obtains an identical solution (see Appendix D).

B. Including the three-photon drive term
We now allow Λ 3 = 0 in Eq. (1). We are still able to exactly solve for the steady state in this case; unless κ 2 = 0, it has a qualitatively different form from the Λ 3 = 0 case. The CQA method proceeds as in Sec. IV A. We again write the two-mode dark state as |ψ = |ψ + |0 − , and the eigenvalue equation again reduces to finding the kernel of a non-Hermitian operatorĤ + : Comparing against Eq. (18), we see that the presence of Λ 3 creates a term proportional toĉ 2 + . Attempting to solve directly for |ψ + in the Fock basis leads a complicated recursion relation, as now we have terms that add a photon (∝ĉ † + ), as well as those that subtract two photons (∝ĉ 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, where 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: Here, ± correspond to the non-unitary gauge choices in which the displaced two-photon drive vanishes (c.f. Eq. (27)): To be manifestly consistent with the solution of the driven Kerr cavity without three-photon terms, we have again written the solution in the plus gauge. Finally, λ 3 , λ 2 , λ 1 , D are the following general complex constants: We have again defined∆ For the case where Λ 3 → 0, these parameters revert to those given before Eq. (22). Note that for vanishing twophoton loss, K is real, and hence Eq. (36) implies that λ 3 = 0. In this case, the three photon drive does not give us anything qualitatively new, as it can be completely eliminated by our displacement transformation. In contrast, for non-zero κ 2 , three photon driving gives rise to genuinely new phenomena.
The most significant consequence of the three photon drive is that it can be used to tune the parameter λ 2 to zero without needing to make an additional non-unitary transformation as in Eq. (24). Indeed, if we pick Λ 3 to satisfy then λ 2 vanishes (c.f. Eq. (36)). One can then take = 0 in Eq. (33), implying that the displaced dark state |ξ + is directly given by a hypergeometric function. While this might seem like a technical nicety, it has pronounced observable consequences; this is discussed in detail in Sec. V C and App. B.

V. PHASE DIAGRAM OF NEW PHENOMENA
We now use our exact solutions in Eqs. (28) and (32) 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 Λ 3 = 0 (no three photon drive), these are: The various drive amplitudes Λ 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 Λ 3 , one has r 1 = ( 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 Λ 1 = ∆ = κ 1 = 0 (i.e. r 1 = r 2 = 0) [10]. 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 Λ 3 = 0. We also highlight the fact that with a nonzero three-photon drive, the observable consequences of the photon blockade and anti-blockade phenomena can be made even more dramatic.

A. Basic intuition
Recall that the steady state is determined by a singlemode pure-state |ψ + (c.f. Fig. (2)), and that further, this state is related to a simpler state |φ via a "non-unitary gauge transformation" (c.f. Eq. (25)). We could always expand the transformed state |φ in the Fock basis as: Defining the scaled Fock state amplitudes (c.f. Eq. (41))  . This termination corresponds to a kind of quantum interference effect, and will be at the heart of the new blockade, antiblockade and bistability phenomena we describe.
Note that we can directly go from the Fock state structure of |φ to the SB wavefunction of the desired, untransformed state |ψ + . For Λ 3 = 0 the SB wavefunction of |ψ + is For the more general case with non-zero Λ 3 , up to a displacement, we have: where ± are defined in Eq. (33). Recall that these SB wavefunctions directly determine the steady-state Wigner function of the physical cavity via Eq. (15).
B. Pure unique steady states: r1 = r2 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, linearlydriven, 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 Λ 3 = 0, this requires tuning the one and two photon drives Λ 1 , Λ 2 so that: 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 the state |ψ + is a coherent state with amplitude 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 γ/ √ 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.
C. Higher-order photon blockade: r1 = n0 Surprising effects also occur when drives and detuning are chosen so that r 1 = n 0 , where n 0 is a nonnegative integer. The recursion relation in Eq. (43) now terminates at m = n 0 : Fock state amplitudes c m vanish for all m ≥ n 0 + 1. This is an example of a generalized strong photon-blockade phenomena: the gaugetransformed steady-state |φ has strictly zero probability to have more than n 0 photons. Unlike standard photon blockade [32], the mechanism here does not require infinitely strong nonlinearity. Also, unlike the so-called "unconventional" photon blockade [33][34][35], the blockade here is complete: there is strictly no probability to have more than n 0 photons in the state.
While the "gauge-transformed" state |φ exhibits blockade, physical phenomena is controlled by the untransformed state |ψ + . Eq. (25) 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 blockadeinduced suppression can be observed by considering how the steady state changes as a function of the single photon drive amplitude Λ 1 (as this tunes r 1 but not r 2 ). From Eq. (39), one sees that blockade occurs periodically as a function of Λ 1 , with the nth-order blockade occuring when where Λ Fig. 4 shows representative results for κ 2 = Λ 3 = 0: the average cavity photon number shows a sharp suppression whenever Λ 1 is tuned to make r 1 a positive integer. We stress that this is a purely quantum phenomena: the semiclassical solution does not exhibit any such oscillations. We also stress that the width of these blockade suppressions (as a function of Λ 1 ) are much smaller than κ 1 .
The results shown in Fig. 4 are for the usual case of no three-photon drive, i.e. Λ 3 = 0. The situation becomes even more interesting with non-zero Λ 3 . By tuning this drive, the generalized photon blockades associated with r 1 = n 0 can be made sharp: after a simple phase-space displacement, the cavity-a steady state has strictly zero probability for having more than n 0 photons in it. This is discussed in more detail in App. B.
Finally, we note that when r 1 = n 0 , the SB wavefunction for the dark state |ψ + (which directly determines the cavity-a Wigner function) reduces to an associated Laguerre polynomial L (α) m (z): 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) [3]. 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 m0 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 photonblockade phenomena, this will also have implications for our physical cavity, via Eq. (25).
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 κ 1 , κ 2 ≥ 0). This remains true even in the presence of a three-photon drive, where r 2 = D, with D given by Eq. (34). Nonetheless, for weak dissipation (namely κ 1 ∆, κ 2 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. 5. Note that this resonance phenomena was observed in Ref. [3], though connections to photon statistics and the properties of the analytic steady-state solution were not discussed.

E. 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 Λ 3 = 0, we have Without the exponential prefactor, this state would exhibit both photon blockade and antiblockade (i.e. the its photon number distribution is cut-off at small and large photon numbers). 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 ≥ 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: 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: This also implies multi-stability for the physical a cavity, which exhibits a two-parameter continuous family of steady states, 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 antiblockade 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 κ 1 , κ 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 κ 1 = 0, κ 2 > 0. This parameter point corresponds to the well-studied cat-state bistability in a two-photon driven Kerr resonator [21].
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 A 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).

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 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. (53).
Suppose we chose parameters that result in (r 1 , r 2 ) being close to integers (n 1 , n 2 ): These small deviations kill the bistability. However, for small δr j the resulting pure steady state of the + mode is a particular linear combination of the states φ j (z) that span the bistable manifold at δ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 In the regime that δr 1 , δr 2 1, we can see that the ratio c m+1 /c m is essentially 1, except for the ratio c n+1 /c n = δr 1 /δr 2 . Therefore, as δr 1 , δr 2 → 0, the unique steady state solution (i.e. solution to the recursion relation) has the limiting form as a superposition of the bistable solutions given in Eqs. (50),(51). Note that in writing this equation, we must pick the overall phase of ψ 2,SB such that the ratio between c n+1 (appearing in ψ 2,SB ) and c n (appearing in ψ 1,SB ) is precisely δr 1 /δr 2 .
As a result, the unique steady state Wigner function of the physical a cavity will be: 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.
A. Cat-state bistability: (r1, r2) = (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: Λ 1 = ∆ = κ 1 = 0. This corresponds to the well-known quantum bistability that occurs in a two-photon driven Kerr resonator [9,11,12], a system where photon number parity is conserved. The two distinct solutions to the recurrence relation in Eq. (43) are c j = δ j,0 and c j = 1−δ j,0 (c.f. Eqs. (50)-(51)). This corresponds to two distinct dark states for the + mode, with SB wavefunctions ψ 1,SB (z) corresponds to a coherent state with amplitude ≡ i √ λ 2 , whereas ψ 2,SB (z) corresponds to an odd cat state (odd superposition of coherent states with amplitude ). Note that we have picked the global phase of ψ 2,SB to be compatible with Eq. (57).
We thus have a direct connection between this paritybased bistability and the photon blockade and antiblockade 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 δr 1 = δr 2 = 0: For small imperfections, we can then use Eq. (57) to give us the steady-state SB wavefunction: where N is a normalization constant, and Eq. (62) directly gives us the Wigner function of the unique steady state via Eq. (15). Each term in Eq. (62) on its own corresponds to a simple coherent state (amplitudes ± ≡ ±i √ λ 2 ); thus, in general we have a steady state that is an incoherent mixture of two coherent states. This equation also reveals something surprising: the localization of the steady state in phase space is a nonmonotonic function of Λ 1 . The state is delocalized both for Λ 1 = 0, and for Λ 1 large enough to make Q 1. Representative results are shown in Figure 6. B. Quantum bistability with a single photon drive: (r1, r2) = (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 ) (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 = n k=1 δ j,k and c j = n k=1 (1 − δ j,k ) (c.f. Eq. (50-51)). These in turn correspond to two distinct +-mode states where Γ(r, z) ≡ ∞ z 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 ∆ = n/2K and a single photon drive Λ 1 = −i(n/2) √ Λ 2 K. If we now include single photon loss (but keep κ 2 = 0), and also shift Λ 1 slightly from the above value, the recurrence parameters are slightly shifted as well: Hence, via Eq. (57), 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 Λ 1 . Note that by picking parameters so that δr 1 = δr 2 , the steady state becomes a coherent state with amplitude γ = −Λ 2 /K, whereas if δr 1 = 0, it has a bimodal form.

C. 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 nonzero eigenvalues of the Liouvillian 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 slowmode from other decay modes. Formally, if we let γ j denote the decay modes of the Liouvillian (i.e. negative real parts of the eigenvalues of L 0 ), and order rates such that γ 1 ≤ γ 2 ≤ ...., then in near-bistable regimes: Note that this hierarchy of dissipative rates has already been described for the more familiar (r 1 , r 2 ) = (0, 0) "catstate" bistable point [12]; 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. 7(a).
For near-bistable parameters, the CQA approach provides insight into the nature of the slow decay mode of L 0 . As one might expect, this mode corresponds to slow relaxation within the bistable manifold of states.
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. (53) (and incoherent mixtures of these states). Density matrices in this bistable manifold lie in  Slow dynamics near generalized bistable regimes. (a) Solid line: ratio of the two smallest relaxation rates (i.e. dissipative rates of the system Liouvillian L0), as a function of κ1. Dashed line: κ1/γ1. κ1 → 0 corresponds to being at a bistable parameter point (r1, r2) = (n, n), either n = 2 (light green) or n = 4 (dark green). Parameters are Λ2 = 6K, Λ3 = κ2 = 0, ∆ = nK/2, and Λ1 = −in √ Λ2/2. One sees that the slow rate γ1 is much slower than κ1, and that there is a pronounced dissipative gap. (b) Solid line: The measure 1 − P (c.f. Eq. (71)) of how closely the slowest system decay mode (with rate γ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 span of the four operators (i, j = 1, 2): By Appendix C, these operators have Wigner transforms with N a normalization constant.
The slow mode (rate γ 1 ) has an associated right eigen-vectorM slow , i.e. L 0Mslow = −γ 1Mslow . If the slow dynamics is entirely in the bistable manifold, thenM slow should lie completely within the span of theM 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: Here M ij is an orthonormal basis for the span ofM ij (obtained via the standard Gram-Schmidt process). As 0 ≤ P ≤ 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. 7(b). One sees that for small κ 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 takinĝ with α j the classical amplitudes, determined by (with K ≡ 1): One sees from Fig. 7(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.

VII. PARITY-CONSERVING DYNAMICS: TRUE QUANTUM BISTABILITY
We now focus on a special case that has received considerable recent attention [9][10][11][12]: a system where κ 1 = Λ 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 [36]. 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. More generally, the CQA method provides a transparent understanding of this quantum bistability that is difficult to achieve using the positive-P solution method.
A. The CQA approach in the limit of vanishing single-photon loss We start by revisiting the CQA method of Sec. III for systems described by Eq. (2) with κ 1 = Λ 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 κ 1 = 0, we only have the a two-photon loss dissipator, given by where again the collectiveĉ ± modes are defined in Eq. (7). There are now two distinct possibilities for a non-trivial dark state: either theĉ − mode is forced to be in vacuum (with the + mode occupied), or theĉ + mode is forced to be in vacuum (with the − mode occupied). The first option is the same as what we did for κ 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: 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 |ψ dk can have both + and − modes occupied.

B. Zero detuning: quantum bistability
Consider first the case ∆ = 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 [11,12]. Including two-photon loss, the dissipative system exhibits true quantum bistability: it admits multiple steady states. 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. (75). To have this state be a steady state, it must also be an eigenstate of the cascaded Hamiltonian (c.f. Eq. (10)). For ∆ = 0, this leads to the equation: with λ 2 = 2Λ 2 /(K − iκ 2 ). SinceĤ casc always adds an excitation to both modesĉ − andĉ + , the only possible energy eigenvalue is E = 0. The equation then decouples into separate equations for |ψ + , |θ − 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: where the coherent state amplitude = i √ λ 2 , as in the previous section. We see that the cascaded twocavity 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. (77), 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: with˜ ≡ / √ 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 µ − = ν + = 0.
The upshot is that we have a steady state manifold for cavity a that is two dimensional, and spanned by the states |±˜ a (in agreement with previous work [10]). In simple terms, the steady-state manifold corresponds to a quantum bit, i.e. a full single-qubit Bloch sphere [36]. This is what we mean by the system exhibiting quantum bistability.

Loss of quantum bistability
We next consider the case of a non-zero detuning ∆. While parity continues to be conserved, this addition has a dramatic consequence for the nature of the bistability. Solving the system again using the CQA method, the requirement of having our general dark state in Eq. (75) be a energy eigenstate of the cascaded Hamiltonian leads to the equationŝ where D = 2∆/(K − iκ 2 ). As before, the equations determining |ψ + and |θ − 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: With the same equation for θ SB , and with λ 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 ≡ 0, the ODE above has no singular points, and thus the standard existence theorem ( §12.22 in [37]) guarantees two independent, analytic solutions. As discussed earlier, this leads to quantum bistability for the physical mode a. However, the term ∝ D∂ 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: 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 |ψ + and |θ − , the most general dark state has the form of Eq. (75), and corresponds to a twodimensional subspace. In what follows, it will be useful to write this general dark state as with 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 steady-state manifold of cavity a.
This procedure leads us to consider four linearlyindependent cavity-a operators (that determine the cavity a density matrix after tracing out cavity b): To understand the structure of these operators, we consider their corresponding Q-functions (easily obtainable using the SB representation): We obtain an important result: these four operators are not all independent. Because of the symmetry of each integral under the mapping u → −u, we havê 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ρ + a,ss ≡M ++ a,sŝ ρ − a,ss ≡M +− a,ss .
withρ ± a,ss 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 where the extremal statesρ 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 ≤ p ≤ 1, which simply corresponds to the dynamicallyconserved probability of having an even photon number parity. In simpler terms, the cavity-a steady-state manifold corresponds to a classical bit [36].

Full characterization of the steady state manifold
To conclude our discussion of ∆ = 0, we use the CQA method to compute exactly each steady state in the bistable manifold. We begin by noting that the states |Φ e/o in Eq. (85) 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ρ e/o :ρ To compute these steady-states, we note that by substituting Eq. (85) into Eq. (92) we can expand, e.g.
where N is just the normalization constant N for the dark state |ψ + , which has the exact expression Inverting the above linear relation, we get ρ + a,ss = 1 2 This equation immediately leads to exact expressions for ρ e/o , which are given in Appendix G.
Incidentally, Eq. (95) tells us which state is selected from the bistable manifold in the presence of an infinitessimal amount of single-photon loss (or any other symmetry-breaking perturbation, e.g. one-or threephoton driving). Indeed, we can rewrite it aŝ ρ a,ss ∼ κ1,Λ1,Λ3→0 We will briefly discuss physical intuition from this formula, as well as the closed form of N . In the weakdriving limit λ 2 → 0, the hypergeometric series defining N collapses to just the first term, and we get N −1 → 1, soρ a,ss ∼ Λ2→0ρ e .
In contrast, in the strong-driving limit N diverges, and thusρ a,ss A final piece of physical intuition: since N is a function only of the modulus |λ 2 |, the relative bias (towards eitherρ o/e ) is independent of the phase φ of the drive λ 2 ≡ e iφ |λ 2 |.

VIII. CONCLUSIONS
In this, work, we have presented a generalization of the coherent quantum absorber method developed by Stannigel et. al. [1] 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 photonblockade 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. Figure 8. 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.
the chirality of one of the waveguides in the absorber setup (see Fig. 8). 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 tô Again, using the absorber method, we can solve this master equation in a manner identical to before, i.e. with |ψ = |ψ + |0 − , except now we have κ In this case Eq. (56) actually becomes a relation for selecting a pure state in the bistable manifold: |ψ + = δr 2 |ψ +,1 + δr 1 |ψ +,2 . (A2) In this case, |ψ +,j 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 Λ 3 = 0, where |z 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 = m), one stabilizes even more exotic states, whose sum may no longer be Gaussian. Appendix B: Using three photon driving to achieve sharp photon blockade and anti-blockade As discussed in Sec. V C, by tuning system parameters, the "gauge-transformed" state |φ (c.f. Eq. (25)) can exhibit sharp cutoffs in its photon number distribution (i.e. photon blockades and anti-blockades). The physical cavity state is however controlled by the untransformed state |ψ + . For = 0, the exponential factor in Eq. 25 tends to smear these blockades / anti-blockades. Remarkably, if one has a non-zero three photon drive and properly tunes its amplitude Λ 3 to satisfy Eq. (38), then λ 2 = 0, implying that the states |φ and |ψ + coincide. For these parameters, there is thus no smearing of phenomena, and the sharp generalized photon blockade exhibited by |φ manifests itself more directly in the physical cavity-a steady state. To be explicit, in this case, the physical cavity-a steady state exhibits a sharp photon blockade (sharp cut-off in its photon number distribution) after a simple phase-space displacement. We discuss this regime further in what follows; closed-form expressions for photon statistics and cavity moments are given in Appendix F.
We start by noting that for Λ 3 = 0, the condition of Eq. (38) for eliminating the smearing parameter reduces to Λ 2 = 0; in this case, we are left with the simplest version of the Kerr system, whose Bessel-function solution (c.f. Eq. (23) does not support any unusual photon statistics. In contrast, if instead satisfies Eq. (38) with nonzero Λ 3 and Λ 2 , we then get something more non-trivial. In particular, we can still tune the parameter r 1 = n 0 so that the displaced state |ξ + (c.f. Eq. (31)) exhibits photon blockade. In this limit, the displaced dark state reduces to an associated Laguerre polynomial: Note that the condition for r 1 = n 0 implies the quantization condition λ 1 = −nλ 3 , and so, at these blockaded points, Λ 1 is also fixed in terms of the strength Λ 3 of the cubic drive: with the first-order blockade occurring at Recall that the state Eq. (B1) determines the physical cavity-a state through a set of two operations: a phasespace displacement as per Eq. (31), and the mixing with vacuum at a 50/50 beamsplitter. This beamsplitter operation is of course unable to add photons to the original state. The displacement operation can be effectively undone by performing a final displacement of the cavity-a state by an amount After this simple displacement operation, the cavity-a steady state will exhibit a sharp cutoff in its photon number distribution, i.e. photon blockade. Phrased equivalently, the cavity-a steady state (for these parameter tunings) is a displaced version of a truly blockaded state. This cutoff in the photon number distribution is also manifest in the exact expression Eq. (F5), which gives us the cavity photon number distribution: where ξ n ≡ ∂ n ξ +,SB (0) is the n-th derivative of the dark state in Eq. (B1), evaluated at the origin z ≡ 0 in phase space. Therefore, if the dark state is algebraic, i.e. a polynomial of some finite degree d, then P (n) vanishes for all n > d. Representative results are given in Figure  9. Eq. (38) and Eq. (B2) together provide the fine-tuned values of Λ 1,2 (as a function of Λ 3 ) needed to stabilize this exotic blockaded state. In summary, we see that the three-photon term adds novel features to the steady-state phase diagram of the Kerr cavity.
Appendix C: Steady-state Wigner function We will now rigorously prove the connection between the Wigner function of the steady-stateρ a,ss of the driven Kerr cavity, and the modulus squared of the SB representation of its purification |ψ + . 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 [38], and we will review the salient results here. Specifically, given a (possibly non-Hermitian operator)Â of a quantum-mechanical mode, define its normally-ordered symbol σ N to be where :Â : is the operatorÂ 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: where :Â : S is the operatorÂ but re-expressed according to the symmetric ordering convention (as defined in [38]). The symmetrically-ordered symbol is proportional to the standard Wigner transform, which can be formally computed via an integral: For positive semi-definite operators, e.g. a density matrix ρ ≡ρ † , the symmetrically-ordered and normally-ordered symbols coincide with the Wigner-and Q-functions respectively: What Glauber and Cahill showed in [38] is that operator symbols corresponding to different ordering conventions are related by the heat semigroup. In particular, we have the following theorem: Theorem ( [38]). LetÂ be a Hilbert-Schmidt operator (i.e. Tr[Â †Â ] < ∞). 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/4, i.e.
where K(t, z, z ) is the heat kernel on C, which can be exactly computed and comes out to 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 |ψ 1 , · · · , |ψ k of the space of dark states, i.e.ĉ − |ψ j = 0.
Clearly, these states form a k-dimensional Bloch sphere, spanned by their outer-products: According to CQA, to obtain the corresponding stationary modesM ij of the physical cavity-a Lindbladian, we must trace-out the ancilla mode: 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ĉ ± as |ψ j = |ψ j,+ |0 − : Letting σ +,ij (z) denote the symmetrically-ordered symbol (i.e. Wigner transform) of the outer-product |ψ i,+ ψ j,+ |, and σ − (z) denote the symmetricallyordered symbol of the vacuum state |0 − 0 − |, let denote the symmetrically-ordered symbol (i.e. the Wigner transform) of the stationary modeM 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 σ − (z) then acts as a Gaussian filter for the symbol σ +,ij (z): where we have defined symmetrized and antisymmetrized phase-space variables z ± ≡ (z ± z )/ √ 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 → √ 2z. Indeed, we can rewrite it in terms of the heat kernel: where σ N +,ij is the normally-ordered symbol of the mode |ψ i,+ ψ 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: where |z denotes a coherent state with amplitude z. By expandingÂ in terms of outer-products aŝ and utilizing the property ψ SB (z) = z * |ψ e −|z| 2 /2 , Bargmann in [16] was able to show that this implies By substituting the exact expression for the normallyordered symbol into Eq. (C15), we can finally state the main result utilized in the main text: 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 where α ij = α * 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 {α ij }, and letting {φ j (z)} denote the Segal-Bargmann representations of the corresponding eigenvectors, the Wigner function can be equivalently written as where normalization forces j p j = 1.
As a simple example, when the dark-state subspace is one-dimensional, the Wigner function is the squaredmodulus of the SB representation of the unique, normalized dark state in that subspace: 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.
However, the results are gauge-invariant, as we can write where ψ +,SB is the dark state Eq. (32) in the main text, and in the first line (Eq. (D2)) we utilized Kummer's transformation (see e.g. Ref. [30]), which is a fundamental symmetry of the confluent hypergeometric differential equation: a. Steady-state density matrix In terms of these Taylor coefficients, the steady state density matrix can be computed in the Fock basis: Using identities of the form 0|ĉ m+l + (ĉ † + ) j |0 = δ m+l,j j!, etc., we get the remarkably simple result: This expression matches similar expressions obtained using positive-P solutions, as we will see later in this section.

b. Cavity moments
We can also express the normally-ordered moments of a driven Kerr cavity exactly in terms of the scaled Fockstate amplitudes ψ l . The calculation is slightly more straightforward: Tr[ρ a,ss (a † ) n a m ] = ψ|(a † ) n a m |ψ where |ψ , as before, is the purification of the density matrix obtained from the absorber method. Expanding the dark state yields Tr[ρ a,ss (a † ) n a m ] Defining a new variable l such that j ≡ n + l, we find that k = m + l, and that furthermore l ≥ 0. So our sum simplifies to 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 [3], using Pfunction methods.
Utilizing the identity (z) m+l = (z) m (z + m) l , the sum closes, and we get m|D αρa,ssD † α |n = ξ m ξ * n √ 2 m+n n!m! · 2 F 2 (n − r 1 , m − r * 1 ; m − r 2 ; n − r * 2 ; |λ 3 | 2 /2) (F5) 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: Appendix G: 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 E on unique steady states: m|ρ a,ss |n = 1 √ 2 m+n m!n! ∞ l=0 ψ m+l ψ * n+l 2 l l! (G1) note that, as a direct consequence of ψ 2l−1 ≡ 0, we have 2j + 1|ρ a,ss |2k = 2j|ρ a,ss |2k + 1 = 0, as each term in the sum over l would identically vanish in these cases. In summary, Π eρa,ssΠo =Π oρa,ssΠe = 0, Therefore, to compute the steady statesρ e/o , it suffices to compute matrix elements ofρ a,ss . We note that this was done in [3] (as this represents the unique steadystate regime κ 1 = 0), and so we're technically done, as we could simply cite the result here.
For completeness, however, we show that the calculation of the expressions on the RHS's of Eqs. (G4-G5) can be reproduced in a straightforward manner within the quantum absorber formalism. Assuming m ≡ 2j, n ≡ 2k are both even, we have Therefore, in total we have Assuming m ≡ 2j + 1, n ≡ 2k + 1 are both odd, we have Therefore, in total we have ∞ l=0 ψ m+l ψ * n+l 2 l l! = ψ 2j ψ * 2k 2 (j + 1 2 )(k + 1 2 )|λ 2 | 2 (j − r 2 )(k − r * 2 ) ; |λ 2 /4| 2 . (G10) In summary, we have the following closed-form forρ a,ss in the Fock basis: For small detuning |D| 1, the CQA solutions above approach even/odd cat states, both of which exhibit Wigner function negativity. In the large detuning limit |D| 1,ρ e/o also approach pure states: the vacuum state and one-photon state respectively, one of which exhibits Wigner function negativity. The exact CQA solutions are validated against master equation numerics in Figure 10.
For completeness, we include here the calculation of N in Eq. (G15) (this expression also shows up in the main text, in Eq. (96), where it controls the average parity of the unique steady state selected when parity-symmetry is spontaneously broken). The derivatives of the dark state, evaluated at the origin in phase space, are whereas the odd derivatives vanish at the origin (ψ 2l−1 = 0). However, note the following identity (2l)! l! = 2 2l ( 1 2 ) l . With this identity, the Fock state amplitudes take the simpler form: with r 1 ≡ −1/2, and r 2 ≡ r 1 + D/2. Having computed the Taylor coefficients ψ l , Eq. (E10) gives us the normalization: (G18) This concludes the calculation of the normalization constants in the expressions Eqs. (G4-G5).