# Supplementary information for:
Driven-dissipative phases and dynamics in nonlinear non-Markovian photonics

Jamison Sloan1,†, Nicholas Rivera2,3,†, and Marin Soljačić1,3 1Research
Laboratory of Electronics, Massachusetts Institute of Technology, Cambridge,
MA 02139, USA.
2Department of Physics, Harvard University, Cambridge, MA 02138, USA.
3Department of Physics, Massachusetts Institute of Technology, Cambridge, MA
02139, USA

###### Abstract

In this Supplementary Information (S.I.), we derive the results of the main
text, related to dispersive coupled resonances in the presence of intensity-
dependent nonlinearity. We use the input-output formalism to show how these
systems can be used to generate quantum states of light when coherently
pumped. We also perform stability analysis to show how these systems can
exhibit modulational instability due to non-Markovian dissipation.

###### Contents

  1. S1 General formalism for non-Markovian dissipation
    1. A Heisenberg-Langevin equations of motion
    2. B Kramers-Kronig relations and non-Markovian correlation functions
    3. C Preservation of commutation relations
  2. S2 Steady-state mean-field behavior
  3. S3 Steady-state quantum noise properties
    1. A General formalism for quantum noise
    2. B Analytical approximations for quantum noise in the adiabatic limit
    3. C Additional results for for generation of highly intensity-squeezed states
  4. S4 non-Markovian loss models
    1. A Friedrich-Wintgen model
    2. B Fano mirror model
  5. S5 Modulational instabilities and phase diagrams

##  S1 General formalism for non-Markovian dissipation

We begin by developing the general formalism for the description of a
nonlinear cavity mode with non-Markovian loss rivera2023creating . Such a non-
Markovian loss arises from coupling to frequency-dependent reservoirs, and
result in frequency-dependent dissipation. In the presence of intensity-
dependent (Kerr) nonlinearity, such frequency-dependent dissipation can lead
to an effective intensity-dependent dissipation.

###  A Heisenberg-Langevin equations of motion

The equation of motion for the annihilation operator of a Kerr nonlinear
resonance $a$ coupled to non-Markovian dissipation channels $i$ is:

$\dot{a}=-i\omega_{a}a-i\beta(a^{\dagger}a)a-\int
dt^{\prime}\,K_{l}(t-t^{\prime})a(t^{\prime})+\sum_{i}\int
dt^{\prime}\,K_{c,i}(t-t^{\prime})s_{i}(t^{\prime}).$ (S1)

Here, $\omega_{a}$ is the bare frequency of the resonance $a$, $\beta$ is the
Kerr frequency shift per photon, $K_{l}(t-t^{\prime})$ is a time-domain loss
Kernel which provides frequency-dependent damping of $a$ through a
convolution, and $K_{c,i}(t-t^{\prime})$ describe the coupling of external
input fields $s_{i}(t)$ in through various channels $i$. Each $i$ can be used
to describe the coupling of the mode $a$ to a different reservoir. In the
frequency-domain, the coupling and loss kernels are associated with frequency-
dependent functions

$K_{c,l}(\omega)=\int dt\,e^{i\omega t}K_{c,l}(t).$ (S2)

These functions can be constructed through theoretical modeling (see Section
S4), or empirically through experimental data.

Additionally, we note that the input fields are normalized so that
$[s_{i}(t),s_{j}^{\dagger}(t^{\prime})]=\delta_{ij}\delta(t-t^{\prime})$,
which means that $|s_{i}(t)|^{2}$ has units of photon flux in photons/s. Going
forward, we will sometimes omit the reservoir index $i$ from the input fields
and coupling kernels if it is assumed that there is only one input channel.

Markovian limit: We briefly describe how the theory reduces to the more
typically used Markovian limit. Physically, the Markovian limit corresponds to
situations where the bandwidth of the reservoir is large compared to that of
the mode under consideration, resulting in a so-called “white noise” coupling.
In the time-domain, this is equivalent to the condition that the cavity mode
and reservoir interact instantaneously: any internal dynamics of the coupling
are fast compared to all other relevant timescales. In this case, the cavity
is described by a frequency independent loss rate $\gamma$. The loss rate
determines the coupling and loss functions as $K_{l}(\omega)=\gamma$, and
$K_{c}(\omega)=\sqrt{2\gamma}$. In the time-domain, the loss and coupling
functions are delta functions (indicating the instantaneous response), so that
for a single Markovian coupling channel, one recovers the equation of motion:

$\dot{a}=-i\omega_{a}a-\gamma a+\sqrt{2\gamma}s(t).$ (S3)

The quantum fluctuations of the source $s$ generate the Langevin force
$F_{\gamma}(t)=\sqrt{2\gamma}\delta s(t)$. In particular, we have decomposed
the input field as $s(t)=\braket{s(t)}+\delta s(t)$. For a vacuum reservoir,
$\braket{F_{\gamma}(t)F_{\gamma}^{\dagger}(t^{\prime})}=2\gamma\delta(t-t^{\prime})$,
while all other quadratic expectation values of $F_{\gamma}$ and
$F_{\gamma}^{\dagger}$ vanish.

Classical limit: In the classical limit, defined by the expectation value
$\braket{\cdot}$ which converts operators to c-numbers, Eq. S1 reduces to a
temporal coupled mode theory (TCMT) equation which includes the time delay
effects needed to describe memory effects. In this case, we obtain the
classical non-Markovian TCMT equation:

$\dot{\alpha}=-i\omega_{a}\alpha-i\beta|\alpha|^{2}\alpha-\int
dt^{\prime}\,K_{l}(t-t^{\prime})\alpha(t^{\prime})+\sum_{i}\int
dt^{\prime}\,K_{c,i}(t-t^{\prime})\braket{s_{i}(t^{\prime})}.$ (S4)

Here, we have used the notation $\alpha=\braket{a}$, and the sources
$\braket{s_{i}(t)}$ are now c-number classical sources. Due to the
nonlinearity, this TCMT equation makes the additional approximation that
$\braket{a^{\dagger}aa}\approx\braket{a^{\dagger}}\braket{a}\braket{a}=|\alpha|^{2}\alpha$.
Such an approximation provides an accurate description of systems with quantum
states which are localized in phase space, or in other words, states which are
well-described as a coherent part, with small added quantum fluctuations. An
equation such as Eq. S4 admits numerical solution with split-step methods, in
which the Kerr nonlinearity is treated in the time domain (due to the time-
local nature of the term), and the non-Markovian terms are treated in the
frequency domain (where the convolutions become products such as
$K_{l}(\omega)a(\omega)$).

If the Markovian and classical limits are taken together, we recover the
standards temporal coupled mode theory expression for a leaky resonance:

$\dot{\alpha}=-i\omega_{a}\alpha-\gamma\alpha+\sqrt{2\gamma}\braket{s(t)}.$
(S5)

###  B Kramers-Kronig relations and non-Markovian correlation functions

The coupling and loss kernels are not independent from one another. This is a
manifestation of the fact that the frequency-dependent damping rate arises
from the frequency-dependent coupling of the resonance $a$ to the external
reservoirs. In fact, the knowledge of one of these functions is sufficient to
determine the other. The loss and coupling functions are related in the
frequency-domain by a Kramers-Kronig relation:

$K_{l}(\omega)=i\int\frac{d\omega^{\prime}}{2\pi}\frac{|K_{c}(\omega^{\prime})|^{2}}{\omega-\omega^{\prime}+i\eta},$
(S6)

with $\eta$ infinitesimal. This is consistent with a non-Markovian
generalization of the so-called “Einstein relations” which link drift and
diffusion terms in Heisenberg-Lanegvin equations in Markovian settings
scully1999quantum .

We will describe some additional properties of the non-Markovian correlation
functions. To do so, we will consider Eq. S1 without any Kerr nonlinearity
$(\beta=0)$. Additionally, we will assume that the input field $s(t)$ has no
mean value, and consists only of operator-valued vacuum fluctuation $\delta
s(t)$. Physically, this corresponds to the case where light is allowed to leak
out of the cavity from an initial state, without additional energy input. Such
an assumption allows us to isolate the action of the non-Markovian reservoir
on the cavity, and enables us to write the equation of motion as:

$\dot{a}=-i\omega_{a}a-\int
dt^{\prime}\,K_{l}(t-t^{\prime})a(t^{\prime})+F(t).$ (S7)

In this expression, the final term is a Langevin force defined as

$F(t)\equiv\int dt^{\prime}\,K_{c}(t-t^{\prime})\delta s(t^{\prime}).$ (S8)

In the time domain, the correlator of the Langevin force is given as

$\braket{F(t)F^{\dagger}(t^{\prime})}=\int
dt^{\prime\prime}\,K_{c}(t-t^{\prime\prime})K_{c}^{*}(t^{\prime}-t^{\prime\prime}).$
(S9)

We thus immediately see that this relation will in general not be a delta
function.

This result can also be used to obtain correlations in the frequency domain,
which have a much more clear interpretation. In the frequency domain, this
input field commutation relations are
$[s(\omega),s^{\dagger}(\omega^{\prime})]=2\pi\delta(\omega-\omega^{\prime})$.
Using this, as well as the time domain correlator for $F$ shown above, the
nonzero correlator of the frequency domain Langevin force $F(\omega)$ can be
evaluated as

$\displaystyle\braket{F(\omega)F^{\dagger}(\omega^{\prime})}$
$\displaystyle=\int dt\,dt^{\prime}\,e^{i\omega
t}e^{-i\omega^{\prime}t^{\prime}}\braket{F(t)F^{\dagger}(t^{\prime})}$ (S10)
$\displaystyle=\int dt\,dt^{\prime}\,dt^{\prime\prime}\,e^{i\omega
t}e^{-i\omega^{\prime}t^{\prime}}K_{c}(t-t^{\prime\prime})K_{c}^{*}(t^{\prime}-t^{\prime\prime})$
(S11) $\displaystyle=2\pi|K_{c}(\omega)|^{2}\delta(\omega-\omega^{\prime}).$
(S12)

As noted earlier, the modulus squared of the coupling function is proportional
to the frequency-dependent loss $K_{l}^{\prime}(\omega)$, and thus the
correlator can also be written as
$\braket{F(\omega)F^{\dagger}(\omega^{\prime})}=2\,\text{Re}\,K_{l}(\omega)\cdot
2\pi\delta(\omega-\omega^{\prime}).$ Thus, we arrive at the intuitive
conclusion that the frequency domain correlations are proportional to the
frequency-dependent loss $\text{Re}\,K_{l}(\omega)$ that results from the
frequency-dependent out-coupling. Moreover, the correlator still retains a
delta function in the frequency domain, indicating that fluctuations at
different frequencies remain uncorrelated. This is the case because even
though the system we describe is dispersive, time-translation invariance is
still retained.

###  C Preservation of commutation relations

Now, we will prove that the model used to describe the non-Markovian damping
of the resonance $a$ preserves the commutation relations $[a,a^{\dagger}]=1$,
which is required for any quantum-mechanically consistent theory of
dissipation. We will continue to neglect nonlinearity: it can be shown that if
the equation without $\beta$ conserves the commutator, then the equation with
nonzero $\beta$ does as well. Since the operator equation of motion for $a$
(Eq. S7) is linear, it can be transformed into the frequency domain without
approximation. By doing this, we find that $a(\omega)$ is related to the
Langevin force as

$a(\omega)=\frac{F(\omega)}{-i(\omega_{a}-\omega)+K_{l}(\omega)}\equiv
i\xi(\omega)F(\omega).$ (S13)

As defined, $\xi(\omega)$ is a response function which governs how the input
field (and its fluctuations) couple into the resonance $a$ as a function of
frequency. We can then express the equal time commutation relation of $a$ as

$\displaystyle[a(t),a^{\dagger}(t)]$
$\displaystyle=\int\frac{d\omega}{2\pi}\frac{d\omega^{\prime}}{2\pi}e^{-i\omega
t}e^{i\omega^{\prime}t}[a(\omega),a^{\dagger}(\omega^{\prime})]$ (S14)
$\displaystyle=\int\frac{d\omega}{2\pi}\frac{d\omega^{\prime}}{2\pi}e^{-i\omega
t}e^{i\omega^{\prime}t}\xi(\omega)\xi^{*}(\omega^{\prime})[F(\omega),F^{\dagger}(\omega^{\prime})]$
(S15)
$\displaystyle=\int\frac{d\omega}{\pi}|\xi(\omega)|^{2}\mathrm{Re}\,K_{l}(\omega).$
(S16)

We note that in the final step, we have taken an expectation value of the
Langevin force commutator with respect to the reservoir. We see then that the
commutation relation of $a$ is expressed as a frequency-dependent integral
containing the frequency-dependent loss rate, as well as the response function
$\xi(\omega)$ which describes the linear coupling of the input field into the
cavity.

Now, we will show that this integral can be written in a simpler and more
fundamental form, using the analytic structure of $\xi$. In particular, we
note that

$|\xi(\omega)|^{2}=\frac{1}{(\omega_{a}+K_{l}^{\prime\prime}(\omega)-\omega)^{2}+K_{l}^{\prime}(\omega)^{2}},$
(S17)

where for brevity we will sometimes express the real and imaginary parts of
$K_{l}$ as $K_{l}=K_{l}^{\prime}+iK_{l}^{\prime\prime}$. We then note that the
imaginary part of $\xi$ is given as

$\text{Im}\,\xi(\omega)=\frac{K_{l}^{\prime}(\omega)}{(\omega_{a}+K_{l}^{\prime\prime}(\omega)-\omega)^{2}+K_{l}^{\prime}(\omega)^{2}}.$
(S18)

Thus, the commutator integrand can be written as

$[a(t),a^{\dagger}(t)]=\int\frac{d\omega}{\pi}|\xi(\omega)|^{2}\mathrm{Re}K_{l}(\omega)=\int\frac{d\omega}{\pi}\text{Im}\,\xi(\omega).$
(S19)

By using this form, the commutation relation is written as a frequency
integral over the imaginary part of a single function.

To conclude the proof, we will show that $\xi$ is a response function which
obeys its own Kramers-Kronig (K.K.) relation, and develop a sum rule which
evaluates the integral of Eq. S19. Recall that $\xi(\omega)$ is defined such
that in the absence of any input field, $a(\omega)=-i\xi(\omega)F(\omega)$.
This frequency-domain multiplicative relation corresponds to a time-domain
convolution $a(t)=-i\int dt^{\prime}\,\xi(t-t^{\prime})F(t^{\prime})$. Since
this system is linear and time invariant, the time domain response function
$\xi$ must be causal. Mathematically stated,
$\xi(\tau)=\theta(\tau)\xi(\tau)$, where $\theta(\tau)$ is the Heaviside step
function. This requirement of causality implies that in the frequency domain,
the response function $\xi(\omega)$ obeys the K.K. relation

$\text{Re}\,\xi(\omega)=-\mathcal{P}\int\frac{d\omega^{\prime}}{\pi}\frac{\text{Im}\,\xi(\omega^{\prime})}{\omega-\omega^{\prime}}.$
(S20)

The relevant sum rule can be developed by considering the asymptotic
$\omega\to\infty$ behavior of both sides of the K.K. relation. For the left
hand side, we have

$\text{Re}\,\xi(\omega)_{\omega\to\infty}\sim\frac{1}{\omega_{a}+K_{l}^{\prime\prime}(\omega)-\omega}\sim-\frac{1}{\omega}.$
(S21)

For the right hand side, we have

$\lim_{\omega\to\infty}-\mathcal{P}\int\frac{d\omega^{\prime}}{\pi}\frac{\text{Im}\,\xi(\omega^{\prime})}{\omega-\omega^{\prime}}\sim-\frac{1}{\omega}\int\frac{d\omega^{\prime}}{\pi}\text{Im}\,\xi(\omega^{\prime}).$
(S22)

By enforcing equality of the two expressions, we obtain the sum rule

$\int d\omega\,\text{Im}\,\xi(\omega)=\pi.$ (S23)

As a result of this sum rule, the commutation relation is now evaluated as

$[a(t),a^{\dagger}(t)]=\int\frac{d\omega}{\pi}\text{Im}\,\xi(\omega)=1.$ (S24)

Thus, we see that the commutation relation is preserved, so long as the
coupling and loss functions are K.K. consistent with one another as to
maintain causality.

##  S2 Steady-state mean-field behavior

In terms of this general formalism, we will now describe how the mean-field
and noise of $a$ behave in response to monochromatic coherent driving via the
input field $s(t)$. Physically, this corresponds to shining light into the
dispersive and nonlinear resonance via the coupling port (i.e., a waveguide),
and examining the corresponding response of the resonance field. We start by
decomposing the input field operator into the sum of a mean-field component,
and operator valued quantum fluctuations as $s(t)=\braket{s(t)}+\delta s(t)$.
Then, we assume that the input field is coherently and monochromatically
driven with amplitude $s_{0}$ and frequency $\omega_{p}$, so that
$\braket{s(t)}=s_{0}e^{-i\omega_{p}t}$. Additionally, we will denote $\langle
a\rangle=\alpha_{0}e^{-i\omega_{p}t}$ as the mean-field component of $a$,
which will also oscillate at the pump frequency. By substituting these
assumptions into Eq. S1, we find that in the steady-state, the mean-field
amplitude obeys:

$\begin{split}\frac{d}{dt}(\alpha_{0}e^{-i\omega_{p}t})=-i\omega_{a}\alpha_{0}e^{-i\omega_{p}t}&-i\beta|\alpha_{0}|^{2}\alpha_{0}e^{-i\omega_{p}t}\\\
&-\int dt^{\prime}K_{l}(t-t^{\prime})\alpha_{0}e^{-i\omega_{p}t^{\prime}}+\int
dt^{\prime}\,K_{c}(t-t^{\prime})s_{0}e^{-i\omega_{p}t^{\prime}}.\end{split}$
(S25)

After performing Fourier manipulations, and cancelling out oscillating terms,
we find the frequency domain steady-state equation

$[i(\omega_{a}-\omega_{p})+i\beta|\alpha_{0}|^{2}+K_{l}(\omega_{p})]\alpha_{0}=s_{0}K_{c}(\omega_{p}).$
(S26)

Thus, we see that the monochromatic drive causes the loss and coupling kernels
to be sampled at the pump frequency $\omega_{p}$. To obtain an equation which
describes the steady-state photon number $n=|\alpha_{0}|^{2}$ rather than just
the field, we note that the loss function $K_{l}$ is complex. We will
specifically write
$K_{l}(\omega)=K_{l}^{\prime}(\omega)+iK_{l}^{\prime\prime}(\omega).$ This
decomposition comes with a clear interpretation, namely that
$K_{l}^{\prime}(\omega)$ encodes the frequency-dependent damping rate of the
resonance $a$, while $K_{l}^{\prime\prime}(\omega)$ encodes the shift of the
resonance frequency $\omega_{a}$ which accompanies the loss. In other words,
the real and imaginary parts of $K_{l}$ respectively correspond to the
dissipative and reactive dynamics induced by the frequency-dependent coupling
to the environment. Using this decomposition, an equation for the steady-state
resonance photon number $n$ can be written as:

$[(\omega_{a}+\beta
n+K_{l}^{\prime\prime}(\omega_{p})-\omega_{p})^{2}+K_{l}^{\prime}(\omega_{p})^{2}]n=|s_{0}|^{2}|K_{c}(\omega_{p})|^{2}.$
(S27)

We note that via Kramers-Kronig relations between the loss and coupling
kernels, $|K_{c}(\omega)|^{2}=2\,\text{Re}\,K_{l}(\omega)$. This equation of
state is consistent with results known in the context of Kerr induced optical
bistability drummond1980quantum . The main difference here is that in-
coupling, damping, and detuning are now frequency-dependent. In general, this
means that the presence and extent of bistable behavior can now depend heavily
on the frequency at which the system is driven.

##  S3 Steady-state quantum noise properties

In addition to the steady-state behavior, the quantum noise properties of this
driven system can also differ substantially from the behaviors seen in
ordinary Kerr bistability. In this section, we will develop the theory of
quantum noise in the non-Markovian driven-dissipative cavity, and develop
analytical approximations which lead to the “sharp loss” interpretation of
intensity noise reduction.

###  A General formalism for quantum noise

First, we develop the general theory of quantum noise in a non-Markovian Kerr
resonator. We do so through a linearization approach, which decomposes the
quantum cavity and input fields into c-number mean values, plus operator
valued fluctuations. For the input field, this is $s=\braket{s}+\delta s$,
while for the cavity field it is $a=\braket{a}+\delta a$. By substituting this
decomposition into Eq. S1, and maintaining terms of only linear order in the
fluctuations, we can obtain linear equations of motion for the fluctuations.
In the frequency domain, the linearized operator-valued fluctuations obey the
equations:

$\begin{pmatrix}\eta(\omega)&i\beta n\\\ -i\beta
n&\eta^{*}(-\omega)\end{pmatrix}\begin{pmatrix}\delta a(\omega)\\\ \delta
a^{\dagger}(-\omega)\end{pmatrix}=\begin{pmatrix}K_{c}(\omega_{p}+\omega)\delta
s(\omega)\\\ K_{c}^{*}(\omega_{p}-\omega)\delta
s^{\dagger}(-\omega)\end{pmatrix}$ (S28)

Here, we have defined $\eta(\omega)\equiv i(\omega_{ap}-\omega+2\beta
n)+K_{l}(\omega_{p}+\omega)$, where $\omega_{ap}\equiv\omega_{a}-\omega_{p}$
is the bare pump-cavity detuning. Additionally, we note that any instance of
$n$ refers to the steady state mean-field dictated by Eq. S27. At a high
level, this result indicates that the combination of frequency-dependent and
Kerr interactions results in a frequency-dependent mixing of the two
quadratures of the resonance $a$. This will result in frequency-dependent
squeezing phenomena, that can result in the reduction of fluctuations below
the shot noise limit. Additionally, we note that if $\delta s$ represents
classical fluctuations from the mean-field steady-state solution, then the
equation above can also be used to describe classical oscillations from
equilibrium. A discussion of these transient behaviors is given in Section S5.

The matrix equation above can be inverted to solve for $\delta a(\omega)$ and
$\delta a^{\dagger}(-\omega)$ in terms of the input fluctuations $\delta
s(\omega)$ and $\delta s^{\dagger}(-\omega)$. We define quantities $p(\omega)$
and $q(\omega)$ such that $\delta a(\omega)=p(\omega)\delta
s(\omega)+q^{*}(-\omega)\delta s^{\dagger}(-\omega)$. In terms of previously
defined quantities,

$\displaystyle p(\omega)$
$\displaystyle=\frac{1}{M(\omega)}\eta^{*}(-\omega)K_{c}(\omega_{p}+\omega)$
(S29) $\displaystyle q(\omega)$ $\displaystyle=\frac{i\beta
n}{M(\omega)}K_{c}(\omega_{p}+\omega)$ (S30)

where we have also defined
$M(\omega)\equiv\eta(\omega)\eta^{*}(-\omega)+(\beta n)^{2}$ as the
determinant of the matrix acting on the fluctuations.

In addition to the field operators $a$ and $a^{\dagger}$, we can also consider
quadrature operators $X=a+a^{\dagger}$ and $Y=-i(a-a^{\dagger})$. This allows
us to analytically express the steady-state variance of the quadrature
operators in terms of the functions $p(\omega)$ and $q(\omega)$ defined above.
For example, the variance of the $X$ quadrature is given by

$\displaystyle\braket{\delta X(t)\delta X(t)}$
$\displaystyle=\int\frac{d\omega}{2\pi}\frac{d\omega^{\prime}}{2\pi}\,e^{-i(\omega+\omega^{\prime})t}\braket{\delta
X(\omega)\delta X(\omega^{\prime})}$ (S31)
$\displaystyle=\int\frac{d\omega}{2\pi}\frac{d\omega^{\prime}}{2\pi}\,e^{-i(\omega+\omega^{\prime})t}\Braket{\left(\delta
a(\omega)+\delta a^{\dagger}(-\omega)\right)\left(\delta
a(\omega^{\prime})+\delta a^{\dagger}(-\omega^{\prime})\right)}$ (S32)
$\displaystyle=\int\frac{d\omega}{2\pi}\frac{d\omega^{\prime}}{2\pi}\,e^{-i(\omega+\omega^{\prime})t}\left(p(\omega)+q(\omega)\right)\left(p^{*}(-\omega^{\prime})+q^{*}(-\omega^{\prime})\right)\Braket{\delta
s(\omega)\delta s^{\dagger}(-\omega^{\prime})}$ (S33) $\displaystyle(\Delta
X)^{2}$ $\displaystyle=\int\frac{d\omega}{2\pi}|p(\omega)+q(\omega)|^{2}.$
(S34)

In going from the second to the third line, we have expanded $\delta
a(\omega)$ and its conjugate in terms of $p(\omega)$ and $q(\omega)$, and
dropped all terms which have $s$-correlators which are zero for a vacuum
reservoir. Using the fact that $\Braket{\delta s(\omega)\delta
s^{\dagger}(\omega^{\prime})}=2\pi\delta(\omega-\omega^{\prime})$ for a vacuum
reservoir allows us to integrate over $\omega^{\prime}$ to achieve the final
result. More general inputs (such as squeezed-state inputs) can be considered
by considering more general input correlations.

After doing the same for $Y$, we find that the quadrature variances can be
written as

$\displaystyle(\Delta X)^{2}$
$\displaystyle=\int\frac{d\omega}{2\pi}|p(\omega)+q(\omega)|^{2}$ (S35)
$\displaystyle(\Delta Y)^{2}$
$\displaystyle=\int\frac{d\omega}{2\pi}|p(\omega)-q(\omega)|^{2}.$ (S36)

At this stage, no approximations have been made. As such, the evaluation of
these frequency integrals result in exact results for the quadrature
variances. Without loss of generaality, we choose the pump phase so that
$\braket{a}$ in the steady-state is real. Then, $(\Delta X)^{2}$ denotes the
amplitude quadrature variance, while $(\Delta Y)^{2}$ denotes the phase
quadrature variance. The commutation relation $[a,a^{\dagger}]$ for the
quantized resonance imposes the uncertainty relation $(\Delta X)(\Delta Y)\geq
1$ on the quadrature variances, which is a manifestation of the amplitude-
phase uncertainty relationship on the cavity field. Thus, the reductions in
intensity noise we consider in this work are necessarily accompanied by an
increase in fluctuations in the phase quadrature. The ultimate limit of this
type of nonlinear dissipation (which is formally no longer fully captured
through this linearization approach) is a Fock state with perfectly defined
amplitude, and a completely undefined phase.

We can rewrite these noise expressions in forms which make them more
interpretable, and also more amenable to analytical approximations. To start,
we will examine the quantity $M(\omega)$, which originates from the
determinant of the linearized fluctuation matrix. By using the definition of
$\eta$, we can write
$M(\omega)=\Omega^{2}(\omega)-\omega^{2}+i\Gamma^{2}(\omega)$, where

$\displaystyle\Omega^{2}(\omega)$ $\displaystyle\equiv 3(\beta
n)^{2}+(\omega_{ap}+\delta_{-})(\omega_{ap}+\delta_{+})+2\beta
n(2\omega_{ap}+\delta_{+}+\delta_{-})+(\delta_{+}-\delta_{-})\omega+\kappa_{+}\kappa_{-}$
(S37) $\displaystyle\Gamma^{2}(\omega)$ $\displaystyle\equiv 2\beta
n(\kappa_{-}-\kappa_{+})+\kappa_{-}(\omega_{ap}+\delta_{+}-\omega)-\kappa_{+}(\omega_{ap}+\delta_{-}+\omega).$
(S38)

We have expressed the loss function in terms of its real and imaginary parts
as $K_{l}(\omega)=K_{l}^{\prime}(\omega)+iK_{l}^{\prime\prime}(\omega)$.
Additionally, we have adopted the shorthand notation $\kappa_{\pm}\equiv
K_{l}^{\prime}(\omega_{p}\pm\omega)$ (and likewise for $\delta$) to denote the
loss (and detuning) evaluated around the pump frequency $\omega_{p}$. As a
reminder, it is understood that $n$ refers to the steady-state value of $n$.

This means then that the denominator of all of these noise and commutator
quantities will take the form

$\frac{1}{|M(\omega)|^{2}}=\frac{1}{\left[\Omega^{2}(\omega)-\omega^{2}\right]^{2}+\Gamma^{4}(\omega)}.$
(S39)

In this case, the $X$-quadrature variance can be expressed as

$(\Delta
X)^{2}=\int\frac{d\omega}{\pi}\frac{K_{+}^{\prime}(\omega)\left[(\omega_{ap}+\beta
n+\omega)^{2}+K_{-}^{\prime
2}(\omega)\right]}{\left[\Omega^{2}(\omega)-\omega^{2}\right]^{2}+\Gamma^{4}(\omega)}.$
(S40)

Similarly, the $Y$-quadrature variance can be expressed as

$(\Delta
Y)^{2}=\int\frac{d\omega}{\pi}\frac{K_{+}^{\prime}(\omega)\left[(\omega_{ap}+3\beta
n+\omega)^{2}+K_{-}^{\prime
2}(\omega)\right]}{\left[\Omega^{2}(\omega)-\omega^{2}\right]^{2}+\Gamma^{4}(\omega)}$
(S41)

###  B Analytical approximations for quantum noise in the adiabatic limit

Now, we focus on developing analytic expressions for the intensity quadrature
variance $(\Delta X)^{2}$. Our approximations are based on the so-called
_adiabatic approximation_ , which holds when the non-Markovian elements
responds quickly in time compared to the resonance $a$. In this case, we find
that the noise spectra are sharply peaked at frequencies
$\omega_{p}\pm\Omega$. Here, $\Omega$ is the so-called relaxation oscillation
frequency of a Markovian Kerr cavity, which is defined as
$\Omega=\sqrt{\Delta^{2}-(\beta n)^{2}}$, where $\Delta=\omega_{ap}+2\beta n$.
The other portions of the noise integral, on the other hand, are slowly
varying, and thus the integral can be approximated as a pair of Lorentzian
integrals peaked at $\pm\Omega$ to a high degree of accuracy. After these
integrations are performed, we find an analytical expression for the $X$
quadrature variance of the cavity field:

$\displaystyle(\Delta X)^{2}$
$\displaystyle\approx\frac{1}{2\Omega\Gamma^{2}}\left[\kappa_{+}((\Delta-\beta
n+\Omega)^{2}+\kappa_{+}^{2})+\kappa_{-}((\Delta-\beta
n-\Omega)^{2}+\kappa_{-}^{2})\right]$ (S42)
$\displaystyle\approx\frac{1}{2\Omega\Gamma^{2}}\left[\kappa_{+}(\Delta-\beta
n+\Omega)^{2}+\kappa_{-}(\Delta-\beta n-\Omega)^{2}\right].$ (S43)

In the above expressions, $\kappa_{\pm}=K_{l}^{\prime}(\omega_{p}\pm\Omega)$
is the dispersive loss evaluated at the pump shifted up and down by the
frequency $\Omega$. Additionally, the linewidth $\Gamma$ is defined as

$\Gamma^{2}=|\Delta(\kappa_{-}-\kappa_{+})-\Omega(\kappa_{-}+\kappa_{+})|$
(S44)

In the second line of the expression for the noise, we make the approximation
that the terms cubic in the loss are small compared to the detunings. This is
well respected across the regimes of parameters that we consider in this work.
The above expressions constitute one of our key theoretical results. As we
will see, this expression for the quadrature variance depends heavily on the
dispersion of the loss. In particular, we see that the dispersive loss
$K_{l}^{\prime}(\omega)$ is evaluated at two key points:
$\kappa_{\pm}=K_{l}^{\prime}(\omega_{p}\pm\Omega)$. As a reminder,
$\omega_{p}$ is the pump frequency, and $\Omega$ is a frequency related to
resonance of fluctuations, related to how the nonlinear detuning changes with
$n$. Additionally, we see that the expression for $\Gamma$ (corresponding to
the damping rate of the integrated peaks), depends on both the difference and
the sum of the loss at the two points $\omega_{p}\pm\Omega$. This is the
mechanism by which the noise reduction depends on the dispersive loss. Roughly
speaking, the sum and difference terms present in $\Gamma$ serve as probes for
the mean and slope of the dispersive loss. Down below, we will show that in
cases of maximal noise reduction near the edge of a bistable branch, these
interpretations become exact, yielding compact and intuitive expressions for
the intensity quadrature noise.

The expression for the noise provided by Eq. S43 can be cast into a more
interpretable form. After algebraic manipulations, we find that it can be
written as

$(\Delta X)^{2}\approx\left(1-\frac{\beta
n}{\Delta}\right)\frac{\left(\frac{\Delta}{\Omega}\right)^{2}+r\left(\frac{\Delta}{\Omega}\right)}{1+r\left(\frac{\Delta}{\Omega}\right)}.$
(S45)

Here, $r\equiv(\kappa_{+}-\kappa_{-})/(\kappa_{+}+\kappa_{-})$ is the ratio
between the difference and sum in the loss probed at $\omega_{p}\pm\Omega$.
Eq. S45 constitutes one of our key theoretical results. It describes how the
amplitude noise of an externally pumped non-Markovian Kerr resonance is
determined entirely by the steady state cavity photon number $n$, detuning
parameters, and the dispersive loss evaluated at sideband frequencies
$\pm\Omega$.

Important limiting cases

Here, we discuss some important limiting cases for out analytical expression
for the amplitude quadrature variance derived above.

Same loss at both frequencies

In this case, we have $\kappa_{+}=\kappa_{-}\equiv\kappa$. When this is the
case, we get $\Gamma^{2}=2\kappa\Omega$. We end up with

$(\Delta X)^{2}=\frac{1}{2}\left[1+\frac{\Delta-\beta n}{\Delta+\beta
n}\right].$ (S46)

The minimum possible value is $1/2$, and occurs when $\Delta=\beta n$, which
occurs when $\beta n=\omega_{p}-\omega_{a}$. When these is no nonlinearity
present, we have $\beta=0$ which immediately gives a quadrature variance of 1.
This is what is to be expected from a linear cavity which is coherently
pumped, as such a system produced a steady-state field which is in a coherent
state whose two quadrature variances are equal and uncertainty limited.

Large $n$ limit

When $n$ becomes large, we additionally have $\Delta\to 2\beta n$, which gives

$(\Delta X)^{2}\to\frac{2}{3}.$ (S47)

This limit is consistent with, for example, what has been noted in
drummond1980quantum  which investigated the quantum statistical properties of
optical bistability in (Markovian) Kerr cavities. In other words, some limited
degree of intensity squeezing is possible in a nondisperisve Kerr cavity. On
physical grounds, this effect results from different parts of the phase space
distribution experiencing different rotation rates due to the intensity-
dependent frequency characteristic of a cavity that contains a Kerr medium.
Thus the phase space distribution of light in a pumped Kerr medium is deformed
away from its symmetric coherent state form, leading to some squeezing. Note
that this limit also effectively falls under the category of the losses at the
two points becoming the same, because $\Omega$ grows linearly with $n$ for
large $n$, and thus the dispersive loss is evaluated far away from the
dispersive part. Within the context of our loss model, this means that both
evaluation points are at opposite ends of the Lorentzian dip, where the
dispersion is flat, and the loss values are the same.

Sharp loss

The expression for amplitude noise given in Eq. S45 allows us to make an
important interpretation related to “sharp loss” when the system is pumped
near a bistable point. In particular, the onset of a bistable point is
associated with the vanishing of $\Omega$. As one approaches the bistable
point, $\Omega$ thus becomes small until its eventual vanishing. When $\Omega$
becomes small relative to the scale of frequency variations in the dispersion
loss $K_{l}^{\prime}(\omega)$, this expression for the amplitude noise probes
the slope of $K_{l}^{\prime}(\omega)$. To see this we note that

$\lim_{\Omega\to 0}\frac{r}{\Omega}=\lim_{\Omega\to
0}\frac{K_{l}^{\prime}(\omega+\Omega)-K_{l}^{\prime}(\omega-\Omega)}{\Omega(K_{l}^{\prime}(\omega+\Omega)+K_{l}^{\prime}(\omega-\Omega))}=\frac{1}{K_{l}^{\prime}(\omega_{p})}\frac{dK_{l}^{\prime}(\omega)}{d\omega}\Bigg{|}_{\omega=\omega_{p}}.$
(S48)

Thus, as one approaches the bistable points, the noise is increasingly
determined by the slope of the loss curve at the pump frequency, relative to
the loss at this point. It is clear, then, that in order to maximize
squeezing, one should: (1) operate with a detuning which is sufficiently large
to induce bistability, (2) operate at a photon number $n$ at the lower region
of the upper bistable branch, which creates a small sideband frequency
$\Omega$, and (3) ensure that the ratio
$\frac{1}{K_{l}^{\prime}(\omega_{p})}\frac{dK_{l}^{\prime}(\omega)}{d\omega}\Big{|}_{\omega=\omega_{p}}$
is as large as possible, by maximizing the steepness of the dispersive loss,
while minimizing all other background losses. Examples of these principles are
explored further in the next section.

###  C Additional results for for generation of highly intensity-squeezed
states

Figure S1: Large intensity squeezing of macroscopic states. (a) Dispersive
loss curve associated with the F.W. model. Curve is marked with four pump
frequencies which are represented in (b). (b) Fano factor as a function of
steady state photon number for different pump frequencies $\omega_{p}$. A
comparison between the exact numerical integral (Eq. S41) and the analytical
approximation (Eq. S45) shows very strong agreement. Parameters used are the
same as those used in main text Figs. 2-3. Figure S2: Near-fock state
generation with the F.W. model. (a) Disperisve loss dip caused by
interference. Values of the bare resonance $\omega_{a}$ and pump $\omega_{p}$
are shown, which are chosen to lie near the lossless mode (BIC). (b) Input-
output curve which exhibits bistability for small photon numbers. (c)
Amplitude noise reduction characterized by the Fano factor. Near-Fock state
generation occurs around $n=42$ at the lower part of the upper bistable
branch. Comparisons to other noise limits are shown. Parameters used are:
$\beta=10^{-4}\omega_{a}$, $\kappa=10^{-4}\omega_{a}$,
$\gamma=10^{-2}\omega_{a}$, $\omega_{d}=\omega_{a}+0.4\gamma$,
$\omega_{p}=1.0042\omega_{a}$.

In this section, we provide additional data for high intensity squeezing that
can be achieved through non-Markovian dissipation. We first provide additional
data for intensity noise reduction for the nonlinear F.W. system operated with
macroscopic photon numbers (the focus of Fig. 3 of the main text). In
particular, we show the amplitude noise results of main Fig. 3a, but for four
different pump frequencies (Fig. S1). By operating at a pump frequency just to
the right of the vanishing loss point created by the BIC, it is possible to
achieve a very high degree of steady state squeezing in the cavity ($>15$ dB
in this case).

We remark on a few additional features of the amplitude noise data. First, we
note that the amplitude noise exhibits divergences at lower photon numbers,
due to a combination of modulational instability and ordinary bistability.
Noise values are not shown in these regions, since there is no steady state
for which to compute the noise. Noise reduction below the Markovian limit of
$1/2$ is seen for the three pump frequencies which lie to the right of the BIC
point. In these case, it is the “sharp loss” mechanism we have described that
stabilizes these intensity-squeezed steady-states. Notably, when the pump
frequency is on the left side of the BIC (where loss decreases with
frequency), no such noise reduction is achieved. This pump frequency is
actually below the detuning required for bistability, so there is only a
single region of diverging noise, which results from modulational instability,
even in the absence of bistability. Finally, all of the models approach a Fano
factor of $2/3$ as $n$ increases, as described in the previous section.

We also provide additional details on the generation of near-Fock cavity
states through the coherent pumping of nonlinear non-Markovian systems. We
show how this can be realized through both the nonlinear F.W. model (Fig. S2),
as well as the Fano mirror model (Fig. S3). For both models, we assume a
nonlinear strength $\beta=10^{-4}\omega_{a}$, which is characteristic of the
extreme nonlinearities found in exciton-polariton condensates
fink2018signatures .

In both models, near-Fock states of large number (order of tens) can be
created by operating around the frequency range where the losses drop due to
interference. In both systems, near-Fock state generation is optimized by
setting the bare cavity frequency $\omega_{a}$ slightly below the loss dip,
while pumping at frequency $\omega_{p}$ slightly to the right of the loss dip.
This ensures that bistable operation can occur, and that the pump frequency
lies just to the right of the loss minimum, where the loss increases as a
function of frequency. This enables one to exploit the “sharp loss” mechanism
described in order to achieve the highest possible noise reduction. The lowest
noise states can be created inside the cavity by operating at the very lowest
part of the upper bistable branch (Figs. S2b, S3b), which can be accessed
through hysteresis on the bistable input-output curve.

For the parameters chosen here, the F.W. model produces a state with $n\approx
42$ and $\Delta n\approx 1$, while the Fano mirror model produces one with
$n\approx 70$ and $\Delta n\approx 1$. The non-Markovian architecture employed
here gives these states a degree of amplitude squeezing which is far beyond
the bounds set by other standard architectures. For starters, these states
quite easily exceed the Fano factor bound $F>1/2$ which is attainable with a
conventional (Markovian) driven-dissipative cavity. Additionally, these states
easily exceed the bound for the amplitude squeezing created by the
displacement of squeezed vacuum ($(\Delta n)^{2}\geq n^{2/3}$)
bondurant1984squeezed , and even the bound for what is possible with more
elaborate Kerr squeezing architectures ($(\Delta n)^{2}\geq n^{1/3}$)
kitagawa1986number . The non-Markovian driven-dissipative architecture thus
presents attractive opportunities for creating high number Fock or near-Fock
states which are impossible to create by more conventional means.

The examples given here, as well as the more general physical insights
provided by our analytical approximations, indicate that the generation of
highly intensity-squeezed (and even near-Fock states) in nonlinear non-
Markovian driven-dissipative systems is a phenomenon that can be realized
quite broadly.

Figure S3: Near-fock state generation with Fano mirror. Similar setup to that
shown in Fig. S2, except that the near-Fock state generation is done with a
Fano mirror model. Parameters used are: $\omega_{a}=1.03\times 10^{15}$ s-1,
$\kappa=10^{-4}\omega_{a}$, $L=5$ $\mu$m, and $\beta=10^{-4}\omega_{a}$, ,
which are chosen to represent systems with strong per-photon nonlinearity
fink2018signatures .

##  S4 non-Markovian loss models

In these sections, we briefly summarize the two models of non-Markovian
dissipation which are used in this work. Detailed derivations are provided in
rivera2023creating .

###  A Friedrich-Wintgen model

The Friedrich-Wintgen (F.W.) model is a temporal coupled mode theory model
which exhibits bound state in the continuum (BIC) formation via two resonances
$a$ and $d$ coupled to a common continuum friedrich1985interfering ;
hsu2016bound . We consider a quantum treatment of the F.W. model, in which one
of the resonances ($a$) is additionally nonlinear. This nonlinear F.W. model
thus comprises two quantized resonances, $a$ and $d$, which are coupled to a
common continuum $s(t)$. The resonance $a$ is assumed to be nonlinear. In this
case, the Heisenberg-Langevin equations of motion for the mode operators are:

$\displaystyle\dot{a}$
$\displaystyle=-i\omega_{a}a-i\beta(a^{\dagger}a)a-\kappa
a-\sqrt{\kappa\gamma}d+\sqrt{2\kappa}s(t)$ (S49a) $\displaystyle\dot{d}$
$\displaystyle=-i\omega_{d}d-\gamma d-\sqrt{\kappa\gamma}a+\sqrt{2\gamma}s(t)$
(S49b)

Here, $\kappa$ and $\gamma$ are the respective decay rates of $a$ and $d$.
Since $d$ obeys linear dynamics, it can be eliminated to yield an equation of
motion for $a$ alone, which takes the form of Eq. S1.

For the F.W. model, the coupling and loss functions for the resonance $a$ are:

$\displaystyle K_{c}(\omega)$
$\displaystyle=\sqrt{2\kappa}\left[1-\frac{\gamma}{i(\omega_{d}-\omega)+\gamma}\right]$
(S50a) $\displaystyle K_{l}(\omega)$
$\displaystyle=\kappa\left[1-\frac{\gamma}{i(\omega_{d}-\omega)+\gamma}\right],$
(S50b)

where $\kappa$ and $\gamma$ are the respective decay rates of $a$ and $d$,
which have frequencies $\omega_{a,d}$.

###  B Fano mirror model

For the Fano mirror model, we use the following coupling and loss functions to
describe the dispersive dissipation of $a$:

$\displaystyle K_{l}(\omega)$
$\displaystyle=2\kappa\left(1-\frac{e^{2i\theta_{2}}}{r_{d}-e^{-i\omega
T}}\right)=2\kappa\left(1-\frac{r_{d}-it_{d}\sigma}{r_{d}-e^{-i\omega
T}}\right)$ (S51) $\displaystyle K_{c}(\omega)$
$\displaystyle=\sqrt{2\kappa}e^{i\theta_{1}}\left[1+\frac{it_{d}e^{i(\theta_{2}-\theta_{1})}}{r_{d}-e^{-i\omega
T}}\right].$ (S52)

Here, $\kappa$ is the decay rate of $a$, $\theta_{1,2}$ are phases associated
with the interference, $\sigma=\pm$ is the parity, $r_{d}$ and $t_{d}$ are the
reflection and transmission coefficients of the interference, and $T$ is the
round trip time of the cavity.

##  S5 Modulational instabilities and phase diagrams

In this section, we discuss the stability properties of non-Markovian driven-
dissipative cavities. We begin with a stability analysis of the nonlinear F.W.
model, and then describe how these conclusions can be more broadly
generalized.

In steady state operation, we can describe the quantum state of the field in
terms of fluctuations about the steady state values. The equations of motion
for the fluctuations associated with Eqs. S49a, S49b are:

$\displaystyle-i\omega_{p}\delta a+\delta\dot{a}$
$\displaystyle=(-i\omega_{a}-\kappa)\delta a-i\beta n(2\delta a+\delta
a^{\dagger})-\sqrt{\kappa\gamma}\delta d+\sqrt{2\kappa}\delta s$ (S53)
$\displaystyle i\omega_{p}\delta a^{\dagger}+\delta\dot{a}^{\dagger}$
$\displaystyle=(i\omega_{a}-\kappa)\delta a^{\dagger}+i\beta n(2\delta
a^{\dagger}+\delta a)-\sqrt{\kappa\gamma}\delta
d^{\dagger}+\sqrt{2\kappa}\delta s^{\dagger}$ (S54)
$\displaystyle-i\omega_{p}\delta d+\delta\dot{d}$
$\displaystyle=(-i\omega_{d}-\gamma)\delta d-\sqrt{\kappa\gamma}\delta
a+\sqrt{2\gamma}\delta s$ (S55) $\displaystyle i\omega_{p}\delta
d^{\dagger}+\delta\dot{d}^{\dagger}$ $\displaystyle=(i\omega_{d}-\gamma)\delta
d^{\dagger}-\sqrt{\kappa\gamma}\delta a^{\dagger}+\sqrt{2\gamma}\delta
s^{\dagger}.$ (S56)

The two resonances can also be described in terms of quadrature operators
$X_{a}=a+a^{\dagger}$, $Y_{a}=-i(a-a^{\dagger})$ (and likewise for $d$). Then
fluctuations in the quadratures from their mean values (for example,
$X_{a}=\braket{X_{a}}+\delta X_{a}$). The linearized equation of motion for
these quadratures in the case of the nonlinear F.W. model are given as:

$\frac{d}{dt}\begin{pmatrix}\delta X_{a}\\\ \delta Y_{a}\\\ \delta X_{d}\\\
\delta Y_{d}\end{pmatrix}=\begin{pmatrix}-\kappa&\omega_{ap}+\beta
n&-\sqrt{\kappa\gamma}&0\\\ -\omega_{ap}-3\beta
n&-\kappa&0&-\sqrt{\kappa\gamma}\\\
-\sqrt{\kappa\gamma}&0&-\gamma&\omega_{dp}\\\
0&-\sqrt{\kappa\gamma}&-\omega_{dp}&-\gamma\end{pmatrix}\begin{pmatrix}\delta
X_{a}\\\ \delta Y_{a}\\\ \delta X_{d}\\\ \delta
Y_{d}\end{pmatrix}+\begin{pmatrix}F_{X}^{a}\\\ F_{Y}^{a}\\\ F_{X}^{d}\\\
F_{Y}^{d}\end{pmatrix}$ (S57)

The source terms on the right hand side are operator-valued Langevin forces,
which are defined as:

$\displaystyle F_{X}^{a}$ $\displaystyle=\sqrt{2\kappa}\left(\delta s+\delta
s^{\dagger}\right)$ (S58) $\displaystyle F_{Y}^{a}$
$\displaystyle=-i\sqrt{2\kappa}\left(\delta s-\delta s^{\dagger}\right)$ (S59)
$\displaystyle F_{X}^{d}$ $\displaystyle=\sqrt{2\gamma}\left(\delta s+\delta
s^{\dagger}\right)$ (S60) $\displaystyle F_{Y}^{d}$
$\displaystyle=-i\sqrt{2\gamma}\left(\delta s-\delta s^{\dagger}\right).$
(S61)

By finding the eigenvalues of the matrix in the above equations of motion, we
can determine the nature of the stability. Although Eq. S57 has been developed
in the full quantum picture, it can also be reduced to a classical limit by
taking expectation values of both sides. In this case, the Langevin forces
(which have zero expectation value) will vanish, leaving only a a classical
linear matrix equation with no sources. Such an equation governs the
linearized dynamics of the quadrature variables away from the mean-field
steady-state behavior.

Figure S4: Stability eigenvalues in the nonlinear F.W. model. (a) Complex
eigenvalues of the nonlinear F.W. model as a function of steady state photon
number $n$. (b) Zoom in of the small region of (a) over which one of the
eigenvalues acquires a small positive real part, leading to modulational
instability. (c) Phase diagram of the nonlinear F.W. model which shows the
amount of MI gain as a function of pump frequency $\omega_{p}$ and photon
number $n$. Gain is seen to be highest for regions of phase space in proximity
to the conventional unstable region.

In Fig. S4, we provide more details about the eigenvalues of the nonlinear
F.W. model, obtained by diagonalizing the matrix in Eq. S57. Each of the
eigenvalues may have an imaginary part $\Omega$, and a real part $\Gamma$,
which dictate the decaying (or growing) oscillations which occur when the
system is perturbed from equilibrium. In particular, $\Gamma>0$ is taken to
refer to damped oscillations, while $\Gamma<0$ is taken to refer to growing
oscillations. Fig. S4a shows the real and imaginary parts of the eigenvalues
as a function of steady state photon number $n$. The region in the middle
marked as “unstable” has an eigenvalue with positive real part, while
$\Omega=0$. This region is unstable due to conventional cavity bistability,
and trajectories in this region will be pushed away from the region. It can
also be seen that one of the eigenvalue pairs have imaginary parts $\pm\Omega$
which approach zero near these boundaries.

The region of photon numbers which give rise to modulational instability is
highlighted in Fig. S4b. This region hosts a small band of photon numbers over
which $\Gamma<0$, corresponding to growth of fluctuations. Fig. S4c shows the
MI gain $-\Gamma$ in the context of the whole phase space of $n$ and
$\omega_{p}$. We see that damping occurs in the majority of the phase diagram
(blue), while modulational instability gain occurs over two bands which emerge
from the conventional unstable region. Moreover, the MI gain is highest in the
vicinity of the conventional unstable region, and tapers off toward zero with
larger detunings from $\omega_{p}$ in either direction.

Figure S5: F.W. model phase diagrams for different adiabatic ratios. Phase
diagrams of the driven nonlinear F.W. model are shown. Parameters used are:
$\beta=10^{-10}\omega_{a}$, $\gamma=10^{-2}\omega_{a}$,
$\omega_{d}=\omega_{a}+\gamma$. Values of $\kappa$ are varied to produce the
adiabatic ratios $\kappa/\gamma=\\{0.1,1,10\\}$ marked in the panels.

In Fig. S5, we show how the phase diagram of the nonlinear F.W. model depends
on the leakage rates $\kappa$ and $\gamma$. For the three panels shown,
$\gamma$ is fixed, while $\kappa$ is varied. All three sets of parameters have
regions of instability due to the traditional bistability, as well as the
self-pulsing induced by modulational instability (MI).

More general models

In the adiabatic limit, we find that the the real part of the larger
eigenvalue is well-approximated by

$\text{Re}\,\lambda=-\frac{\kappa_{+}+\kappa_{-}}{2}-\frac{\Delta}{\Omega}\frac{\kappa_{+}-\kappa_{-}}{2}.$
(S62)

We thus see that this can vanish, and thus change sign, when

$1+\frac{\Delta}{\Omega}\frac{\kappa_{+}-\kappa_{-}}{\kappa_{-}+\kappa_{+}}=0.$
(S63)

This is exactly the condition see in the denominator of the noise expressions
in the adiabatic limit (see Eq. S45). Thus, this expression can be used to
determine the boundaries of modulational instability, even in more general
cases, so long as the adiabatic approximation is respected. This is how we
have determined the phase diagram in the Fano mirror model (Fig. 4 of the main
text). We find the MI boundaries predicted by this method to be in strong
agreement with those which are exactly determined by observing where the
amplitude and phase noise diverge using the exact expression for the noise
(Eq. S41), as evaluated by numerical integration.

## References

  * [1] Nicholas Rivera, Jamison Sloan, Yannick Salamin, John D Joannopoulos, and Marin Soljačić.  Creating large fock states and massively squeezed states in optics using systems with nonlinear bound states in the continuum.  Proceedings of the National Academy of Sciences, 120(9):e2219208120, 2023.
  * [2] Marlan O Scully and M Suhail Zubairy.  Quantum optics, 1999.
  * [3] PD Drummond and DF Walls.  Quantum theory of optical bistability. i. nonlinear polarisability model.  Journal of Physics A: Mathematical and General, 13(2):725, 1980\.
  * [4] Thomas Fink, Anne Schade, Sven Höfling, Christian Schneider, and Ataç Imamoglu.  Signatures of a dissipative phase transition in photon correlation measurements.  Nature Physics, 14(4):365–369, 2018.
  * [5] Roy S Bondurant and Jeffrey H Shapiro.  Squeezed states in phase-sensing interferometers.  Physical Review D, 30(12):2548, 1984.
  * [6] M Kitagawa and Y Yamamoto.  Number-phase minimum-uncertainty state with reduced number uncertainty in a kerr nonlinear interferometer.  Physical Review A, 34(5):3974, 1986.
  * [7] H Friedrich and D Wintgen.  Interfering resonances and bound states in the continuum.  Physical Review A, 32(6):3231, 1985.
  * [8] Chia Wei Hsu, Bo Zhen, A Douglas Stone, John D Joannopoulos, and Marin Soljačić.  Bound states in the continuum.  Nature Reviews Materials, 1(9):1–13, 2016.