# Single-photon nonlinearities and blockade from a strongly driven photonic
molecule

Davide Nigro  Dipartimento di Fisica, Università di Pavia, via Bassi 6,
I-27100 Pavia, Italy  Marco Clementi  Dipartimento di Fisica, Università di
Pavia, via Bassi 6, I-27100 Pavia, Italy  Photonic Systems Laboratory (PHOSL),
Swiss Federal Institute of Technology Lausanne (EPFL), CH-1015 Lausanne,
Switzerland  Camille-Sophie Brès  Photonic Systems Laboratory (PHOSL), Swiss
Federal Institute of Technology Lausanne (EPFL), CH-1015 Lausanne, Switzerland
Marco Liscidini  Dipartimento di Fisica, Università di Pavia, via Bassi 6,
I-27100 Pavia, Italy  Dario Gerace  Dipartimento di Fisica, Università di
Pavia, via Bassi 6, I-27100 Pavia, Italy

###### Abstract

Achieving the regime of single-photon nonlinearities in photonic devices just
exploiting the intrinsic high-order susceptibilities of conventional materials
would open the door to practical semiconductor-based quantum photonic
technologies. Here we show that this regime can be achieved in a triply
resonant integrated photonic device made of two coupled ring resonators,
without necessarily requiring low volume confinement, in a material platform
displaying an intrinsic third-order nonlinearity. By strongly driving one of
the three resonances of the system, a weak coherent probe at one of the others
results in a strongly suppressed two-photon probability at the output,
evidenced by antibunched second-order correlation function at zero-time delay
under continuous wave driving.

##  I Introduction

Reaching the ultimate limit of single-photon nonlinearity in conventional
semiconductor materials is the key to a number of applications in quantum
technology based on integrated photonic devices.[1] The regime of single-
photon nonlinearity would allow to control the response to an external
driving, e.g., through the presence of a single photon in a quantum optical
device[2, 3]. The extreme realization of such a single-photon switch occurs in
the so-called “photon blockade” (PB) regime, originally proposed for a cavity
enhanced optical nonlinearity[4]: the presence of two radiation quanta within
the resonator produces a nonlinear shift that is larger than the cavity mode
linewidth, which prevents a second photon from entering the resonator until
the first one has been released. The effect was first observed by measuring
antibunched second-order correlation of transmitted photons from an optical
cavity containing streams of single atoms in strong coupling with the cavity
field,[5] and later in semiconductor cavities coupled to single quantum dots
[6, 7]. More recently, signatures of antibunching in the transmitted radiation
have also been probed in polariton-confined microcavity samples [8, 9]. All
these demonstrations of PB rely on the presence of a strong dipolar
interaction of an optically active transition and a tightly confined
electromagnetic field. However, it would also be desirable to achieve the
regime of strong photon nonlinearities in purely passive materials, only
exploiting the intrinsic bulk nonlinear susceptibilities,[10, 11] although no
significant experimental result has been reported, to date. These approaches
typically require engineering low-volume cavities with large quality
factors,[12] such that the resonance shift associated with the effective
nonlinear response is of the order of its linewidth [10]. We hereby explore a
different approach to induce a single-photon nonlinearity on a weak probe
beam. The scheme relies on a strong coherent driving of a triply resonant
nonlinear system using an intense pump beam, in which enhanced four-wave
mixing occurs if the three resonances are equally spaced in frequency, a
condition easily met in whispering gallery mode resonators. At difference with
previous PB proposals, here we do not require the resonant device to be made
of small volume cavities. In fact, a possible experiment is represented in
Fig. 1: an asymmetric photonic molecule made of two ring resonators with
different radii results in an inhomogeneous frequency comb with at least three
equally-spaced resonances that can be spectrally isolated from the other
nearly resonant modes. A similar device has been proposed for enhanced
classical nonlinear processes,[13, 14] or the generation of squeezed
light.[15]

We show that this device can implement an effective second-order nonlinear
response on a weak probe, whose strength depends on the intensity of the pump,
and it thus allows one to produce PB on the probe according to the very same
mechanism studied before.[11] A similar model has been previously proposed for
dipolaritons in a planar microcavity[16]. Our results apply instead to passive
nonlinear devices, which can be implemented via integrated photonic
technologies in material platforms that not only possess a strong Kerr
nonlinearity – such as Lithium Niobate (LiN),[17], Silicon Nitride (SiN) and
Silicon-Rich Nitride (SRN),[18, 19, 20], Silicon Carbide (SiC)[21], and III-V
semiconductors[22, 23, 24] – but are also operated in their transparency
region, with negligible nonlinear absorption. This is crucial to avoid
possible heating and saturation of the pump, which would prevent reaching the
strong single-photon nonlinear regime. In addition, our theoretical proposal
acquires significant experimental relevance in light of the recent
demonstrations of strong coupling[25] and strong nonlinear coupling[26] in
ring resonators.

Figure 1: (a) Sketch of the proposed experiment: a single-mode waveguide,
driven by a strong and weak coherent field, the Pump and Seed sources
respectively, is evanescently coupled to an asymmetric photonic molecule with
Kerr-type nonlinearity. (b) Distribution of resonances and transitions (curly
arrows) involved in the effective $\chi^{(2)}$ nonlinear process.

##  II Model and results

The device hereby proposed is schematically depicted in Fig. 1(a). It consists
of a photonic molecule realized by two evanescently coupled ring resonators in
the “snowman” configuration, the larger having radius $R$ and the smaller
$R/2$. The two rings are evanescently side-coupled to a ridge waveguide, which
serves as an input-output channel to either inject excitations in the systems
or probe their quantum dynamics via output correlation measurements. In our
picture, we assume the “pump”, “seed”, and “idler” resonance frequencies,
respectively $\omega_{p,\,s,\,i}$, to fulfill the four-wave mixing condition
$2\omega_{s}=\omega_{p}+\omega_{i}$, and to correspond to those obtained by
accounting for any self- and cross-phase modulation correction induced by the
presence of a strong continuous-wave (CW) driving field at frequency
$\Omega_{p}$, which is used to induce a significant photon occupation of the
resonance at $\omega_{p}$. In this respect, the snowman configuration allows
for both spectrally isolating these three resonances from the rest of the
comb, and to realize, at the same time, a subset of resonances fulfilling
four-wave mixing condition, as detailed in appendix A. In fact, if the system
is driven only by the Pump CW source, since the resonances of the coupled
rings are not equally spaced, the $\Omega_{p}$ source will not induce any
significant population transfer from the $\omega_{p}$ resonance to nearby
ones. It is then reasonable to assume that in the absence of the Seed driving
the $p$-resonance is the only one macroscopically populated, with a steady-
state characterized by a mean number of photons $|\alpha_{p}|^{2}$
proportional to the input power $P_{p}$. Hence, starting from this stationary
condition for the $p$-resonance, the presence of a second weak CW source at
$\Omega_{s}\neq\Omega_{p}$ will not induce any significant depletion of the
$\omega_{p}$ resonance. Nevertheless, for $\Omega_{s}\sim\omega_{s}$, due to
the resonant condition, as schematically represented in Fig. 1(b), such CW
source will induce a net population of the seed resonance. We will show that
the latter undergoes PB, by characterizing the output statistics at the seed
frequency (assuming the pump photons to be fully filtered out).
In the rotating frame defined by the two CW source frequencies, for $P_{p}\gg
P_{s}$, the system Hamiltonian reads (see appendix B and appendix D for
derivation details):

$\hat{\mathcal{H}}=\hbar\Delta(\hat{a}_{s}^{\dagger}\hat{a}_{s}+2\hat{a}_{i}^{\dagger}\hat{a}_{i})+i\hbar\sqrt{\gamma}(\,F_{s}\,\hat{a}_{s}^{\dagger}-F^{*}_{s}\,\hat{a}_{s})+\hat{\mathcal{H}}_{nl},$
(1)

with $\hat{a}_{\sigma}$ and $\hat{a}_{\sigma}^{\dagger}$ representing
annihilation and creation (bosonic) operators at frequency $\omega_{\sigma}$
($\sigma=s,i$), with $\Delta=\omega_{s}-\Omega_{s}$, $F_{s}$ and $\gamma$
denoting the frequency detuning, the field amplitude of the CW source at
$\Omega_{s}$, and molecule-waveguide input coupling rate (assumed to be equal
for the three modes), respectively. Finally, in Eq. (9) the nonlinear
Hamiltonian reads

$\hat{\mathcal{H}}_{nl}=\hbar\tilde{g}_{nl}(\hat{a}_{s}^{\dagger\,2}\hat{a}_{i}+\hat{a}_{i}^{\dagger}\hat{a}_{s}^{2}),$
(2)

which describes the effective $\chi^{(2)}$ nonlinear interaction mediated by
the $p-$resonance, where $\tilde{g}_{nl}=g_{nl}|\alpha_{p}|$ denotes the
effective nonlinear coupling, dependent on both the bare
$\chi^{(3)}$-nonlinearity $g_{nl}$ and the pump mean-occupation via
$|\alpha_{p}|$. To account for optical lossess, we address the convergence
towards the stationary condition induced by the CW field at $\Omega_{s}$ by
means of a master equation in the Lindblad form [27], solving for the unique
steady-state solution for the system density matrix (see e.g. appendix C for a
discussion about uniqueness), $\rho_{ss}=\lim_{t\to+\infty}\rho(t)$, where

$\frac{d}{dt}\rho(t)=-\frac{i}{\hbar}\left[\hat{\mathcal{H}},\,\rho(t)\right]+\sum_{\sigma=s,i}\Gamma_{\sigma}\left[\hat{a}_{\sigma}\rho(t)\hat{a}_{\sigma}^{\dagger}-\frac{1}{2}\left\\{\hat{a}_{\sigma}^{\dagger}\hat{a}_{\sigma},\rho(t)\right\\}\right]$
(3)

with $\Gamma_{\sigma}$ denoting the total linewidth of the $\sigma$-resonance.
One can express $\Gamma_{\sigma}=\Gamma^{(0)}_{\sigma}+\gamma$, with
$\Gamma^{(0)}_{\sigma}$ being the intrinsic resonance linewidth. In the
following, we will assume, for the sake of clarity and without loss of
generality, that $\Gamma^{(0)}_{s}=\gamma$ (critical coupling condition), so
that $\Gamma_{s}=2\Gamma^{(0)}_{s}$.

The observables of interest are the steady-state average occupation of the
seed resonance,
$\langle\hat{n}_{s}\rangle\equiv\langle\hat{a}^{\dagger}_{s}\hat{a}_{s}\rangle=\mbox{Tr}[\hat{a}^{\dagger}_{s}\hat{a}_{s}\rho_{ss}]$,
and the equal time steady-state second-order autocorrelation function, that is
$g^{(2)}(0)$, whose expression reads

$g^{(2)}(0)=\frac{\langle\hat{a}^{\dagger}_{s}\hat{a}^{\dagger}_{s}\hat{a}_{s}\hat{a}_{s}\rangle}{\langle\hat{a}^{\dagger}_{s}\hat{a}_{s}\rangle\langle\hat{a}^{\dagger}_{s}\hat{a}_{s}\rangle}=\frac{\mbox{Tr}[\hat{a}^{\dagger}_{s}\hat{a}^{\dagger}_{s}\hat{a}_{s}\hat{a}_{s}\rho_{ss}]}{\mbox{Tr}[\hat{a}^{\dagger}_{s}\hat{a}_{s}\rho_{ss}]^{2}}.$
(4)

Numerical results obtained for such quantities as a function of the
dimensionless nonlinearity,
$\tilde{g}_{nl}/\Gamma_{s}=g_{nl}|\alpha_{p}|/\Gamma_{s}$, are reported in
Fig. 2. Some further results and comments concerning the numerical simulations
are reported in appendix E
In the absence of nonlinearity, since the seed and idler resonances are
decoupled, the seed corresponds to the stationary state of a standard driven-
dissipative oscillator. This is confirmed by the behavior displayed by both
the mean occupation and the autocorrelation function, the former being
characterized by a standard Lorentzian lineshape, i.e.,
$\langle\hat{n}_{s}\rangle={\gamma|F_{s}|^{2}}/\left[{\Delta^{2}+\left(\Gamma_{s}/2\right)^{2}}\right]$
and $g^{(2)}(0)=1$. The latter condition is the fingerprint of the steady-
state being in a pure coherent state, whose properties are uniquely determined
by the balance between driving and losses. Note that in the present case, due
to the weak-pumping of the seed resonance, there is no appreciable resonance
splitting, in contrast to what has been reported in the strong-nonlinear
regime [26]. The latter is indeed evidenced at stronger seed driving, as
shown, e.g., in appendix E.
On increasing the effective nonlinearity, the system approaches a stationary
configuration characterized by a non-trivial degree of correlation, crossing
the threshold value $g^{(2)}(0)=0.5$ (i.e., corresponding to the $n=2$ Fock
state) for effective nonlinearities $\tilde{g}_{nl}/\Gamma_{s}>0.325$. We
notice that the average photon occupation of the seed mode is only slightly
affected by the nonlinearity (see inset), while the different autocorrelation
curves display a strong dependence on $\tilde{g}_{nl}/\Gamma_{s}$. Notably,
the pronounced dip in the vicinity of $\Delta=0$ is a signature of a strong
antibunching regime, induced by the effective $\chi^{(2)}$-nonlinearity. In
other words, our numerical results show that the effective
$\chi^{(2)}$-nonlinearity is responsible for a clear suppression of the
occupation of any state containing more than a single seed-photon.

Figure 2: (a): Behavior of the steady-state mean occupation $\langle\hat{n}_{s}\rangle$ as a function of $\Delta/\Gamma_{s}$. The inset shows the behavior of the maximum seed occupation $\mbox{Max}(\langle\hat{n}_{s}\rangle)$ as a function of $\tilde{g}_{nl}/\Gamma$. (b): Steady-state equal time autocorrelation $g^{(2)}(0)$ as a function of $\Delta/\Gamma_{s}$. The different curves correspond to the different values of $\tilde{g}_{nl}/\Gamma_{s}$ reported in panel $(\textbf{a})$. The other simulation parameters are: $\Gamma_{i}=\Gamma_{s}$, $\gamma/\Gamma_{s}=0.5$ and $F_{s}/\sqrt{\Gamma_{s}}=0.1$. Material |  $n$ |  $n_{2}$ |  $Q_{0}$ |  $R$ |  $\gamma_{nl}$ | $\tilde{g}_{nl}/\Gamma_{s}$
---|---|---|---|---|---|---
|  $\rm{(m^{2}W^{-1})}$ |  $\times 10^{6}$ |  $\rm{(m)}$ |  $\rm{(m^{-1}W^{-1})}$ |  $P_{p}=$0.1\text{\,}\mathrm{W}$$ |  $P_{p}=$1\text{\,}\mathrm{W}$$ |  $P_{p}=$10\text{\,}\mathrm{W}$$
$\mathrm{LiNbO_{3}}$ [17] |  $2.24$ |  $1.8\times 10^{-19}$ |  $4.4$ |  $80\times 10^{-6}$ |  $0.7$ |  $5.7\times 10^{-4}$ |  $1.8\times 10^{-3}$ |  $5.7\times 10^{-3}$
$\mathrm{SRN}$ [19] |  $3.1$ |  $2.8\times 10^{-17}$ |  $0.06$ |  $100\times 10^{-6}$ |  $550$ |  $7.3\times 10^{-4}$ |  $2.3\times 10^{-3}$ |  $7.3\times 10^{-3}$
$\mathrm{SiC}$ [21] |  $2.6$ |  $6.9\times 10^{-19}$ |  $1.1$ |  $50\times 10^{-6}$ |  $3.2$ |  $3.8\times 10^{-4}$ |  $1.2\times 10^{-3}$ |  $3.8\times 10^{-3}$
$\mathrm{GaP}$ [22] |  $3.3$ |  $1.1\times 10^{-17}$ |  $0.25$ |  $50\times 10^{-6}$ |  $250$ |  $4.1\times 10^{-3}$ |  $1.3\times 10^{-2}$ |  $4.1\times 10^{-2}$
$\mathrm{Si_{3}N_{4}}$ [18] |  $2.0$ |  $2.5\times 10^{-19}$ |  $37$ |  $115\times 10^{-6}$ |  $0.7$ |  $1.2\times 10^{-2}$ |  $3.9\times 10^{-2}$ |  $0.12$
$\mathrm{AlGaAs}$ [23] |  $3.3$ |  $2.6\times 10^{-17}$ |  $3.14$ |  $28\times 10^{-6}$ |  $390$ |  $0.23$ |  $0.73$ |  $2.3$
Table 1: State of the art parameters and estimated nonlinear coupling for
microring resonators in different material platforms, comparing linear ($n$)
and nonlinear refractive indices ($n_{2}$), intrinsic quality factor in the
resonators ($Q_{0}$), radius ($R$), and nonlinear parameter ($\gamma_{nl}$).
The nonlinear coupling, $\tilde{g}_{nl}$, is estimated for different values of
input pump power ($P_{p}$), and assuming a working frequency
$\omega_{s}=$193\text{\,}\mathrm{THz}$$, with $\Gamma_{s}=\omega_{s}/Q_{0}$.
The values of group velocity ($v_{g}$) and effective area ($A_{eff}$) required
to estimate $\gamma_{nl}$ were numerically calculated via eigenmode expansion
using Ansys Lumerical software.

##  III Discussion

We now consider the accessibility of the experiment described above by
integrated photonic technologies at the state of the art and in different
material platforms of current interest. First, we directly compare different
platforms typically used in nonlinear integrated photonics in Table 1. The
choice has been made based on those materials presenting significant nonlinear
response, low propagation losses, and absence of two-photon absorption in the
typical telecom band. This explains the exclusion of Silicon, for example.
Based on these figures of merit, the same platforms are commonly used for the
generation of optical frequency combs. From the reported values of the
(intrinsic) quality factor $Q_{0}=\omega_{s}/\Gamma^{(0)}_{s}$, we estimate
the coupling rate according to the following expression [26]
$g_{nl}/\Gamma_{s}=\hbar\omega_{s}v_{g}^{2}\gamma_{nl}/(2\pi R\Gamma_{s})$,
where $v_{g}$ is the group velocity and
$\gamma_{nl}=\frac{\omega_{s}n_{2}}{A_{eff}c}$ is the nonlinear parameter of
the waveguide related to its $\chi^{(3)}$ response, which accounts for the
enhancement effect due to spatial confinement, where $n_{2}$ is the nonlinear
refractive index, and the effective mode area calculated for the given ring
resonator is defined as
$A_{eff}=\left(\int\left|E\right|^{2}dS\right)^{2}/\int_{ring}\left|E\right|^{4}dS$,
in which the steady state electric field profile is integrated on a surface
domain encompassing the ring cross-section (the denominator integral is only
performed on the nonlinear material domain, i.e. the ring cross section). From
the above expressions, it is possible to estimate the nonlinear coupling rate
as:

$\frac{\tilde{g}_{nl}}{\Gamma_{s}}=\frac{g_{nl}}{\Gamma_{s}}\sqrt{\frac{P_{p}}{\hbar\omega_{s}\Gamma_{s}}}=\frac{\gamma_{nl}v_{g}^{2}}{2\pi
R}\sqrt{\frac{\hbar\omega_{s}P_{p}}{\Gamma_{s}^{3}}}$ (5)

We report in Table 1 realistic values for the relevant parameters appearing in
Eq. 5, as well as a few estimates for the effective nonlinerarity at different
pump powers ($P_{p}=0.1W,\,1W,\,10W$), for different material platforms. The
expected behavior for the effective nonlinearity in the same platforms as a
function of driving powers in the $1\mu W$ \- $10\,W$ range is shown in Fig.
3.

Figure 3: Behavior of the effective nonlinearity $\tilde{g}_{nl}/\Gamma_{s}$
(logscale) as a function of the pump power $P_{p}\,[mW]$ for the material
platforms reported in Table 1. The dashed horizontal line locates the
threshold value $\tilde{g}_{nl}/\Gamma_{s}=0.325$, leading to a
$g^{(2)}(0)=0.5$.

First of all, we note that all of the values listed in the Table 1 are within
3 orders of magnitude from the “nonclassical limit”
$\tilde{g}_{nl}/\Gamma_{s}=0.325$ (threshold value corresponding to the dashed
horizontal line in Fig. 3), whereas the highest calculated values largely
exceed the best result achieved in high-$Q$, $\chi^{(2)}$ resonators[28] even
at relatively low pump power. The threshold for single-photon antibunching at
$g^{(2)}(0)=0.5$ is clearly reached and overcome, for the input pump powers
considered here, in the AlGaAs platform [23], which successfully combines a
high nonlinearity with a tight temporal and spatial confinement. Further
improvements on the other material platforms are expected by reducing the
microring radii while keeping large quality factors. These results suggest
that single photon blockade in passive $\chi^{(3)}$ devices might
realistically be achieved in current fabrication technology, upon appropriate
engineering of the nonlinear device. More in general, a high figure of merit
is attained in III-V alloys [22, 24], for which bandgap engineering allows to
maximize the nonlinear index while avoiding two-photon absorption. In this
perspective, the recently developed frequency combs in bichromatic photonic
crystal resonators [29, 30, 24] represent a promising approach to obtain a
triplet of equally-spaced resonances, as compared to the ring resonator layout
described in this work, where the intrinsically low mode volume would
contribute to increase the nonlinear coupling. For completeness, an estimate
of the nonlinear parameters for such microresonators has been numerically
obtained from data reported in [24], resulting in, e.g., in
$\tilde{g}_{nl}/\Gamma_{s}\simeq 0.02$ (0.064) at $P_{p}=0.1$ W (1 W).
Furthermore, a similar approach could also be investigated for intrinsic
$\chi^{(2)}$ processes in doubly-resonant photonic crystal resonators [31].
The SiC [21] and SiN [18] platforms have recently gained popularity in
integrated photonics, thanks to their excellent optical and nonlinear
properties, and to CMOS compatibility. They also represent two promising
candidates for quantum technology applications relying on photon blockade,
although significant improvements in respectively $Q$ factor and nonlinearity
would be needed to reach a performance comparable to the one of III-V devices.
From a technological perspective, the performance of the SiN platform could be
dramatically improved by using non-stoichiometric depositions (SRN)[20], for
which improvements in the $n_{2}$ index exceeding two orders of magnitude have
been demonstrated [19], without incurring in nonlinear absorption. Last but
not least, the LiN thin-film platform [17] represents a promising candidate
for applications in integrated quantum- and nonlinear-photonics. While in this
case the nonlinearity is quite small as compared to III-V semiconductors, the
large second-order nonlinear response of this material might be exploited to
achieve more complex interactions involving both the second- and third-order
susceptibilities. As a final comment, let us point out that the pump power
used to calculate the nonlinear coupling might differ in practical devices,
and could be further increased in some of these materials before incurring in
the damage threshold. Furthermore, we note that the present treatment, done
for a CW pump, could be extended to the pulsed regime, thus allowing to reach
higher peak power and therefore increasing the effective nonlinearity
achievable.

## References

  * Wang  _et al._ [2020a] J. Wang, F. Sciarrino, A. Laing, and M. G. Thompson, Integrated photonic quantum technologies, Nature Photonics 14, 273 (2020a).
  * Chang  _et al._ [2007] D. E. Chang, A. S. Sørensen, E. A. Demler, and M. D. Lukin, A single-photon transistor using nanoscale surface plasmons, Nature Physics 3, 807 (2007), arXiv:0706.4335 .
  * Gerace  _et al._ [2009] D. Gerace, H. E. Türeci, A. Imamoğlu, V. Giovannetti, and R. Fazio, The quantum-optical Josephson interferometer, Nature Physics 5, 281 (2009).
  * Werner and Imamoğlu [1999] M. J. Werner and A. Imamoğlu, Photon-photon interactions in cavity electromagnetically induced transparency, Physical Review A 61, 011801 (1999), arXiv:9902005 [quant-ph] .
  * Birnbaum  _et al._ [2005] K. M. Birnbaum, A. Boca, R. Miller, A. D. Boozer, T. E. Northup, and H. J. Kimble, Photon blockade in an optical cavity with one trapped atom, Nature 436, 87 (2005).
  * Faraon  _et al._ [2008] A. Faraon, I. Fushman, D. Englund, N. Stoltz, P. Petroff, and J. Vučković, Coherent generation of non-classical light on a chip via photon-induced tunnelling and blockade, Nature Physics 4, 859 (2008), arXiv:0804.2740 .
  * Reinhard  _et al._ [2012] A. Reinhard, T. Volz, M. Winger, A. Badolato, K. J. Hennessy, E. L. Hu, and A. Imamoğlu, Strongly correlated photons on a chip, Nature Photonics 6, 93 (2012).
  * Delteil  _et al._ [2019] A. Delteil, T. Fink, A. Schade, S. Höfling, C. Schneider, and A. İmamoğlu, Towards polariton blockade of confined exciton–polaritons, Nature Materials 18, 219 (2019).
  * Muñoz-Matutano  _et al._ [2019] G. Muñoz-Matutano, A. Wood, M. Johnsson, X. Vidal, B. Q. Baragiola, A. Reinhard, A. Lemaître, J. Bloch, A. Amo, G. Nogues, B. Besga, M. Richard, and T. Volz, Emergence of quantum correlations from interacting fibre-cavity polaritons, Nature Materials 18, 213 (2019).
  * Ferretti and Gerace [2012] S. Ferretti and D. Gerace, Single-photon nonlinear optics with Kerr-type nanostructured materials, Physical Review B 85, 033303 (2012), arXiv:1201.5072 .
  * Majumdar and Gerace [2013] A. Majumdar and D. Gerace, Single-photon blockade in doubly resonant nanocavities with second-order nonlinearity, Physical Review B 87, 235319 (2013), arXiv:1306.2183 .
  * Choi  _et al._ [2017] H. Choi, M. Heuck, and D. Englund, Self-Similar Nanocavity Design with Ultrasmall Mode Volume for Single-Photon Nonlinearities, Physical Review Letters 118, 223605 (2017).
  * Gentry  _et al._ [2014] C. M. Gentry, X. Zeng, and M. A. Popović, Tunable coupled-mode dispersion compensation and its application to on-chip resonant four-wave mixing, Optics Letters 39, 5689 (2014), arXiv:1406.2750 .
  * Heuck  _et al._ [2019] M. Heuck, J. G. Koefoed, J. B. Christensen, Y. Ding, L. H. Frandsen, K. Rottwitt, and L. K. Oxenløwe, Unidirectional frequency conversion in microring resonators for on-chip frequency-multiplexed single-photon sources, New Journal of Physics 21, 033037 (2019), arXiv:1811.11741 .
  * Zhang  _et al._ [2021] Y. Zhang, M. Menotti, K. Tan, V. D. Vaidya, D. H. Mahler, L. G. Helt, L. Zatti, M. Liscidini, B. Morrison, and Z. Vernon, Squeezed light from a nanophotonic molecule, Nature Communications 12, 2233 (2021), arXiv:2001.09474 .
  * Kyriienko and Liew [2014] O. Kyriienko and T. C. H. Liew, Triggered single-photon emitters based on stimulated parametric scattering in weakly nonlinear systems, Physical Review A 90, 063805 (2014).
  * He  _et al._ [2019] Y. He, Q.-F. Yang, J. Ling, R. Luo, H. Liang, M. Li, B. Shen, H. Wang, K. Vahala, and Q. Lin, Self-starting bi-chromatic LiNbO3 soliton microcomb, Optica 6, 1138 (2019).
  * Ji  _et al._ [2017] X. Ji, F. A. S. Barbosa, S. P. Roberts, A. Dutt, J. Cardenas, Y. Okawachi, A. Bryant, A. L. Gaeta, and M. Lipson, Ultra-low-loss on-chip resonators with sub-milliwatt parametric oscillation threshold, Optica 4, 619 (2017).
  * Wang  _et al._ [2015] T. Wang, D. K. T. Ng, S.-K. Ng, Y.-T. Toh, A. K. L. Chee, G. F. R. Chen, Q. Wang, and D. T. H. Tan, Supercontinuum generation in bandgap engineered, back-end CMOS compatible silicon rich nitride waveguides, Laser & Photonics Reviews 9, 498 (2015).
  * Bucio  _et al._ [2020] T. D. Bucio, C. Lacava, M. Clementi, J. Faneca, I. Skandalos, A. Baldycheva, M. Galli, K. Debnath, P. Petropoulos, and F. Gardes, Silicon Nitride Photonics for the Near-Infrared, IEEE Journal of Selected Topics in Quantum Electronics 26, 1 (2020).
  * Guidry  _et al._ [2020] M. A. Guidry, K. Y. Yang, D. M. Lukin, A. Markosyan, J. Yang, M. M. Fejer, and J. Vučković, Optical parametric oscillation in silicon carbide nanophotonics, Optica 7, 1139 (2020), arXiv:2004.13958 .
  * Wilson  _et al._ [2020] D. J. Wilson, K. Schneider, S. Hönl, M. Anderson, Y. Baumgartner, L. Czornomaz, T. J. Kippenberg, and P. Seidler, Integrated gallium phosphide nonlinear photonics, Nature Photonics 14, 57 (2020), arXiv:1808.03554 .
  * Xie  _et al._ [2020] W. Xie, L. Chang, H. Shu, J. C. Norman, J. D. Peters, X. Wang, and J. E. Bowers, Ultrahigh-Q AlGaAs-on-insulator microresonators for integrated nonlinear photonics, Optics Express 28, 32894 (2020), arXiv:2004.14537 .
  * Marty  _et al._ [2021] G. Marty, S. Combrié, F. Raineri, and A. De Rossi, Photonic crystal optical parametric oscillator, Nature Photonics 15, 53 (2021), arXiv:2011.13289 .
  * Guo  _et al._ [2016] X. Guo, C.-L. Zou, H. Jung, and H. X. Tang, On-Chip Strong Coupling and Efficient Frequency Conversion between Telecom and Visible Optical Modes, Physical Review Letters 117, 123902 (2016).
  * Ramelow  _et al._ [2019] S. Ramelow, A. Farsi, Z. Vernon, S. Clemmen, X. Ji, J. E. Sipe, M. Liscidini, M. Lipson, and A. L. Gaeta, Strong Nonlinear Coupling in a Si3N4 Ring Resonator, Physical Review Letters 122, 153906 (2019).
  * Breuer and Petruccione [2002] H. Breuer and F. Petruccione,  _The Theory of Open Quantum Systems_ (Oxford University Press, 2002).
  * Lu  _et al._ [2020] J. Lu, M. Li, C.-L. Zou, A. Al Sayem, and H. X. Tang, Toward 1lithium niobate microring resonators, Optica 7, 1654 (2020).
  * Alpeggiani  _et al._ [2015] F. Alpeggiani, L. C. Andreani, and D. Gerace, Effective bichromatic potential for ultra-high q-factor photonic crystal slab cavities, Applied Physics Letters 107, 261110 (2015), https://doi.org/10.1063/1.4938395 .
  * Clementi  _et al._ [2020] M. Clementi, A. Barone, T. Fromherz, D. Gerace, and M. Galli, Selective tuning of optical modes in a silicon comb-like photonic crystal cavity, Nanophotonics 9, 205 (2020).
  * Wang  _et al._ [2020b] J. Wang, M. Clementi, M. Minkov, A. Barone, J.-F. Carlin, N. Grandjean, D. Gerace, S. Fan, M. Galli, and R. Houdré, Doubly resonant second-harmonic generation of a vortex beam from a bound state in the continuum, Optica 7, 1126 (2020b).
  * Nigro [2019] D. Nigro, On the uniqueness of the steady-state solution of the lindblad–gorini–kossakowski–sudarshan equation, Journal of Statistical Mechanics: Theory and Experiment 2019, 043202 (2019).
  * Spohn [1976] H. Spohn, Approach to equilibrium for completely positive dynamical semigroups of n-level systems, Reports on Mathematical Physics 10, 189 (1976).
  * Drummond and Walls [1980] P. D. Drummond and D. F. Walls, Quantum theory of optical bistability. i. nonlinear polarisability model, Journal of Physics A: Mathematical and General 13, 725 (1980).
  * Meaney  _et al._ [2014] C. H. Meaney, H. Nha, T. Duty, and G. J. Milburn, Quantum and classical nonlinear dynamics in a microwave cavity, EPJ Quantum Technology 1, 10.1140/epjqt7 (2014).
  * Bartolo  _et al._ [2016] N. Bartolo, F. Minganti, W. Casteels, and C. Ciuti, Exact steady state of a kerr resonator with one- and two-photon driving and dissipation: Controllable wigner-function multimodality and dissipative phase transitions, Phys. Rev. A 94, 033841 (2016).

##  Appendix A Resonance spectrum of the photonic molecule

In this section, we discuss the motivation for using a device consisting of
two coupled ring resonators, such as the photonic molecule considered in the
main text, rather than a simpler one consisting of a single ring resonator
coupled to a waveguide, in light of engineering strong nonlinear interactions
between two given resonances of the structure, mediated by the presence of a
macroscopic population of a third resonance.
The starting point of our discussion concerns the distribution of resonances
in a ring resonator. In an ideal circular (i.e. rotationally invariant)
structure, the frequency distribution of azimuthal modes is uniquely
determined by $(i)$ the effective index $n_{\rm eff}$ of the waveguide, and
$(ii)$ the geometry of the resonator, that is, its radius $R$. By requiring
the optical path difference traveled by light along one round-trip of the
resonator, that is $n_{\rm eff}2\pi R$, to be a multiple of the wavelength,
and the effect of dispersion, it is straightforward to obtain that resonances
are distributed according to:

$\omega_{m}=m\delta\omega=m\frac{c}{n_{\rm eff}R},$ (6)

with the azimuthal number $m$ being an integer used to label the set of
resonances $\\{\omega_{m}\\}$, and with $\delta\omega$ denoting the Free
Spectral Range (FSR) of the ring. In such a scenario, due to the equal spacing
between consecutive resonances, the maximum population that can be stored
within a single resonance is severely limited by four-wave mixing processes
converting two photons at the frequency of interest, let us say
$\omega=\omega_{p}$, into a pair of photons corresponding to the adjacent
resonances, that is $(\omega_{p}+\delta\omega,\,\omega_{p}-\delta\omega)$.
In the case of rings with non-commensurate FSRs, a similar limitation would
occur also in the device reported in Fig. 4. On the contrary, if the radii of
the two resonators are commensurate, as sketched in panel (a) and panel (b) of
Fig. 5, some of the resonances (namely those at $\bar{\omega}_{p}$ and
$\bar{\omega}_{i}$) may share the same resonance frequency, and, due to the
interaction between the rings, get hybridized, generating new doublets of
resonances. With reference to panel c of Fig. 5, in the case of perfectly
aligned resonances, the position of such doublets is given by:

$\omega^{\pm}_{\sigma}=\bar{\omega}_{\sigma}\pm J,\,\,\sigma=p,\,i$ (7)

with $J$ denoting the coupling between the rings (which is assumed to be of
the same intensity for the resonances at $\bar{\omega}_{p}$ and
$\bar{\omega}_{i}$). In particular, as explicitly shown in Eq. 7, the
hybridization of such resonances leads to $(i)$ doublets which are less
sensitive to four-wave mixing processes (since adjacent resonances are not
evenly spaced) and $(ii)$ to the natural formation of isolated sets of three
equally spaced resonances from the comb of each single ring resonator. Indeed,
it is easy to verify that (see e.g. panel $(\textbf{c})$ of Fig. 5):

$\omega^{-}_{p}+\omega^{+}_{i}=2\omega_{s}\,\quad\mbox{and}\quad\omega^{+}_{p}+\omega^{-}_{i}=2\omega_{s}.$
(8)

Although in the main text we set $\omega_{p}=\omega^{-}_{p}$ and
$\omega_{i}=\omega^{+}_{i}$, the pair of resonances
$(\omega_{p},\,\omega_{i})$ that can be used to generate an effective strong
$\chi^{(2)}$ interaction may also correspond to
$(\omega^{+}_{p},\,\omega^{-}_{i})$. In particular, the analysis shown in the
main text would remain the same. In addition, we stress that for the purposes
of discussing effective $\chi^{(2)}$ nonlinearities induced by a strong
continuous-wave (CW) source driving one of the resonances in the $p$-doublet,
one can safely neglect nonlinear effects involving the other resonance. For
the sake of clarity, in order to show that we can arbitrarily choose one of
the two $p$-resonances in the vicinity of $\bar{\omega}_{p}$, let us make the
opposite choice with respect to the manuscript and assume
$\omega_{p}=\omega^{+}_{p}$ and $\omega_{i}=\omega^{-}_{i}$. In the absence of
any other CW source, since the resonances of the coupled rings are not equally
spaced, the $\Omega_{p}$ source acting in the proximity of $\omega_{p}$ will
not induce any significant population transfer from this resonance to nearby
ones. Therefore, it is reasonable to assume that in the stationary regime
induced by the $\Omega_{p}$ source such a resonance is the only one
macroscopically populated, with a state characterized by a mean amplitude
$\alpha_{p}$, such that $|\alpha_{p}|^{2}\propto P_{p}$, the latter being the
optical power of the strong CW source in the input waveguide. Once such a
stationary regime is reached, the presence of a second weak CW source at
$\Omega_{s}\neq\Omega_{p}$ will not induce any significant depletion of the
$\omega_{p}$ resonance. Nevertheless, for $\Omega_{s}\sim\omega_{s}$, such CW
source would induce a net population of the seed resonance. In parallel, the
nonlinear processes corresponding to the curly arrows in panel c of Fig. 5,
would scatter photons out of the seed resonance. Notice, however, that due to
the presence of the strong driving on the pump resonance at $\omega^{+}_{p}$,
nonlinear effects connecting states in the
$(\omega^{+}_{p},\,\omega_{s},\,\omega^{-}_{i})$ triplet will be more relevant
than those connecting states in the
$(\omega^{-}_{p},\,\omega_{s},\,\omega^{+}_{i})$ triplet. As a consequence,
these latter contributions, which are controlled by the intrinsic $\chi^{(3)}$
nonlinearity of the material and convert two seed photons into a
$(\omega^{-}_{p},\,\omega^{+}_{i})$ pair, can be neglected while describing
the dynamics induced by the weak CW source at $\Omega_{s}$. One comes to same
conclusions by considering $\omega_{p}=\omega^{-}_{p}$ and
$\omega_{i}=\omega^{+}_{i}$ as in the main text.

Figure 4: Proposed experimental setup to realize photon blockade in a photonic
integrated chip: a single-mode waveguide is evanescently coupled to an
asymmetric photonic molecule (coupled ring resonators) with Kerr-type
nonlinearity.  Figure 5: Distribution of the resonances of (a) the uncoupled
small ring (radius $R/2$), (b) the uncoupled large ring (radius $R$) and of
(c) the photonic molecule, with $J$ denoting the coupling strength between the
rings. For $J=0$, see panel (a) and panel (b), the resonances in each ring are
evenly spaced in frequency and due to the ratio between the radii those at
$\bar{\omega}_{p}$ and $\bar{\omega}_{i}$ are aligned. For $J\neq 0$, such
aligned resonances give rise to the doublets $(\omega^{-}_{p},\omega^{+}_{p})$
and $(\omega^{-}_{i},\omega^{+}_{i})$.

##  Appendix B Dynamical model of the photonic molecule

In the absence of any source term and any coupling to the surrounding
environment, that is $F_{s}=F_{p}=0$ and
$\Gamma_{\sigma}=0(\sigma=p,\,s,\,i)$, we assume the three resonances of the
photonic molecule to be governed by the following Hamiltonian operator

$\hat{\mathcal{H}}=\sum_{\sigma=p,\,s,\,i}\hbar\tilde{\omega}_{\sigma}\hat{a}_{\sigma}^{\dagger}\hat{a}_{\sigma}+\sum_{k,\,l,\,m,\,n=p,\,s,\,i}\left(\hbar
g_{klmn}\hat{a}_{k}^{\dagger}\hat{a}_{l}^{\dagger}\hat{a}_{m}\hat{a}_{n}+h.c.\right).$
(9)

In our scheme, $\tilde{\omega}_{\sigma}$ denotes the energy of the
$\sigma-$resonance of the cold system, that is in the absence of driving
terms. The $\omega_{\sigma}$ reported in the main text describe instead the
loaded-system resonances. The coupling constant $g_{klmn}$ accounts for
nonlinear effects that convert photon pairs occupying the $m$ and $n$ modes to
photon pairs in the modes $k$ and $l$, and it is proportinal to the component
of the susceptibility tensor $\chi^{(3)}_{klmn}$ of the material used to made
the ring resonators. The operators $\hat{a}_{\sigma}^{\dagger}$ and
$\hat{a}_{\sigma}$ correspond to creation and annihilation bosonic operators
of the $\sigma-$resonance, which satisfy the canonical commutation rules:

$[\hat{a}_{\alpha},\,\hat{a}_{\beta}]=0,\quad[\hat{a}_{\alpha},\,\hat{a}^{\dagger}_{\beta}]=\delta_{\alpha\,\beta},\quad\,\alpha,\beta=p,\,s,\,i$
(10)

with $\delta_{\alpha\,\beta}$ denoting the Kronecker delta, i.e.
$\delta_{\alpha\,\beta}=1\,$ if and only if $\alpha=\beta$.
In this framework, once the continuous-wave sources are turned on, some extra
terms must be included into the model, and the total Hamiltonian, which
becomes explicitly time-dependent, reads:

$\hat{\mathcal{H}}_{tot}=\hat{\mathcal{H}}+\hat{\mathcal{H}}_{F},$ (11)

with

$\hat{\mathcal{H}}_{F}=i\hbar\sqrt{\gamma}(\,F_{p}\,\hat{a}_{p}^{\dagger}e^{-i\,\Omega_{p}\,t}-F^{*}_{p}\,\hat{a}_{p}e^{i\,\Omega_{p}\,t})+i\hbar\sqrt{\gamma}(\,F_{s}\,\hat{a}_{s}^{\dagger}e^{-i\,\Omega_{s}\,t}-F^{*}_{s}\,\hat{a}_{s}e^{i\,\Omega_{s}\,t}),$
(12)

where the first term describes the driving of the $p$-resonance by a
monochromatic source at frequency $\Omega_{p}$, and the latter accounts for
the driving of the seed resonance with a monochromatic source at frequency
$\Omega_{s}$. In particular, for the sake of clarity, it is worth noting that
since $\gamma$ has the dimensions of a rate, that is
$\left[\gamma\right]=s^{-1}$, and since
$[\sqrt{\gamma}F_{p}]=[\omega_{\sigma}]=s^{-1}$, the field amplitudes
$F_{\sigma}$ do have the dimensions of the square root of a rate, i.e.
$\left[F_{\sigma}\right]=s^{-1/2}$.
In a realistic setup, photons can escape from the resonators. As a
consequence, the Hamiltonian model in Eq. 9 is not sufficient to describe the
time evolution of the photonic system, and one needs to introduce losses in
the theoretical description. A possible strategy is to use Lindblad formalism,
that is to cast the time evolution of the system density operator $\rho(t)$ in
the following form

$\frac{d}{dt}\rho(t)=\mathcal{L}[\rho(t)]\equiv-i\left[\hat{\mathcal{H}}_{tot}/\hbar,\,\rho(t)\right]+\sum_{\sigma=p,\,s,\,i}\Gamma_{\sigma}\,\mathbb{D}[\hat{a}_{\sigma};\,\rho(t)],$
(13)

with $\mathcal{L}$ being the Lindbladian (super-)operator, and with

$\mathbb{D}[\hat{L};\,\rho(t)]\equiv\hat{L}\rho(t)\hat{L}^{\dagger}-\frac{1}{2}\left\\{\hat{L}^{\dagger}\hat{L},\,\rho(t)\right\\}$
(14)

denoting the dissipator associated to a decay channel governed by the Lindblad
operator $\hat{L}$. In particular, in the present case we assume two-photon
losses to be negligible and we only consider non-unitary processes associated
to one-photon losses, that is those induced by a linear coupling between the
resonances and the surrounding environment. Such effects are captured by
Lindblad operators proportional to annihilation operators of the molecule
modes. In particular, in our model the parameters
$\Gamma_{\sigma}=\Gamma^{(0)}_{\sigma}+\gamma$ account for the total decay
rate of the resonances, being $\Gamma^{(0)}_{\sigma}$ the resonance linewidth
in the absence of coupling with the waveguide.

##  Appendix C Steady-state configuration and driving protocols

In the continuous-wave regime, the probability distributions of interest are
those corresponding to the steady-state configurations. Such configurations,
that can be formally defined as follows

$\lim_{t\to+\infty}\rho(t)=\rho^{(n)}_{ss}\quad n=1,\cdots,\,N_{ss},$ (15)

and that correspond to the $N_{ss}$ configurations annihilated by the
Lindbladian, that is

$\frac{d}{dt}\rho^{(n)}_{ss}=0=\mathcal{L}[\rho^{(n)}_{ss}],$ (16)

are those states approached by the system when there is a perfect balance
between unitary processes and non-unitary ones accounted by the dissipator. In
particular, once such configurations are determined, one has access to
expectation values by simply computing the following quantities:

$\langle\hat{O}\rangle^{(n)}=\mbox{Tr}[\hat{O}\rho^{(n)}_{ss}],$ (17)

with $\hat{O}$ being an operator and with “$\mbox{Tr}[A]$” denoting the trace
operation over the quantity $A$.
As specified discussed in the main text, our final goal is to characterize the
stationary regime reached in the presence of a strong driving on the
$p$-resonance and a weak one acting on the seed. In terms of the Hamiltonian
parameters, such regime corresponds to $|F_{p}|^{2}\gg|F_{s}|^{2}$.
In general, one might legitimately wonder if such a regime defines a unique
equilibrium distribution. In other words, one might wonder if different
driving protocols corresponding to the same continuous-wave regime, that is
equivalent after some reference time $t=t_{c}$, may lead to different steady-
states. In order to clarify this point, please consider Fig. 6, where we
sketched two different protocols that consider different driving combinations
for $t<t_{c}$, but defining the same continuous-wave regime after for larger
times. In both cases, let us assume the system to be in the vacuum state at
$t=0$. In the “protocol A”, see Fig. 6(a), we first turn on the pump-driving
field at $t=t_{1}$ and then at later time we also drive the population of the
$s$-resonance. When using the “protocol B” we do consider the opposite
scenario, so that we first start to drive the seed and after $t=t_{c}$ we also
inject photons in the $p$-resonance. As a result, since the protocols are
different, it would not be surprising to have different configurations at
$t=t_{c}$. Nevertheless, what should we expect for $t\gg t_{c}$?

(a) $protocol\,A$

(b) $protocol\,B$

Figure 6: Different driving protocols leading to the same continuous-wave
regime for $t>t_{c}$. Fig. 6(a): the pump field (with intensity
$\propto|F_{p}|^{2}$) is turned on at $t_{1}$ and then at later time $t_{c}$
we start to drive also the population of the $s$-resonance (with intensity
$\propto|F_{s}|^{2}$). Fig. 6(b): we first turn on the driving of the
$s$-resonance at $t_{1}$ (with intensity $\propto|F_{s}|^{2}$) and then at
later time we start to drive also the population of the $p$-resonance (with a
field intensity $\propto|F_{p}|^{2}$).

The answer to this question is that the faith of the (open) quantum system
depends on the number of steady-states supported by the master equation for
$t>t_{c}$. In general, if different steady-states exist, then the future
evolution of the quantum system is defined by its initial configuration.
Indeed, similarly to what happens in the case of closed quantum systems, the
existence of many stationary regimes is related to the existence of at least
one (strong) $symmetry$ preserved by the Lindbladian. Let us denote such a
conserved quantity by $\hat{Q}$. If the two states generated by the two
driving protocols, let us say $\rho_{A}(t_{c})$ and $\rho_{B}(t_{c})$, are
characterised by two different values of such symmetry, that is
$q_{A}=\mbox{Tr}[\hat{Q}\,\rho_{A}(t_{c})]$ and
$q_{B}=\mbox{Tr}[\hat{Q}\,\rho_{B}(t_{c})]$, then due to the conservation of
$\hat{Q}$ one can conclude that

$\rho^{(A)}_{ss}\neq\rho^{(B)}_{ss}.$ (18)

On the other hand, if the steady-state is unique, the dynamics is _relaxing_ ,
that is

$\lim_{t\to+\infty}\rho(t)=\rho_{ss}\quad\forall\,\rho_{0}=\rho(0),$ (19)

or in other words, if the steady-state is unique, then any initial
configuration will converge to it if we wait a sufficiently long time.
Therefore, for our purposes, it becomes of primary interest to understand if
the master equation in Eq. 13 does support a unique steady-state in the
continuous-wave regime for $t>t_{c}$. In order to verify that this is indeed
the case, it is sufficient to consider the results reported in Ref. [32],
where it is shown that the presence of one-photon losses acting on each
resonance of a quantum system is a sufficient condition to conclude that the
master equation has a unique steady-state $\rho_{ss}$. Furthermore, by means
of the results derived in Ref. [33], one can prove that the dynamics encoded
into Eq. 13 is also _relaxing_.
In other words, these facts ensure that:

  * 1)

the stationary state for $t>t_{c}$, that is the steady-state $\rho_{ss}$, is
_unique_ ;

  * 2)

any initial configuration at $t=t_{c}$ will converge to $\rho_{ss}$ for
sufficiently long times;

  * 3)

the dynamics for $t<t_{c}$ is irrelevant for determining the long time
behavior of the system, in the sense that any deviation introduced by a
different choice in the driving protocol for $t<t_{c}$, for instance “protocol
A” versus “protocol B”, will eventually disappear for large enough time.

##  Appendix D Derivation of the effective two-resonance model

In this section we derive an effective model describing the dynamics of the
seed-idler subsystem within the full dynamical regime defined in the previous
Section. As a first step, we look for the physical configuration that
guarantees the dynamics of such a subsystem to be described by the master
equation reported in the main text. To this purpose, let us consider the
generic rotating frame defined by the following unitary operator:

$R(t)=e^{i(\tilde{\Omega}_{p}\hat{a}^{\dagger}_{p}\hat{a}_{p}+\tilde{\Omega}_{s}\hat{a}^{\dagger}_{s}\hat{a}_{s}+\tilde{\Omega}_{i}\hat{a}^{\dagger}_{i}\hat{a}_{i})t}\equiv\exp(i\hat{\Lambda}t)$
(20)

In this frame, the rotated density operator is given by

$\tilde{\rho}=R(t)\rho R^{\dagger}(t),$ (21)

and its time evolution is given by the following rotated master equation

$\frac{d}{dt}\tilde{\rho}=i\,R(t)[\hat{\Lambda},\,\rho]R^{\dagger}(t)+R(t)\frac{d}{dt}\rho(t)R^{\dagger}(t)=iR(t)[\hat{\Lambda},\,\rho]R^{\dagger}(t)+R(t)\mathcal{L}[\rho(t)]R^{\dagger}(t).$
(22)

In order to proceed further, we need an explicit expression for the rotated
bosonic operators, that is
$\tilde{a}_{\sigma}(t)=R(t)a_{\sigma}R^{\dagger}(t)$. By considering that
$\tilde{a}_{\sigma}(t)$ is characterized by the following time evolution

$\frac{d}{dt}\tilde{a}_{\sigma}(t)=R(t)\left[i\hat{\Lambda},\,\hat{a}_{\sigma}\right]R^{\dagger}(t)=-i\tilde{\Omega}_{\sigma}\,\tilde{a}_{\sigma}\rightarrow\tilde{a}_{\sigma}(t)=e^{-i\tilde{\Omega}_{\sigma}t}\tilde{a}_{\sigma}(0)=e^{-i\tilde{\Omega}_{\sigma}t}\hat{a}_{\sigma},$
(23)

one obtains that the rotated density operator evolves through the following
master equation

$\frac{d}{dt}\tilde{\rho}=-i[\hat{\mathcal{H}}_{tot}(t)-\hat{\Lambda},\,\tilde{\rho}]+\sum_{\sigma}\Gamma_{\sigma}\mathbb{D}[\hat{a}_{\sigma};\,\tilde{\rho}],$
(24)

where

$\begin{split}\hat{\mathcal{H}}_{tot}(t)-\hat{\Lambda}&=(\tilde{\omega}_{p}-\tilde{\Omega}_{p})\hat{a}^{\dagger}_{p}\hat{a}_{p}+(\tilde{\omega}_{s}-\tilde{\Omega}_{s})\hat{a}^{\dagger}_{s}\hat{a}_{s}+(\tilde{\omega}_{i}-\tilde{\Omega}_{i})\hat{a}^{\dagger}_{i}\hat{a}_{i}+\\\
&+i\sqrt{\gamma}(\,F_{p}\,\hat{a}_{p}^{\dagger}-F^{*}_{p}\,\hat{a}_{p})+i\sqrt{\gamma}(\,F_{s}\,\hat{a}_{s}^{\dagger}e^{i(\tilde{\Omega}_{s}-\,\Omega_{s})\,t}-F^{*}_{s}\,\hat{a}_{s}e^{-i(\tilde{\Omega}_{s}\,-\Omega_{s})\,t})+\\\
&+\sum_{k,\,l,\,m,\,n=p,\,s,\,i}\left(\hbar
g_{klmn}\hat{a}_{k}^{\dagger}\hat{a}_{l}^{\dagger}\hat{a}_{m}\hat{a}_{n}e^{i(\tilde{\Omega}_{k}+\tilde{\Omega}_{l}-\tilde{\Omega}_{m}-\tilde{\Omega}_{n})t}+h.c.\right).\end{split}$
(25)

Let us now consider the rotating frame defined by the following constraints:

$\tilde{\Omega}_{s}=\Omega_{s}\,,$ (26)

and

$\tilde{\Omega}_{i}=2\tilde{\Omega}_{s}-\tilde{\Omega}_{p}\,,$ (27)

the former corresponding to choose a reference frame rotating at $\Omega_{s}$
for the $s$-resonance, and the latter imposing energy conservation for the
nonlinear terms converting two seed photons into a pump-idler pair, i.e., the
processes considered in the main text. In this frame, by making a standard
rotating-wave approximation (which corresponds to neglecting any off-resonant
interaction term), the previous Hamiltonian model reduces to

$\begin{split}\hat{\mathcal{H}}_{tot}(t)-\hat{\Lambda}&\approx(\tilde{\omega}_{p}-\tilde{\Omega}_{p})\hat{a}^{\dagger}_{p}\hat{a}_{p}+(\tilde{\omega}_{s}-\Omega_{s})\hat{a}^{\dagger}_{s}\hat{a}_{s}+(\tilde{\omega}_{i}+\tilde{\Omega}_{p}-2\Omega_{s})\hat{a}^{\dagger}_{i}\hat{a}_{i}+\\\
&+i\sqrt{\gamma}(\,F_{p}\,\hat{a}_{p}^{\dagger}-F^{*}_{p}\,\hat{a}_{p})+i\sqrt{\gamma}(\,F_{s}\,\hat{a}_{s}^{\dagger}-F^{*}_{s}\,\hat{a}_{s})+\\\
&+\sum_{\sigma=p,\,s,\,i}\hbar
g_{\sigma\sigma\sigma\sigma}\hat{a}^{\dagger\,2}_{\sigma}\hat{a}^{2}_{\sigma}+\sum_{k\neq
l;k,l=p,\,s,\,i}\hbar
g_{klkl}\hat{a}^{\dagger}_{k}\hat{a}_{k}\hat{a}^{\dagger}_{l}\hat{a}_{l}+\hbar
g_{nl}(\hat{a}^{\dagger\,2}_{s}\hat{a}_{i}\hat{a}_{p}+h.c)+\hat{V}(t)\\\
\end{split}$ (28)

where we set $g_{ssip}\equiv g_{nl}$ for consistency with the notation used in
the main text.
Evidently, the model still depends on the $p$-resonance. In order to get rid
of any dynamical dependence on the $p$-operators, we apply the following
argument. According to the discussion made in the previous section, the
steady-state reached by the system in the continuous wave regime is
independent of the particular protocol used for driving the system, i.e.,
protocol “A” or “B”. Nevertheless, since we are looking for the stationary
regime defined by

$|F_{p}|^{2}\gg|F_{s}|^{2},$ (29)

we do believe that, in order to build some intuition about the steady-state
properties, it is convenient to consider a protocol of type “A”. If this is
the case, at $t=t_{1}$ the system is in vacuum, and by switching the $F_{p}$
field on, it will be driven towards the steady-state of the master equation
defined in Eq. 24 for $F_{s}=0$. We argue that, due to the structure of the
Hamiltonian, such a state will be a product state with the $p$-resonance
somehow populated, and with the seed-idler subsystem empty. Indeed, no term in
the Hamiltonian (due to the shifts between levels) is capable of transfering
population from the $p$-resonance to the others. In other words, with
reference to Fig. 6(a), for $t\approx t_{c}\gg t_{1}$ the system density
operator will be $\tilde{\rho}(t)=\rho_{p;\,ss}\otimes\rho^{vacuum}_{s,i}$,
with $\rho^{vacuum}_{s,i}$ denoting the vacuum state density operator of the
seed-idler subsystem. In this scenario, we have that $\rho_{p,\,ss}$ is the
solution of

$\tilde{\mathcal{L}}[\rho_{p;\,ss}]=-i[\hbar(\tilde{\omega}_{p}-\Omega_{p})\hat{a}_{p}^{\dagger}\hat{a}_{p}+\hbar
g_{pppp}\hat{a}_{p}^{\dagger\,2}\hat{a}_{p}^{2}+i\sqrt{\gamma}(\,F_{p}\,\hat{a}_{p}^{\dagger}-F^{*}_{p}\,\hat{a}_{p}),\,\rho_{p;\,ss}]+\Gamma_{p}\mathbb{D}[\hat{a}_{p};\,\rho_{p;\,ss}]=0\,.$
(30)

Hence, by switching the second driving on, since we have that
$|F_{s}|^{2}\ll|F_{p}|^{2}$, it is reasonable to assume that the $p$-resonance
will only be slightly affected by the photons injected in the seed-idler
subsystem. As a consequence, the steady-state approached in the continuous-
wave regime for $t\gg t_{c}$ can be approximated by a density operator that is
still in factorized form. More precisely, we do expect the system to approach
a configuration $\tilde{\rho}_{ss}\approx\rho_{p;\,ss}\otimes\rho_{s,i;\,ss}$,
with $\rho_{s,i;\,ss}$ corresponding to the steady-state density operator that
is solution of the following master equation:

$\frac{d}{dt}\rho_{s,\,i}=-i\left[\hat{\mathcal{H}}_{s,\,i},\,\rho_{s,\,i}\right]+\Gamma_{s}\mathbb{D}[\hat{a}_{s};\,\rho_{s,\,i}]+\Gamma_{i}\mathbb{D}[\hat{a}_{i};\,\rho_{s,\,i}],$
(31)

with

$\begin{split}\hat{\mathcal{H}}_{s,\,i}=&(\tilde{\omega}_{s}-\Omega_{s}+\delta_{XPM,\,p,\,s})\hat{a}^{\dagger}_{s}\hat{a}_{s}+(\tilde{\omega}_{i}+\Omega_{p}-2\Omega_{s}+\delta_{XPM,\,p,\,i})\hat{a}^{\dagger}_{i}\hat{a}_{i}+\\\
&+i\sqrt{\gamma}(\,F_{s}\,\hat{a}_{s}^{\dagger}-F^{*}_{s}\,\hat{a}_{s})+\hbar
g(\alpha_{p}\hat{a}^{\dagger\,2}_{s}\hat{a}_{i}+\eta^{*}_{p}\hat{a}^{\dagger}_{i}\hat{a}^{2}_{s})+\\\
&+\sum_{\sigma=\,s,\,i}\hbar
g_{\sigma\sigma\sigma\sigma}\hat{a}^{\dagger\,2}_{\sigma}\hat{a}^{2}_{\sigma}+\hbar
g_{sisi}\hat{a}^{\dagger}_{s}\hat{a}_{s}\hat{a}^{\dagger}_{i}\hat{a}_{i},\\\
\end{split}$ (32)

where $\delta_{XPM,\,p,\,\sigma}=g_{p\sigma
p\sigma}\langle\hat{a}^{\dagger}_{p}\hat{a}_{p}\rangle_{p,\,ss}$ denotes the
frequency shift related to the cross-phase modulation effects induced by the
$p$-resonance on the $\sigma$-resonance and with
$\alpha_{p}=\langle\hat{a}_{p}\rangle_{p,\,ss}$ denoting the field amplitude
in the $p$-resonance. Before going any further, it is worth spending some
words on the effective parameters entering the Hamiltonian model. On the one
hand, by construction,
$\langle\hat{a}^{\dagger}_{p}\hat{a}_{p}\rangle_{p,\,ss}$ is a real quantity
and there is no ambiguity in using its value in the Hamiltonian; on the one
other hand, the field amplitude $\alpha_{p}$ is, in general, complex, i.e.
$\alpha_{p}=|\alpha_{p}|e^{i\phi}$, with $|\alpha_{p}|$ and $e^{i\phi}$
denoting the absolute value and the phase of $\alpha_{p}$. Such a phase is
irrelevant for the dynamics. Indeed, such a phase can always be removed from
the master equation by simply redefining the $\hat{a}_{i}$ operator, i.e., by
using $\hat{a}_{i}\to e^{i\phi}\hat{a}_{i}$.
In particular, in the previous expressions,
$\langle\hat{O}\rangle_{p,\,ss}\equiv\mbox{Tr}[\hat{O}\rho_{p;\,ss}]$ denotes
the expectation value of the operator $\hat{O}$ on the state $\rho_{p;\,ss}$,
that, interestingly, have some closed expression in terms of hypergeometric
functions for the Lindblad evolution in Eq. 30, see e.g. [34, 35, 36].
If we set $\omega_{\sigma}=\tilde{\omega}_{\sigma}+\delta_{XPM,\,p,\,\sigma}$,
and suppose, without loss of generality, that

$\omega_{i}+\omega_{p}=2\omega_{s},$ (33)

we have that the Hamiltonian model in Eq.32, up to corrections accounting for
self-phase modulation and cross-phase modulation effects between the seed and
ilder resonances, reduces to the one reported in the main text. Notice,
however, that since such corrections are proportional to the bare coupling
constant $g_{ssss}$, $g_{iiii}$ and $g_{sisi}$ which are by nature orders of
magnitude smaller than the effective coupling $g_{nl}$, they can be safely
neglected.

##  Appendix E Steady-state observables for the seed-idler subsystem

In this section we report some details about the numerical methods used to
determine the steady-state configuration and show some numerical results
describing the behavior of the steady-state mean occupation of the seed
resonance for values of the parameters not displayed in the main text.

###  E.1 Convergence to the steady-state expectation values: finite-$N$
simulations

In the absence of losses and, in general, in the absence of any decay
mechanism, boson systems in second quantization are described by the infinite
dimensional space spanned by Fock states, i.e. number states. In the presence
of losses, since high occupation states decay in favour of lower number
states, one has that density operators can be approximate by means of a finite
number of Fock states. To the best of our knowledge, such a number cannot be
determined _a priori_ , and in order to characterize the system properties in
the stationary regime one needs to look at the scaling of observables while
increasing the number of states used in the representation of the Lindblad
operator. Results shown in the main text and in this supplementary file have
been obtained by adopting the following prescription: we truncate the Fock
basis of the idler resonance at $N_{max}$ excitations, and the seed one at
$2N_{max}$ excitations, with $N_{max}$ being a positive integer. The
motivation for choosing a different truncation on the two resonances is
related to the structure of the Hamiltonian. Indeed, since we considered only
the effects determined by a direct driving on the $\omega_{s}$ resonance and
since population is transferred to the idler resonance only via the nonlinear
interaction (which converts two seed photons into an idler photon), in
stationary conditions an idler state with $m$ photons is occupied if and only
if a seed state with $2m$ photons is populated too. In particular, by
proceeding this way one has that the number of states used to represent
density operators of the seed-idler subsystem is given by
$M(N_{max})\equiv(2N_{max}+1)(N_{max}+1)$, where the “+1” comes from counting
also the vacuum state. For any given $N_{max}$ and any point in the parameter
space, we then determined, via standard linear algebra algorithm, the unique
configuration annihilated by the Lindblad superoperator, that is the fixed
point of Eq.31 (that is by solving the linear system
$\mathcal{L}[\rho_{ss}]=0$).
An example of scaling observed while increasing $N_{max}$ is reported in Fig.
7. There, we show some numerical results obtained for the steady-state mean
occupation of the seed resonance, that is $\langle\hat{n}_{s}\rangle$, as a
function of the detuning $\Delta/\Gamma_{s}$, for different values of the
effective nonlinear interaction $\tilde{g}_{nl}/\sqrt{\Gamma}$ (specified in
each panel). In particular, results in Fig 7(a) and Fig 7(b) describe the
convergence to the steady-state expectation value for
$F_{s}/\sqrt{\Gamma_{s}}=1.0$ and $F_{s}/\sqrt{\Gamma_{s}}=2.5$. As suggested
by the inset plots in the panels describing the behavior of
$Max(\langle\hat{n}_{s}\rangle)$, due to the weaker driving rate, data for
$F_{s}/\sqrt{\Gamma_{s}}=1.0$ can be considered at convergence even at
$N_{max}\sim 5$. On the contrary, for $F_{s}/\sqrt{\Gamma_{s}}=2.5$, data can
be considered at convergence for $N_{max}\sim 8-9$.

(a) $F_{s}/\sqrt{\Gamma_{s}}=1.0$

(b) $F_{s}/\sqrt{\Gamma_{s}}=2.5$

Figure 7: Behavior of the steady-state mean occupation of the seed resonance
$\langle\hat{n}_{s}\rangle$ as a function of the detuning $\Delta/\Gamma_{s}$,
for different values of the effective nonlinear interaction
$\tilde{g}_{nl}/\sqrt{\Gamma}$ (specified in each panel) and for
$F_{s}/\sqrt{\Gamma_{s}}=1.0$ (7(a)) and for $F_{s}/\sqrt{\Gamma_{s}}=2.5$
(7(b)). The different curves in the main panels describe results obtained for
different values of $N_{max}$ reported in the top legend. In each panel, we
report also the behavior of $Max(\langle\hat{n}_{s}\rangle)$ as a function of
$N_{max}$.

###  E.2 Behavior of the seed mean occupation at increasing power

In this section we show some numerical results describing the behavior at
increasing $F_{s}$ of the mean occupation of the seed resonances. In
particular, in Fig. 8 we show the behavior of such quantity (data displayed
correspond to $N_{max}=10$) divided by $|F_{s}|^{2}/\Gamma_{s}$ (which in
general provides an estimate of the mean number of photons stored within a
linear driven-dissipative resonator). As it is possible to see, on the
increasing of the input power acting on the seed resonance
($\propto|F_{s}|^{2}$), the system displays a pronounced line-splitting. At
low driving, the mean-steady state occupation is characterized by a single
Lorentzian-peak, located at $\Delta=0$, i.e. $\Omega_{s}=\omega_{s}$. On the
increasing of the driving intensity of CW at $\Omega_{s}$, in agreement with
results previously appeared in literature, see e.g. [26], two peaks appear in
the seed response, which is a signature of the presence of a strong nonlinear
interaction acting in the seed-idler subsystem.

Figure 8: Behavior of the ratio between the steady-state mean occupation of
the seed resonance, that is $\langle\hat{n}_{s}\rangle$, and the normalized
driving intensity $|F_{s}|^{2}/\Gamma_{s}$, as a function of the detuning
$\Delta/\Gamma_{s}$. The different panels correspond to different values of
the effective effective nonlinearity $\tilde{g}_{nl}/\Gamma_{s}$.