# Fock lasers based on deep-strong coupling of light and matter

Nicholas Rivera1, Jamison Sloan2, Ido Kaminer3, and Marin Soljačić1,2
1Department of Physics, MIT, Cambridge, MA 02139, USA.
2Research Laboratory of Electronics, MIT, Cambridge, MA 02139, USA.
3Department of Electrical and Computer Engineering, Technion, Haifa 32000,
Israel.

###### Abstract

Light and matter can now interact in a regime where their coupling is stronger
than their bare energies. This deep-strong coupling (DSC) regime of quantum
electrodynamics promises to challenge many conventional assumptions about the
physics of light and matter. Here, we show how light and matter interactions
in this regime give rise to electromagnetic nonlinearities dramatically
different from those of naturally existing materials. Excitations in the DSC
regime act as photons with a linear energy spectrum up to a critical
excitation number, after which, the system suddenly becomes strongly
anharmonic, thus acting as an effective intensity-dependent nonlinearity of an
extremely high order. We show that this behavior allows for _N-photon
blockade_ (with $N\gg 1$), enabling qualitatively new kinds of quantum light
sources. For example, this nonlinearity forms the basis for a new type of gain
medium, which when integrated into a laser or maser, produces large Fock
states (rather than coherent states). Such Fock states could in principle have
photon numbers orders of magnitude larger than any realized previously, and
would be protected from dissipation by a new type of equilibrium between
nonlinear gain and linear loss. We discuss paths to experimental realization
of the effects described here.

Recent successes in the coupling of matter and light now make it possible to
realize regimes of light-matter interactions in which the coupling between
light and matter can be much stronger than in established optical
technologies. Because of the central role the physics of light and matter
plays in many fields, these new coupling regimes are being intensely explored.
One such example is the ultra-strong coupling regime, where the coupling
energy is within an order of magnitude of the bare energies of the light and
matter kockum2019ultrastrong . Such regimes promise to give rise to new
chemical processes hutchison2012modifying ; flick2018strong ;
ruggenthaler2018quantum ; thomas2019tilting , strong modifications of
transport and thermodynamic properties of materials orgiu2015conductivity ;
paravicini2019magneto , new phases of matter, as well as advances in quantum
simulators and quantum technologies more broadly kockum2019ultrastrong ;
forn2019ultrastrong .

Taking these ideas to their logical extreme is the so-called _deep-strong
coupling regime_ (DSC), where the strength of the coupling is _greater_ than
the bare energies of the light and matter. In the past few years, the first
experiments in this regime have emerged yoshihara2017superconducting ;
forn2017ultrastrong . Much of the interest in ultra-strong and deep-strong
coupling is focused on the properties of the ground state of either one or
many emitters coupled to a cavity mode, leading to many interesting new
phenomena such as light-matter decoupling de2014light ; rivera2019variational
, population collapses and revivals casanova2010deep , large Lamb shifts
leading to inversion of qubit energy levels yoshihara2017superconducting ;
yoshihara2018inversion , and renormalization of qubit energy levels by a
photonic continuum forn2017ultrastrong . Likely, many of the potential
applications of this regime have yet to be identified.

Here, we consider the opportunities afforded to us by the excited states of a
DSC system, which are important from the perspective of quantum and nonlinear
optics. For example, the emission of light in such systems probes the excited
states. First, we show that deep-strong coupling of a two-level system to a
resonant cavity leads to the formation of excitations (“photonic
quasiparticles”, which we refer to as “DSC photons”) with nonlinear properties
much different than those in any known system. Then, we analyze the coupling
of an emitter to this nonlinear photonic quasiparticle. We find that the
coupling of an excited two-level system to this nonlinear system enables a
phenomenon of $N$-photon blockade in which $N$ excitations can be populated,
but $N+1$ cannot. We show that a laser or maser based on stimulated emission
of DSC photons behaves fundamentally differently from a conventional maser or
laser. Specifically, this maser produces close approximations to Fock states
in its steady-state, rather than coherent states, as in conventional lasers.
They could have a few hundred photons, thus being orders of magnitude larger
than any Fock states realized thus far. Moreover, Fock states produced by this
mechanism are stable against dissipation as they arise from a new type of
equilibrium between nonlinear gain and linear loss. Our results may thus help
to address the long-standing problem in quantum science of generating Fock
states. Finally, we discuss how the concept developed here can be implemented
in superconducting qubit platforms.

##  I Nonlinear photonic quasiparticles based on deep strong light-matter
coupling

Figure 1: High-order nonlinearities in deep strong coupling of light and
matter. (a) Schematic of a two-level system coupled to a single resonator
mode, as in circuit or cavity QED. (b) Spectrum of the system from weak
($g=\tilde{g}/\omega\ll 1$) to deep-strong coupling ($g\gg 1$). Here,
$\lambda=0$. (c) Successive excitation energies for a single spin sector for
different coupling values. For $g\gg 1$ the excitation energies are constant,
as for a bare photon. At large photon number, they deviate rapidly and
nonlinearly from harmonicity, akin to a photon with a strongly intensity-
dependent nonlinearity.

Fundamental to our results is the spectrum of a two-level system (qubit)
interacting with a single-mode cavity (schematically illustrated in Fig. 1a),
which we review here irish2007generalized . The Hamiltonian, referred to as
the (generalized) Rabi Hamiltonian, is given by

$H_{\text{Rabi}}/\hbar=\frac{1}{2}\left(\omega_{0}\sigma_{z}+\lambda\sigma_{x}\right)+\omega
a^{\dagger}a+\tilde{g}\sigma_{x}(a+a^{\dagger}).$ (1)

Here, $\omega_{0}$ is the transition frequency of the two-level system,
$\sigma_{x,z}$ are the $x$ and $z$ Pauli matrices, $\omega$ is the cavity
frequency, $a^{(\dagger)}$ is the cavity annhilation (creation) operator, and
$\tilde{g}$ is the Rabi frequency. It will be convenient to non-dimensionalize
the coupling as $g=\tilde{g}/\omega$. We have also generalized the standard
Rabi Hamiltonian by including a term $\lambda\sigma_{x}$ which is relevant in
contexts of superconducting qubits with applied bias fluxes
yoshihara2017superconducting . For simplicity of presentation, we consider the
case of $\lambda=0$, which leads to approximately degenerate spin states (and
in which case the qubit frequency is $\omega_{0}$). In the SI, and in various
numerical results, we do consider the effect of a finite $\lambda$ term, which
yields the same qualitative conclusions.

While the Rabi Hamiltonian cannot be analytically diagonalized in general, an
approximate spectrum can be found for the regime $g\gg 1$, which forms the
basis for our analytical theory. In the SI, we show that the approximate
eigenstates are labeled by an oscillator quantum number $n=0,1,2,\dots$ and a
spin quantum number $\sigma=\pm 1$. These eigenstates $|n\sigma\rangle$ and
corresponding energies $E_{n\sigma}$, for $g\gg 1$, are given by

$\displaystyle|n\sigma\rangle$
$\displaystyle=\frac{1}{\sqrt{2}}\left(D^{\dagger}(g)|n,x+\rangle+\sigma
D^{\dagger}(-g)|n,x-\rangle\right)$ $\displaystyle E_{n\sigma}/\hbar$
$\displaystyle=\omega(n+\frac{\sigma}{2}e^{-2g^{2}}L_{n}(4g^{2})),$ (2)

where $D(z)\equiv\exp\left[z(a^{\dagger}-a)\right]$ is the displacement
operator, and $L_{n}$ is the Laguerre polynomial of order $n$. The state
$|n\rangle$ on the right-hand side refers to the Fock basis of the cavity,
while the states $|x\pm\rangle$ refer to the $x$-polarized spin states of the
qubit. The spectrum is plotted in Fig. 1b (adding an $g$-dependent offset
$\hbar\omega g^{2}$ for convenience). As seen in Eq. (2), the spectrum in the
DSC regime is organized into two oscillator-like ladders (one for each spin).
Moreover, for large $g$, the spectrum appears almost completely harmonic,
indicating the existence of an effective photon (or photonic quasiparticle,
which we will sometimes call a DSC photon). To understand this, we note that
for $g\gg 1$, the $\sigma_{z}$ acts as a perturbation to the remaining
Hamiltonian, $H_{\text{DSC}}/\hbar\equiv\omega
a^{\dagger}a+g\sigma_{x}(a+a^{\dagger})=\omega(b^{\dagger}b-g^{2})$, where
$b=D^{\dagger}(g\sigma_{x})aD(g\sigma_{x})=a+g\sigma_{x}$. This approximate
Hamiltonian admits a harmonic spectrum, in which the new oscillator variables
$b$ obey canonical commutation relations $[b,b^{\dagger}]=1$, and excitations
are constructed by applying further $b^{\dagger}$ operators. In other words,
the energy eigenstates of $H_{\text{DSC}}$ are Fock states of $b$, or
equivalently, displaced Fock states of $a$.

The $\sigma_{z}$ term breaks the even spacing of the ladder, leading to an
anharmonicity (equivalently, nonlinearity) which we now quantify. Without loss
of generality, we will focus on the lower-energy $\sigma=-1$ ladder, enabling
us to omit the spin index in our notation. We assess the “harmonicity” of the
spectrum by plotting successive excitation energies $E_{n+1}-E_{n}$ as a
function of $n$, as in Fig. 1c (in units of $\hbar\omega$). We will refer to
$n$ as the “photon number.” For strong and ultrastrong coupling, the spectrum
is anharmonic at the level of a single photon, leading to the familiar
phenomenon of photon blockade. For deep-strong coupling, the behavior is quite
different: the spectrum is harmonic up to some critical excitation number
($n_{c}\sim g^{2}$), and then rapidly becomes anharmonic. This may be seen
directly from the properties of $L_{n}(x)$.

To understand the relation of this strong anharmonicity to existing nonlinear
optical systems, recall that a single-mode cavity with an embedded Kerr
nonlinear medium can be described by a Hamiltonian of the form
$H_{\text{Kerr}}=\hbar\omega\left(a^{\dagger}a+\beta a^{\dagger
2}a^{2}\right)$ drummond1980quantum ; walls2007quantum , with $\beta$ a
(typically small) dimensionless coefficient. In such a system, the energy to
add an excitation is $E_{n+1}-E_{n}=\hbar\omega(1+2\beta n)$, meaning that the
deviation from harmonic behavior is linear in the intensity (proportional to
photon number). Thus, the plots of 1(c), for a photon in a Kerr medium, would
be straight lines with slope $2\beta\hbar\omega$. This linearity arises from a
low-order expansion of the nonlinear medium polarization in the cavity
electric field: in this case third order, leading to a refractive index
dependent on intensity. Here, in the case of DSC, the difference is that the
excitation energies are not linear in intensity, but instead are high-order
near the critical photon number (more or less growing exponentially before
oscillating, as in Fig. 1c), as if the effective polarization had a non-
perturbative, or infinite-order, dependence on intensity (as if the system
were described in terms of $\chi^{(n)}$s where $n\gg 1$).

##  II Light emission in the deep-strong coupling regime

Figure 2: Fock lasing due to equilibrium between high-order nonlinearity and
dissipation. (a) Light emission of DSC photons can be understood in terms of
the coupling of an emitter (e.g., a probe qubit) weakly coupled to the DSC
system, as might be realized by coupling a superconducting qubit to a flux-
qubit-LC-resonator system. The probability to stimulatedly emit DSC photons
scales as $n+1$ for small $n$, and then sharply decreases due to the sudden
anharmonicity for $n>n_{c}\sim g^{2}$. “TLS” denotes two-level system. (b)
This behavior leads to a gain medium whose gain coefficient (green lines) is
highly nonlinear. The quantum state of DSC photons will depend on how this
nonlinear gain comes into equilibrium with the loss (red lines). (c) Steady-
state intensity and power fluctuations of lasers in different coupling regimes
as a function of pump intensity. For the “harmonic” regimes (weak, and deep-
strong), a rapid growth in intensity at threshold is seen. In contrast to the
weak coupling regime (as in a “normal” laser; light blue curve), a laser
operating in the deep-strong coupling regime has its intensity saturate, and
its fluctuations vanish at high pump, converging to a high-number Fock state
(dark blue and purple curves), leading to Fock-like statistics (right). (c,
bottom) Statistics for different pump strengths for a single coupling, showing
evolution from thermality to Fock-like statistics.

While the nonlinearity perspective is largely based on the known spectrum of
DSC systems, we now use this intuition to develop the main new results of this
paper. Specifically, we study how light emission is modified by these photonic
quasiparticles. Unlike most studies of light emission with photonic
quasiparticles (reviewed for example in rivera2020light ), here we look at the
unique modifications coming from the nonlinear properties. Consider an
external qubit (denoted ‘em’, for emitter) coupled to this DSC photon. The
exact form of the coupling depends on the circuit implementation. To keep the
discussion concrete, we will consider a simple coupling Hamiltonian of the
type

$\displaystyle H$
$\displaystyle=H_{\text{Rabi}}+\frac{\omega^{\text{em}}_{0}}{2}\sigma^{\text{em}}_{z}+V$
$\displaystyle V/\hbar$
$\displaystyle=\epsilon\sigma^{\text{em}}_{x}(b+b^{\dagger})\approx\epsilon(\sigma^{\text{em}}_{+}b+b^{\dagger}\sigma^{\text{em}}_{-}),$
(3)

which couples the emitter directly to the DSC photon. Regarding the assumed
form of the Hamiltonian, we note that our conclusions are not particularly
sensitive to the exact form of the interaction 111The term
$\sigma^{\text{em}}_{x}(b+b^{\dagger})$ contains an interaction between the
dipole moment of the emitter and that of the qubit. Such interactions are to
be generically expected, as especially emphasized in recent works on
superradiant phase transitions, as well as gauge invariance in ultrastrong
coupling cavity and circuit quantum electrodynamics jaako2016ultrastrong ;
de2018breakdown ; de2018cavity ; di2019resolution ; settineri2019gauge . We
could write the term in question as $\alpha\sigma^{\text{em}}_{x}\sigma_{x}.$
Here, $\alpha=2\epsilon g$. Because this dipole-dipole term leads only to
changes in spin quantum number, and not changes in excitation number (see SI),
and because the spins are not separated by $\omega$, these terms have little
effect on the dynamics of the photon number probabilities that we consider.
For example, we find that ignoring this term altogether leads to the same
conclusions. Hence, for the purposes of the manuscript, we have taken a simple
coupling that illustrates the physics best (emission of a “b” particle by an
emitter).. What we do assume however is that $\epsilon$ is small, so that the
coupling of the external emitter to the DSC system is weak
($\epsilon\ll\omega$). Thus, the system in mind is a single resonator coupled
to two qubits, one with weak coupling and one with deep strong coupling, as
illustrated in Fig. 2a.

To understand emission and absorption of DSC photons, consider the case in
which the qubit is in its excited state $|e\rangle$ and there are $n$ DSC
photons present of spin $-1$ (e.g., occupying the state $|n,-1\rangle$ of Eq.
(2)). If the qubit is at frequency $\omega$ (same as in Eq. (2)), then the
qubit transition will be nearly resonant with the transition $n\rightarrow
n+1$ of the DSC photon, provided $n\lesssim n_{c}$. The dynamics can be
restricted to the subspace $\\{|e,n\rangle,|g,n+1\rangle\\}$, and the
probability of emission $P(n+1)$ is simply given by

$\displaystyle P(n+1)$
$\displaystyle=\frac{(n+1)\epsilon^{2}}{\Delta_{n+1}^{2}+(n+1)\epsilon^{2}}\sin^{2}\left(\sqrt{\Delta_{n+1}^{2}+(n+1)\epsilon^{2}}t\right)$
$\displaystyle\Delta_{n+1}$
$\displaystyle=\frac{\omega}{2}\left(\delta-\frac{1}{2}e^{-2g^{2}}(L_{n}(4g^{2})-L_{n+1}(4g^{2}))\right).$
(4)

Here, $\delta$ is the dimensionless detuning of the emitter and $\omega$, such
that $\omega^{\text{em}}_{0}-\omega\equiv\omega\delta$. Eq. (4) is the direct
consequence of the Jaynes-Cummings dynamics of a two-level system (emitter)
with a boson (DSC photon) with some detuning. The detuning depends on
excitation number due to the nonlinearity of the DSC photon, and the detuning
sharply rises near $n_{c}$ (Fig. 1c). In Fig. 2a, we plot the stimulated
emission probability as a function of $n$ after a small evolution time
$t\ll\epsilon^{-1}$ and for $\delta=0$. For $n<n_{c}$, $\Delta_{n+1}\approx
0$, probability is simply $(n+1)(\epsilon t)^{2}$, corresponding to stimulated
emission proportional to $n+1$, as expected for conventional photons. For
$n\gtrsim n_{c}$, the emission probability drops rapidly, because of the
corresponding rapid increase in $\Delta_{n}$. This can be understood as a type
of $N$-photon blockade, in which a system can readily accept $N$ excitations,
but not $N+1$. For $N=1$, this corresponds to the conventional photon blockade
observed and discussed extensively in strong coupling cavity QED.

###  II.1 A new type of laser

Figure 3: Uniform and tunneling states arising at the harmonic-to-anharmonic
crossover. (a) The function $(1+G_{n}))^{-1}$ which propagates the photon
probability distribution from $n$ to $n+1$, plotted as a function of coupling
strength and photon number ($G_{n}=F_{n}/\Gamma$). A harmonic-to-anharmonic
crossover occurs for a maximum photon number $n_{\text{max}}\sim g^{2}$ for
which the propagation function goes to zero. When this happens, the
probability of having photons larger than $n_{\text{max}}$ vanishes. Near
threshold (where the effective temperature of the photon goes to infinity),
this leads to nearly uniform states of the electromagnetic field sharply
cutoff at the maximum photon number (right panel). (b) When the decay rate of
the gain medium is large, the anharmonic region becomes narrower (bottom
left), and for sufficient pump intensity, the photon distribution can “tunnel”
through the barrier, evolving effectively as a coherent state. In this
tunneling regime, the distribution becomes bimodal, taking on the
characteristics of the Fock and coherent states for some pump parameters
(bottom right).

Eq. (4) displays one of the main results: that the high-order nonlinearities
arising from non-perturbative quantum electrodynamical coupling lead to a type
of gain (stimulated emission) that is correspondingly non-perturbative in
intensity. One may imagine that this type of nonlinear stimulated emission
would have implications for lasers – or in this case, masers, given that the
most imminent implementations, based on circuit QED, would be at microwave
frequencies. We will stick to the term “laser” since it has largely subsumed
masers. In this section, we show that the DSC-based gain discussed before
creates lasing into high-order Fock states (rather than coherent states).

We now show how the nonlinear gain provided by the coupling of an excited two-
level system to DSC photons can result in a laser with new steady state photon
statistics. To capture the resulting lasing dynamics in a quantum mechanical
way, we shall find an equation of motion for the reduced density matrix $\rho$
of the DSC photon (tracing out the gain medium). This equation takes into
account both the stimulated emission dynamics and the loss dynamics associated
with, for example, leakage from the cavity (which we take here for simplicity
as the primary loss mechanism for the DSC photon). In the SI, we derive the
equation using several methods, all in agreement with each other. Here, we
focus on the equation for the DSC photon occupation probabilities,
$\rho_{nn}$. Assuming that excited states of the gain medium are pumped at
rate $r$, the equation of motion for the DSC photon density matrix is found to
be:

$\dot{\rho}_{nn}=R_{n}n\rho_{n-1,n-1}-\left(R_{n+1}(n+1)\rho_{nn}+\kappa_{n}n\rho_{nn}\right)+\kappa_{n+1}(n+1)\rho_{n+1,n+1}$
(5)

Here, $R_{n}=\frac{2r\epsilon^{2}}{\Gamma^{2}+F_{n}}$ is the stimulated
emission coefficient, with $F_{n}\equiv
4n\epsilon^{2}+\frac{1}{4}\omega^{2}e^{-4g^{2}}(L_{n}(4g^{2})-L_{n-1}(4g^{2}))^{2}$.
The $R_{n}$ are plotted (green curves) in Fig. 2b for different coupling
strengths. For weak coupling, it is simply saturable gain $R(n)\sim
1/(1+n/n_{s})$ with $n_{s}$ the saturation photon number. For DSC, we see that
the gain coefficient is given by the standard saturable form for $n<n_{c}$ and
then rapidly decays for $n\gtrsim n_{c}$ (with occasional oscillations arising
from the oscillatory behavior of the Laguerre polynomials). Here,
$\kappa_{n}=\kappa|\langle n-1|a+a^{\dagger}|n\rangle|^{2}$, with $\kappa$ the
decay rate of the cavity in the absence of DSC (see SI for derivation). We
note that for simplicity, the gain medium has been taken to have population
and coherence decay rates arising from the same source (so that
$\Gamma=1/T_{1}=2/T_{2}$). This simplifies the calculations but does not
qualitatively change our conclusions.

The steady state photon probability distribution is entirely different from
that of a traditional laser, which produces a dephased coherent state. To
quantify this, we solve a recursion relation to obtain the steady-state
probability distribution
$\rho_{n,n}=Z^{-1}\prod\limits_{m=1}^{n}R_{m}/\kappa_{m}$, with $Z$ a
normalization coefficient enforcing $\sum_{n}\rho_{n,n}=1$. In Fig. 2c, we
show the intra-cavity photon number and photon fluctuations for DSC in
comparison with weak coupling. We also present the corresponding photon
statistics. In the weak coupling regime, the photon number as a function of
pump follows the canonical “S-curve” relating the input pump and output
intensity of a laser. The output intensity grows sharply for pump beyond the
threshold pump level, $r_{\text{th}}=\kappa\Gamma^{2}/2\epsilon^{2}$. The
fluctuations below threshold are essentially those of a thermal state, and far
above threshold, grow according to shot-noise (as $\sqrt{n}$, as for a
Poissonian distribution corresponding to a randomly-phased coherent state):
this is the textbook result of the laser theory of Lamb and Scully. In
contrast, the “Fock laser” ($g=5,10,18$), saturates (at $n_{c}\sim g^{2}$),
and the photon number fluctuations go to zero, leading to the quantum
statistics of a Fock state (Fig. 2c, right) as the pump increases. Fig. 2c
(bottom) further shows how the photon statistics evolve with pump and coupling
(taken for $g=10$; additional results shown in SI). Beyond threshold, the
distribution of photons (for DSC) approaches that of a thermal state of
negative temperature. Such states, as the pump is increased (and $T\rightarrow
0^{-}$), approach states where only the highest-most level is filled, with
minimal spread, which closely approximates a Fock state of $n_{c}$ DSC
photons.

To understand this Fock lasing effect, it is helpful to refer to the gain and
loss curves of Fig. 2b, as well as the steady-state distribution
$\rho_{n,n}=Z^{-1}\prod\limits_{m=1}^{n}R_{m}/\kappa_{m}$. The steady-state
distribution has the property that the probabilities are maximized where gain
meets loss, and probabilities are suppressed if one of gain or loss far
exceeds the other. In particular, the larger the angle between the gain and
loss curves (at the crossing point), the tighter the concentration of
probabilities about the mean. Increasing the pump rate $r$ further will scale
the gain curve up, leading to a steeper slope and further suppression of
photon number fluctuations, leading asymptotically to a Fock state.

Beyond these close approximations to high-photon-number Fock states, other
unusual states can arise from the equilibrium between gain and loss, due to
this sudden anharmonicity for $n\gtrsim n_{c}$. For example, near threshold,
where the number fluctuations increase dramatically, the resulting
distribution is nearly step-like, going to zero rapidly for $n_{c}$. This
anharmonicity provides a “wall” for the photon probability distribution that
is too hard to pass through, even as the fluctuations get very large near
threshold. These effects also depend on the decay rate of the gain medium: if
the decay rate is high, then it provides gain over a large bandwidth, and so
changes in the DSC photon frequency have a reduced effect on the stimulated
emission rate $R_{n}$. As a result, for increasing pump, the distribution can
“tunnel” through the wall, leading to states that interpolate between Fock and
coherent states, as well as pure coherent states for large enough pump.

##  III Discussion and outlook

Recent work on realizing deep-strong coupling of superconducting qubits to a
microwave (LC) resonator, as in yoshihara2017superconducting ;
yoshihara2018inversion , provides a path to observing the effect predicted
here. It is already possible to have control over $g$ from weak coupling to a
value of nearly 2. With a $g$ of 2, one can see from Fig. 1 that a Fock state
of three or four excitations could be pumped. For smaller $g$, in the ultra-
strong coupling regime where $0.1<g<1$, only one excitation can be created, as
a manifestation of the conventional photon blockade effect ridolfo2012photon ;
le2016fate . Thus, the behavior of our model from weak to (modest) deep strong
coupling can already be realized.

Regarding the gain medium, it is important to point out that while a typically
gain medium consisting of many emitters, the physics can also be realized by a
gain medium consisting of a single qubit. The qubit should be weakly coupled
to the same cavity as the strongly coupled qubit, and will lase, provided that
the gain from this one qubit is above threshold scully1999quantum . Single-
qubit gain is responsible for much of the exciting experiments on “one-atom
lasers” (in real mckeever2003experimental ; an1994microlaser  and artificial
atoms liu2015semiconductor ; stehlik2016double ; burkard2020superconductor ),
in which a single atom or artificial atom provides enough gain to lase.

Thus, a conceptually simpler $-$ and perhaps more attractive $-$ approach to
realize our predictions is to consider a gain medium consisting of a
continuously pumped superconducting qubit which is _weakly_ coupled to the
same resonance as the strongly coupled qubit (which for example happens if
$\epsilon\ll\kappa$). In Fig. 2, we took $\epsilon=10^{-5}\omega$, and
$\Gamma=10^{-3}\omega$. Thus, for a single gain qubit, threshold is reached
provided the quality factor of the resonator is above $5\times 10^{6}$. There
are two advances that would support reaching larger g values: the rapidly
increasing coupling constants that have been realized with superconducting
qubits (see Fig. 1 of forn2019ultrastrong ), and early estimates in this field
suggesting the possibility of $g$ values of roughly 20 devoret2007circuit .

Summarizing, we have shown a physical principle – using non-perturbative
photonic nonlinearity – which could enable lasers that produce deterministic,
macroscopic quantum states of light, such as Fock states. Part of the new
physics uncovered here, related to lasing in systems with sharply nonlinear
gain, could in principle also be extended into the optical regime. In fact, in
rivera2021macroscopic  $-$ inspired by the developments in this manuscript $-$
we discuss how trying to mimic the new “Fock lasers” predicted here, but at
_optical frequencies_. Thus, the principles established here may also give
rise to a class of new physical systems and experiments in the optical domain.

##  Appendix A Fock lasing based on deep-strong light-matter coupling

In this Supplement, we derive and extend the results of the main text.
Consider a system involving matter coupled to a cavity mode very strongly, so
that the system is in the ultra- or deep-strong coupling regime. This system
is described by the Rabi Hamiltonian of Eq. (1) of the main text (Hamiltonian
and variables re-defined here for self-containedness):

$H_{\text{Rabi}}/\hbar=\frac{1}{2}\left(\omega_{0}\sigma_{z}+\lambda\sigma_{x}\right)+\omega
a^{\dagger}a+\tilde{g}\sigma_{x}(a+a^{\dagger}),$ (6)

Here, $\omega_{0}$ is the transition frequency of the two-level system,
$\sigma_{x,z}$ are the $x$ and $z$ Pauli matrices, $\omega$ is the cavity
frequency, $a^{(\dagger)}$ is the cavity annhilation (creation) operator, and
$\tilde{g}$ is the Rabi frequency. We also non-dimensionalize the coupling as
$g=\tilde{g}/\omega$.

Let us now transfer energy into this system by means of external emitters,
treated as two-level systems of energy $\omega_{0}$. Let us assume the emitter
is primarily interacting with the cavity (as it is too far for direct
interactions with the dipole of the matter). Let us then take the full
Hamiltonian describing the coupling of one emitter to the light-matter system
as

$H/\hbar=\frac{\omega^{\text{em}}_{0}}{2}\sigma^{\text{em}}_{z}+H_{\text{Rabi}}+\epsilon\sigma_{x,\text{em}}(b+b^{\dagger}),$
(7)

which couples the emitter directly to the DSC photon. We can also consider
interactions solely between the emitter and the resonator field, replacing
$b\rightarrow a$. We consider this case as well, to show that the exact nature
of the emitter-qubit dipole-dipole coupling does not qualitatively change our
conclusions.

If the emitter is in the excited eigenstate $\ket{e}$, and it is resonant with
a transition of the Rabi Hamiltonian, the emitter can transfer energy to the
light-matter system. Upon interaction with a second emitter, if the next
transition of the Rabi model has nearly the same frequency, the system can get
further excited. A key observation is that in the deep-strong coupling regime
$g\gg\omega$, the eigenstates are approximately equally spaced, and the
excitations are oscillator-like, quite similarly to the zero-coupling case.
This should allow the possibility of reaching a very high excitation number in
the presence of many emitters, based on stimulated emission of these
oscillator modes (We will call them DSC photons). When the coupling is not
infinite, as in a realistic case, the levels are no-longer fully equally
spaced. This detuning is photon-number dependent, thus acting as a
nonlinearity which may qualitatively change the steady-state of this type of
laser.

To begin, we need to derive simple forms for the eigenstates of the Rabi
Hamiltonian in the deep-strong coupling limit. Then, we will consider their
coupling to external emitters, and write a coarse-gained equation of motion
for the density matrix of the DSC bosons, and then solve it.

###  A.1 Eigenstates of the Rabi Hamiltonian

In what follows, we will take $\omega_{0}=\omega$ (resonant) and $\lambda=0$.
In later subsections, we will analytically and numerically consider the case
of a finite $\lambda$, which is found to preserve our main findings.

In the deep-strong coupling regime, we can treat the matter term in the Rabi
Hamiltonian as a perturbation to the remainder of the Hamiltonian. The
remainder of the Hamiltonian (divided by $\hbar$), which we call
$H_{\text{DSC}}$ is

$H_{\text{DSC}}=\omega
a^{\dagger}a+\tilde{g}\sigma_{x}(a+a^{\dagger})=\omega(a^{\dagger}+g\sigma_{x})(a+g\sigma_{x})-\omega
g^{2},$ (8)

where $g\equiv\tilde{g}/\omega$ is a dimensionless measure of the coupling
strength. Introducing the displacement operator
$D(g\sigma_{x})=\exp\left[g\sigma_{x}(a^{\dagger}-a)\right]$, where we’ve
taken $g$ real without loss of generality, we have

$H_{\text{DSC}}=\omega D^{\dagger}(g\sigma_{x})a^{\dagger}aD(g\sigma_{x}),$
(9)

where we’ve omitted the overall constant $-\omega g^{2}$. From here, we can
easily see that the eigenstates of this Hamiltonian are of the form
$D^{\dagger}(\pm g)\ket{\pm x,n}$, where $\ket{x}$ denotes the x-spin basis,
and $n$ is a Fock state. In other words, the eigenstates involve the spin
being _x_ -polarized (rather than _z_ -polarization), and the photon being in
a _displaced_ Fock state (rather than just a Fock state). Clearly,

$H_{\text{DSC}}D^{\dagger}(\pm g)\ket{\pm x,n}=\omega D^{\dagger}(\pm
g)a^{\dagger}aD(\pm g)D^{\dagger}(\pm g)\ket{\pm x,n}=n\omega D^{\dagger}(\pm
g)\ket{\pm x,n}.$ (10)

Clearly then, in this limit, the eigenstates are evenly spaced, and doubly
degenerate. In fact, it can be seen as a system of two non-interacting bosons
(“DSC photons”). Introducing $b_{\sigma}=a+g\sigma_{x}$ we can write the
Hamiltonian as

$H_{\text{DSC}}=\omega b^{\dagger}b.$ (11)

It can also be easily seen that $[b,b^{\dagger}]=1$.

The degeneracy of the DSC Hamiltonian is split by the matter Hamiltonian. We
can find the resulting eigenstates and eigenenergies using degenerate first-
order perturbation theory. The “good” eigenbasis of the problem is

$\ket{n,\sigma}=\frac{1}{\sqrt{2}}\left(D^{\dagger}\ket{+x,n}+\sigma
D\ket{-x,n}\right),$ (12)

where $\sigma=\pm 1$, and a displacement operator without an argument implies
that the argument is $g$. The energies of the resulting states are

$\displaystyle E_{n\sigma}$
$\displaystyle=\frac{\omega}{2}\bra{n,\sigma}\sigma_{z}\ket{n,\sigma}$
$\displaystyle=\frac{\omega}{4}\left(\bra{+x,n}D+\sigma\bra{-x,n}D^{\dagger}\right)\sigma_{z}\left(D^{\dagger}\ket{+x,n}+\sigma
D\ket{-x,n}\right)$
$\displaystyle=\sigma\frac{\omega}{4}\left(\bra{n}D^{2}\ket{n}+\bra{n}D^{\dagger
2}\ket{n}\right)$
$\displaystyle=\sigma\frac{\omega}{2}\bra{n}D^{2}\ket{n}\equiv\sigma\frac{\omega}{2}D_{n},$
(13)

where $D_{n}=\bra{n}D^{2}(g)\ket{n}=\bra{n}D(2g)\ket{n}$. These eigenstates
and energies are sufficiently accurate, even for $g=2$ or $g=3$.

####  A.1.1 Evaluation of $D_{n}$.

Let us evaluate the $D_{n}$. To do so, we write:

$\displaystyle D_{n}(2z)$
$\displaystyle=\bra{n}DD\ket{n}=\frac{1}{n!}\bra{0}a^{n}DDa^{\dagger
n}\ket{0}=\frac{1}{n!}\bra{0}D(D^{\dagger}a^{n}D)(Da^{\dagger
n}D^{\dagger})D\ket{0}$
$\displaystyle=\frac{1}{n!}\bra{-z}(a+z)^{n}(a^{\dagger}-z^{*})^{n}\ket{z}.$
(14)

To proceed, insert a “complete” set of states using the over-completeness of
the coherent states. That leaves us with

$D_{n}(2z)=\frac{1}{\pi n!}\int
d^{2}\alpha\bra{-z}(a+z)^{n}\ket{\alpha}\bra{\alpha}(a^{\dagger}-z^{*})^{n}\ket{z}=\frac{1}{\pi
n!}\int
d^{2}\alpha(\alpha+z)^{n}(\alpha^{*}-z^{*})^{n}\langle\alpha|z\rangle\langle-z|\alpha\rangle.$
(15)

Using the rule for the overlap of two coherent states, we have

$D_{n}(2z)=\frac{1}{\pi n!}\int d\alpha
d\alpha^{*}(\alpha+z)^{n}(\alpha^{*}-z^{*})^{n}e^{-\alpha\alpha^{*}-zz^{*}+\alpha^{*}z-z^{*}\alpha},$
(16)

where we have written things this way to emphasize that $\alpha$ and
$\alpha^{*}$ are independent variables. We can now write this as

$D_{n}(2z)=\frac{e^{zz^{*}}}{\pi n!}\int d\alpha
d\alpha^{*}(\alpha+z)^{n}(\alpha^{*}-z^{*})^{n}e^{-\alpha\alpha^{*}+z(\alpha^{*}-z^{*})-z^{*}(\alpha+z)}.$
(17)

Then, we transform variables as $\alpha\rightarrow\alpha-z$ and
$\alpha^{*}\rightarrow\alpha^{*}+z^{*}$ to get

$D_{n}(2z)=\frac{e^{zz^{*}}}{\pi n!}\int d\alpha
d\alpha^{*}~{}\alpha^{n}\alpha^{*n}e^{-(\alpha-z)(\alpha^{*}+z^{*})+z\alpha^{*}-z^{*}\alpha}=\frac{e^{2zz^{*}}}{\pi
n!}\int d\alpha
d\alpha^{*}~{}\alpha^{n}\alpha^{*n}e^{-\alpha\alpha^{*}+2z\alpha^{*}-2z^{*}\alpha}.$
(18)

This can be generated from simpler integrals by differentiation, as:

$D_{n}(2z)=(-1)^{n}\frac{e^{2zz^{*}}}{\pi
n!}\frac{\partial^{2n}}{\partial^{n}(2z)\partial^{n}(2z^{*})}\int d\alpha
d\alpha^{*}~{}e^{-\alpha\alpha^{*}+2z\alpha^{*}-2z^{*}\alpha}.$ (19)

Completing the square in the remaining integral gives

$\displaystyle D_{n}(2z)$ $\displaystyle=(-1)^{n}\frac{e^{2zz^{*}}}{\pi
n!}\frac{\partial^{2n}}{\partial^{n}(2z)\partial^{n}(2z^{*})}e^{-4zz^{*}}\int
d\alpha d\alpha^{*}~{}e^{-\alpha\alpha^{*}+2z\alpha^{*}-2z^{*}\alpha+4zz^{*}}$
$\displaystyle=(-1)^{n}\frac{e^{2zz^{*}}}{\pi
n!}\frac{\partial^{2n}}{\partial^{n}(2z)\partial^{n}(2z^{*})}e^{-4zz^{*}}\int
d\alpha d\alpha^{*}~{}e^{-(\alpha-2z)(\alpha^{*}+2z^{*})}.$ (20)

Shifting variables as $\alpha\rightarrow\alpha+2z$ and
$\alpha\rightarrow\alpha-2z^{*}$, and performing the final Gaussian integral,
we have

$D_{n}=(-1)^{n}\frac{e^{2zz^{*}}}{n!}\frac{\partial^{2n}}{\partial(2z)^{n}\partial(2z^{*})^{n}}e^{-4zz^{*}}.$
(21)

These are related to Laguerre polynomials. To see this, take the derivative
with respect to $z^{*}$. We will also use the notation $z\rightarrow x/2$ and
$z^{*}\rightarrow y/2$ for clarity.

$D_{n}(2z)=(-1)^{n}\frac{e^{xy/2}}{n!}\frac{\partial^{2n}}{\partial
x^{n}\partial
y^{n}}e^{-xy}=\frac{e^{xy/2}}{n!}\frac{\partial^{n}}{\partial(xy)^{n}}(xy)^{n}e^{-(xy)}=e^{-xy/2}\left(\frac{e^{xy}}{n!}\frac{\partial^{n}}{\partial(xy)^{n}}(xy)^{n}e^{-(xy)}\right).$
(22)

From the Rodrigues formula for the Laguerre polynomials, we then have

$D_{n}(2z)=e^{-2|z|^{2}}L_{n}(4|z|^{2}).$ (23)

Given that the $z=g$, we have then that the the level splitting is given by

$E_{n\sigma}=n\omega+\sigma\frac{\omega}{2}e^{-2|g|^{2}}L_{n}(4|g|^{2})$ (24)

###  A.2 Time-evolution of the coupled system

With the approximate eigenstates of the Rabi Hamiltonian, we now want to
understand the full dynamics of $H$. We will take advantage of the fact that
for a laser, $\epsilon$ is small, and in particular, $\epsilon\ll\omega$, so
that the rotating-wave approximation is valid. In this system, the rotating
wave approximation consists of only considering the dynamics within degenerate
subspaces of the unperturbed Hamiltonian

$H_{0}=\frac{\omega_{0}}{2}\sigma_{z,\text{em}}+H_{\text{Rabi}}.$ (25)

The eigenstates of the problem are $\ket{k}\ket{n\sigma}$, where now $k=0,1$
denotes emitter states (ground is zero, excited is one). The energies of such
states are (up to a shift)

$E_{kn\sigma}=(n+k(1+\delta))\omega+\sigma\frac{\omega}{2}D_{n}.$ (26)

where we have taken $\omega_{0}=(1+\delta)\omega$. From here on out, let us
assume $\delta\ll\omega$. In that case, it is easy to see that the following
four states form our nearly degenerate subspace:

$\\{\ket{1,n-1,+},\ket{0,n,+},\ket{1,n-1,-},\ket{0,n,-}\\}.$ (27)

We now need to understand the action of the interaction Hamiltonian
$V\equiv\epsilon\sigma_{x,\text{em}}(b+b^{\dagger})$ on this subspace. First
of all,

$\bra{k^{\prime}n^{\prime}\sigma^{\prime}}V\ket{kn\sigma}=0~{}\text{if}~{}k=k^{\prime}.$
(28)

For $k=-k^{\prime}$, we have

$\displaystyle\bra{-kn^{\prime}\sigma^{\prime}}V\ket{kn\sigma}$
$\displaystyle=\frac{\epsilon}{2}\left(\bra{+x,n^{\prime}}D+\sigma^{\prime}\bra{-x,n^{\prime}}D^{\dagger}\right)D^{\dagger}(g\sigma_{x})(a+a^{\dagger})D(g\sigma_{x})\left(D^{\dagger}\ket{+x,n}+\sigma
D\ket{-x,n}\right)$
$\displaystyle=\epsilon\frac{1+\sigma\sigma^{\prime}}{2}(\sqrt{n}\delta_{n^{\prime},n-1}+\sqrt{n+1}\delta_{n^{\prime},n+1}).$
(29)

From these matrix elements, we see that: if the pseudo-spin ($\sigma$) is
conserved, then a non-zero matrix element occurs only when the boson number
changes by 1. Noting that $b=a+g\sigma_{x}$, we can also readily describe
interactions using $a$ as $\langle
n^{\prime}\sigma^{\prime}|a+a^{\dagger}|n\sigma\rangle=\sqrt{n}\delta_{n^{\prime},n-1}+\sqrt{n+1}\delta_{n^{\prime},n+1}-2g\frac{1-\sigma\sigma^{\prime}}{2}\delta_{nn^{\prime}}$,
such that: when the pseudo-spin changes, non-zero matrix elements occur only
when the boson number is conserved. When the spin is conserved, the matrix
elements are the same as for $b+b^{\dagger}$. Since only states with different
photon number differ appreciably in frequency (and in particular, will be
resonant with the emitter we introduce), the interactions are effectively the
same whether we describe $a$ or $b$. This is also to say that any modification
in the coefficient of the dipole-dipole interaction between emitter and qubit
will lead to the same result insofar as DSC photon dynamics are concerned. In
Fig. S1, we show the matrix elements of $a$ and $b$ between adjacent states of
the same spin, as well as $a^{\dagger}a$ and $b^{\dagger}b$. For $n<n_{c}\sim
g^{2}$, they behave as one might expect for an oscillator.

We should note that beyond $n_{c}$, these states and matrix elements that we
calculate based on degenerate perturbation theory are expected to change
significantly. However, the approximate result turns out to describe the
system well because the probabilities to find photon numbers beyond $n_{c}$
are strongly suppressed in the Fock laser, rendering the description
relatively insensitive to these details.

Figure 4: Matrix elements of $a^{\dagger}a$, $b^{\dagger}b$ and $a,b$, showing
that $b$ counts excitations of the DSC system over the full range of
eigenstates. However, up to $n_{c}$, $a$ and $b$ act similarly.

Based on these considerations, we see that the Hamiltonian in the degenerate
subspace may be written as:

$\displaystyle H$
$\displaystyle=(\delta\omega+\frac{\omega}{2}D_{n-1})\ket{1,n-1,+}\bra{1,n-1,+}+(\frac{\omega}{2}D_{n})\ket{0,n,+}\bra{1,n,+}$
$\displaystyle+\epsilon\sqrt{n}\ket{1,n-1,+}\bra{0,n,+}+\text{h.c}$
$\displaystyle+(\delta\omega-\frac{\omega}{2}D_{n-1})\ket{1,n-1,-}\bra{1,n-1,-}+(-\frac{\omega}{2}D_{n})\ket{0,n,-}\bra{0,n,-}$
$\displaystyle+\epsilon\sqrt{n}\ket{1,n-1,-}\bra{0,n,-}+\text{h.c}.$ (30)

As we can see here, there are two independent blocks of the Hamiltonian (for
each pseudo-spin) and we can thus study them separately. Let us assume that
we’re at zero temperature, and so the ground state has the $-$ pseudo-spin,
which we assume to be conserved for all times. In that case, we can work with
the simple 2x2 Hamiltonian

$H_{\text{eff}}=\omega(\delta-\frac{1}{2}D_{n-1})\ket{1,n-1}\bra{1,n-1}-\frac{1}{2}\omega
D_{n}\ket{0,n}\bra{0,n}+\epsilon\sqrt{n}\ket{1,n-1}\bra{0,n}+~{}\text{h.c},$
(31)

where the pseudo-spin label has been dropped. This can be written in terms of
Pauli matrices as

$H_{n}=\frac{\omega}{2}(\delta-\frac{1}{2}(D_{n-1}+D_{n}))I+\frac{\omega}{2}(\delta-\frac{1}{2}(D_{n-1}-D_{n}))\sigma_{z}+\epsilon\sqrt{n}\sigma_{x}.$
(32)

Introducing $E_{n}=\frac{\omega}{2}(\delta-\frac{1}{2}(D_{n-1}+D_{n}))$,
$\Delta_{n}=\frac{\omega}{2}(\delta-\frac{1}{2}(D_{n-1}-D_{n}))$, we have very
simply

$H_{n}=E_{n}I+\Delta_{n}\sigma_{z}+\epsilon\sqrt{n}\sigma_{x}.$ (33)

As we will see in the next section, we need to know how states of the form
$\ket{1,n-1}$ evolve over time. Thus we need

$e^{-it(\Delta_{n}\sigma_{z}+\epsilon\sqrt{n}\sigma_{x})}\equiv
e^{-i|U_{n}|t(\hat{U}_{n}\cdot\sigma)},$ (34)

with
$\hat{U}_{n}=\frac{(\Delta_{n},0,\epsilon\sqrt{n})}{\sqrt{\Delta_{n}^{2}+n\epsilon^{2}}}$
and $|U_{n}|=\sqrt{\Delta_{n}^{2}+n\epsilon^{2}}$. Since
$\hat{U}\cdot\sigma=1$, we have

$e^{-i|U_{n}|t(\hat{U}_{n}\cdot\sigma)}=\cos(|U_{n}|t)-i(\hat{U}_{n}\cdot\sigma)\sin(|U_{n}|t).$
(35)

Therefore

$e^{-iHt/\hbar}\ket{1,n-1}=(\cos(|U_{n}|t)-i\hat{U}_{nz}\sin(|U_{n}|t))\ket{1,n-1}-i\hat{U}_{nx}\sin(|U_{n}|t)\ket{0,n},$
(36)

So, the probability of remaining in the same state is

$P(1,n-1)=\cos^{2}(|U_{n}|t)+\frac{\Delta_{n}^{2}}{\Delta_{n}^{2}+n\epsilon^{2}}\sin^{2}(|U_{n}|t)=1-\frac{n\epsilon^{2}}{\Delta_{n}^{2}+n\epsilon^{2}}\sin^{2}(|U_{n}|t),$
(37)

while the probability of transitioning is

$P(0,n)=\frac{n\epsilon^{2}}{\Delta_{n}^{2}+n\epsilon^{2}}\sin^{2}(|U_{n}|t).$
(38)

##  Appendix B Equation of motion for DSC photons

Now we consider the description of laser action. To do so, we formulate an
equation for how the density matrix of the DSC photon changes due to
stimulated emission by the emitter. The method of analysis presented closely
follows the coarse-grained density matrix technique used to describe
conventional lasers. It is described in many books, such as scully1999quantum
; walls2007quantum . Suppose we have our emitter coupled to the light-matter
(DSC) system. The emitter unit starts in the state $|i\rangle$ and the DSC
system is taken to have a density matrix $\rho_{\text{DSC}}$, so that the
initial density matrix of the total system $\rho_{\text{tot}}$ is given by
$\rho_{\text{tot}}(t)=|i\rangle\langle i|\rho_{\text{DSC}}(t)$. Let us look
for an equation describing only the evolution of the DSC system. Assuming the
interaction over a time $T$ corresponds to the evolution operator $U$, we have
that

$\rho_{\text{tot}}(t+T)=U(T)|i\rangle\langle
i|\rho_{\text{DSC}}(t)U^{\dagger}(T).$ (39)

Let us express all operators in terms of their matrix elements, writing the
above equation as

$\rho_{\text{tot}}(t+T)=\sum\limits_{ff^{\prime},mm^{\prime},nn^{\prime}}\rho_{\text{DSC},nn^{\prime}}(t)\langle
fm|U(T)|in\rangle\langle
in^{\prime}|U^{\dagger}(T)|fm^{\prime}\rangle|fm\rangle\langle
f^{\prime}m^{\prime}|,$ (40)

where $\rho_{\text{DSC},nn^{\prime}}$ denotes the matrix elements of the DSC
system, for a fixed pseudo-spin. More compactly,

$\rho_{\text{tot}}(t+T)=\sum\limits_{ff^{\prime},mm^{\prime},nn^{\prime}}\rho_{\text{DSC},nn^{\prime}}(t)U_{fm,in}(T)U^{*}_{f^{\prime}m^{\prime},in^{\prime}}(T)|fm\rangle\langle
f^{\prime}m^{\prime}|,$ (41)

The DSC photon density matrix,
$\rho_{\text{DSC}}=\text{tr}_{\text{em}}\rho_{\text{tot}}$, can then be
expressed as

$\rho_{\text{DSC}}(t+T)=\sum\limits_{k}\langle
k|\rho_{\text{DSC}}(t+T)|k\rangle=\sum\limits_{k,mm^{\prime},nn^{\prime}}\rho_{\text{DSC},nn^{\prime}}(t)U_{km,in}(T)U^{*}_{km^{\prime},in^{\prime}}(T)|m\rangle\langle
m^{\prime}|.$ (42)

From here, a number of approaches can be followed. If there is no loss in the
system, then the density matrix of the total system upon the next iteration is
simply $\rho_{\text{DSC}}(t+T)=|i\rangle\langle i|\rho_{\text{DSC}}(t+T)$ and
this procedure can be iterated in a discrete fashion. The evolution can also
be seen as continuous if, over time $T$, the change in the density matrix is
small. This doesn’t describe the early stages of the evolution, but it can
describe later stages once there are many bosons in the system. If there is a
steady state, then the continuous evolution must describe the run-up to the
steady state, as changes get smaller over time. In such a case, we have

$\dot{\rho}_{mm^{\prime}}=r\left(\sum\limits_{k,nn^{\prime}}U_{km,in}(T)U^{*}_{km^{\prime},in^{\prime}}(T)-\delta_{n,m}\delta_{n^{\prime},m^{\prime}}\right)\rho_{nn^{\prime}},$
(43)

where $r=N/T$ is the number of excited emitters introduced into the system in
time $T$. We have also dropped the “DSC” subscript for the DSC photon for
brevity. These terms in the evolution of the density matrix describe the gain
in the system. In addition, since there are losses associated with the cavity,
the emitter, and the matter coupled to the cavity, we need to describe those.
For simplicity, we will assume the emitter has loss, and so does the DSC
photon, but not the matter (qualitatively similar results arise if the matter
has loss).

###  B.1 Lindblad terms

Here we describe the effect of dissipation of the DSC photon on the equation
of motion for its density matrix. Let’s assume for simplicity that the cavity
loss the primary source of dissipation in the problem. For weak coupling, the
standard prescription is to add a Lindblad term to the Liouvillian which
prescribes the evolution of the density matrix. The Lindblad term would be (at
zero temperature)
$\mathcal{D}[a]\rho\equiv-\frac{\kappa}{2}(a^{\dagger}a\rho+\rho
a^{\dagger}a-2a\rho a^{\dagger})$. As is well known from studies of
dissipation in ultra-strong coupling of light and matter, the use of the
standard Lindblad term leads to unphysical excitations (in the energy
eigenbasis), even at zero temperature, and zero pumping
beaudoin2011dissipation . Part of the issue is that in the USC regime, the $a$
operator can create excitations in the eigenbasis, clearly not representing
dissipation. Framed in terms of the standard derivation, the issue could be
said that the interaction picture $a$ operator has negative frequencies, and
the use of white noise (with frequencies $-\infty$ to $\infty$) introduces
contributions from these negative frequencies ridolfo2012photon ;
kockum2019ultrastrong . The issue can be rectified by keeping in mind the
positive-frequency nature of the reservoir.

We now use this procedure to describe dissipation in the deep-strong coupling
regime. Although the technique has been worked out for ultra-strong coupling,
there is a commonly used assumption in the final result that all transitions
have different frequencies, which does not necessarily hold in DSC, when the
energy ladder is quasi-harmonic. Interestingly, as we will show from a
physical dissipator, the issues described above create much less error in the
DSC regime, and the use of an operator like $a$ or $b$ produces a similar
result to a proper positive-frequency jump operator, as their negative
frequency parts get exponentially suppressed.

Let us consider the Lindblad term arising from a system-bath coupling of the
form

$V=J\sum\limits_{k}(V_{k}b_{k}+V^{*}_{k}b_{k}^{\dagger}),$ (44)

where $J$ is a DSC system operator (e.g., $a+a^{\dagger}$ or $b+b^{\dagger}$),
and the $b_{k}$ are the bath operators, satisfying
$[b_{k},b_{k^{\prime}}^{\dagger}]=\delta_{kk^{\prime}}$. The couplings $V_{k}$
between system and bath are weak. To isolate the positive-frequency parts of
$J$, we express it in its energy eigenbasis as
$J=\sum\limits_{n>m}J_{mn}T_{mn}+\sum\limits_{m>n}J_{mn}T_{mn}+\sum\limits_{n}J_{nn}T_{nn}\equiv
J^{(+)}+J^{(-)}+J^{0}$ , with $J_{mn}=\langle m|J|n\rangle$ and
$T_{mn}=|m\rangle\langle n|$.

In what follows, we will consider the bath to be concentrated around $\omega$,
but broadband enough that the white-noise approximation may be made for any
transitions we consider. For example, a bath with a half-bandwidth of 10% of
$\omega$ would be sufficient for the values of $g,V_{k}$ we consider. It would
include all active transitions of the form $n\rightarrow n+1$, but would not
include higher transitions (though the matrix elements for them are small
anyway), and in the presence of a $\lambda$ term, it would also not include
transitions that only change spin (for $\lambda=0$, the two spins are very
nearly degenerate and so the argument should be treated with more care).
Therefore, we may describe the interaction of Eq. (39) within the rotating
wave approximation, instead considering

$V\approx\sum\limits_{k}(V_{k}b_{k}J^{(+)}+V^{*}_{k}b_{k}^{\dagger}J^{(-)}).$
(45)

We note that it is not necessary to take the RWA at this stage, but it makes
the subsequent manipulations simpler.

Thus, we may approximate the evolution of the reduced density matrix of the
DSC system (in the interaction picture) to second-order in time-dependent
perturbation theory, as:

$\dot{\rho}_{\text{DSC},I}=-i\text{tr}_{b}\left([V_{I}(t),\rho(0)]\right)-\int\limits_{0}^{t}dt^{\prime}~{}\text{tr}_{b}\left(\left[V_{I}(t),\left[V_{I}(t^{\prime}),\rho_{I}(t^{\prime})\right]\right]\right),$
(46)

where $\rho_{I}$ is the system-bath density matrix, $\rho_{\text{DSC},I}$ is
the system density matrix, $V_{I}$ is the system-bath coupling in the
interaction picture, and $\text{tr}_{b}$ denotes the partial trace with
respect to the bath. For simplicity, we will consider the bath at zero
temperature. Upon taking the trace with respect to the bath, the term which is
linear in $V_{I}$ will vanish, and the equation of motion becomes

$\dot{\rho}_{\text{DSC},I}=-\int\limits_{0}^{t}dt^{\prime}~{}\text{tr}_{b}\left(V_{I}(t)V_{I}(t^{\prime})\rho_{I}(t^{\prime})+\rho_{I}(t^{\prime})V_{I}(t^{\prime})V_{I}(t)-V_{I}(t)\rho_{I}(t^{\prime})V_{I}(t^{\prime})-V_{I}(t^{\prime})\rho_{I}(t^{\prime})V_{I}(t)\right).$
(47)

The first term may be simplified, taking the trace with respect to the bath
variables, as

$-\int\limits_{0}^{t}dt^{\prime}\int\limits_{0}^{\infty}d\omega~{}D(\omega)|V(\omega)|^{2}e^{i\omega(t^{\prime}-t)}J^{(-)}_{I}(t)J^{(+)}_{I}(t^{\prime})\rho_{\text{DSC}}(t^{\prime}),$
(48)

where $D(\omega)$ is the density of bath states, and we have replaced the sum
over $k$ by an integral over bath frequencies. Since an operator of the form
$J^{(+)}$ is a pure de-excitation operator, no spurious excitations are
introduced, and the integration limits may be extended to $-\infty$. Doing so,
and making the white noise approximation, one immediately finds that the term
evaluates to
$-\frac{\kappa}{2}J^{(-)}_{I}(t)J^{(+)}_{I}(t)\rho_{\text{DSC}}(t)$, where
$\kappa=2\pi\rho|V|^{2}$. A similar manipulation for the remaining terms
yields that the free dissipation dynamics of the DSC Hamiltonian are governed
by

$\dot{\rho}_{\text{DSC},I}=-\frac{\kappa}{2}\left(J_{I}^{(-)}J_{I}^{(+)}\rho_{\text{DSC}}+J_{I}^{(-)}J_{I}^{(+)}\rho_{\text{DSC}}-2J_{I}^{(+)}\rho_{\text{DSC}}J_{I}^{(-)}\right).$
(49)

Let us use this to find the contribution of dissipation to the equation of
motion for the populations, $\rho_{nn}$. From here on out, we will suppress
the “DSC” subscript. We will ignore the spin degree of freedom (and restrict
the dynamics to a single spin ladder). Although this is not rigorous, one
expects this to capture well the dynamics of the DSC photon number as, for
$\lambda=0$, one will just expect the nearly degenerate spins to be mixed,
with little change of the oscillator quantum numbers. We validate this
numerically. For finite $\lambda$ the spin ladders can be split appreciably,
and so they will decouple. Consider a $J$ of the form $b+b^{\dagger}$. As
discussed in the main text, $b$ is a pure de-excitation operator, and
$b^{\dagger}$ is a pure creation operator. Therefore, $J^{(+)}=b$. Using the
fact that $\langle
n^{\prime}\sigma|b|n\sigma\rangle=\sqrt{n}\delta_{n^{\prime},n-1}$ and
$b^{\dagger}b|n\sigma\rangle=n|n\sigma\rangle$, one immediately arrives at

$\dot{\rho}_{nn}=-\kappa n\rho_{nn}+\kappa(n+1)\rho_{n+1,n+1},$ (50)

which is similar to the form one would expect for damping of a conventional
photon. This is perhaps unsurprising in light of the fact that the DSC photon
is essentially harmonic up to $n_{c}\sim g^{2}$. It is worth noting that the
matrix elements derived for $a,b$ in Eq. (24) are based on first-order
degenerate perturbation theory. Beyond $n_{c}$, these approximations do not
hold up and the states and matrix elements change significantly. However, the
approximate result turns out to describe the system well because the
probabilities to find photon numbers beyond $n_{c}$ are strongly suppressed.
It is also worth noting that if we chose $a$ instead of $b$ as the jump
operator, when we neglect spin, the matrix elements are the same. Numerically,
we find that whether we choose $a$ or $b$ as the jump operator, negligible
levels of excitations are created in the ground state, and the steady-state of
the Fock laser we describe is not qualitatively changed. These numerics are
shown in the last section.

####  B.1.1 Rate equations

We will now obtain a closed set of equations for the diagonals of the DSC
density matrix, to get the probability of different Fock state occupations of
the DSC photons. Setting $m=m^{\prime}$, we have

$\dot{\rho}_{mm}=r\left(\sum\limits_{k,nn^{\prime}}U_{km,in}(T)U^{*}_{km,in^{\prime}}(T)-\delta_{n,m}\delta_{n^{\prime},m}\right)\rho_{nn^{\prime}}.$
(51)

The set of equations for the coarse grained density matrix is only closed when
$U_{km,in}(T)U^{*}_{km,in^{\prime}}(T)$ is zero unless $n=n^{\prime}$. In that
case, we have

$\dot{\rho}_{mm}=r\left(\sum\limits_{k,n}|U_{km,in}(T)|^{2}-\delta_{n,m}\delta_{n^{\prime},m}\right)\rho_{nn^{\prime}}=r\left(\sum\limits_{k,n}|U_{km,in}(T)|^{2}\rho_{nn}-\rho_{mm}\right).$
(52)

We can now note the conditions under which the equations for the populations
become closed. We require $U_{km,in}(T)U^{*}_{km,in^{\prime}}(T)$ is zero
unless $n=n^{\prime}$. This is equivalent to saying that a transition
$in\rightarrow km$ and $in^{\prime}\rightarrow km$ are not simultaneously
possible. Supposing $i$ is also an eigenstate of the light-matter system, and
that we are in the RWA, this statement appears to amount to energy
conservation, as transitions are assumed to be only efficient if they are
resonant, so that $E_{i}+E_{n}=E_{k}+E_{m}$. Therefore the condition that
$U_{km,in}(T)U^{*}_{km,in^{\prime}}(T)\neq 0$ for $n\neq n^{\prime}$ requires
$E_{n}=E_{n^{\prime}}$, which, for a single oscillator, requires
$n=n^{\prime}$.

In the weak coupling regime then, we have (adding in the photon losses)

$\dot{\rho}_{mm}=r\sum\limits_{k,n}|U_{km,in}(T)|^{2}\rho_{nn}-r\rho_{mm}+\kappa(m+1)\rho_{m+1,m+1}-\kappa
m\rho_{mm}.$ (53)

Let’s now consider the case of the emitter coupled to our light-matter system.
Since we inject emitters in the excited state, we have $i=1$. The state $1n$
couples only to $1n$ and $0(n+1)$. So, the sum over probabilities leaves only
the scattering matrix coefficients $U_{1m,1m}$ and $U_{0m,1(m-1)}$. Therefore,
the coarse-grained equation simplifies to:

$\dot{\rho}_{mm}=r(|U_{0m,1(m-1)}(T)|^{2}\rho_{m-1,m-1}+|U_{1m,1m}(T)|^{2}\rho_{mm})-r\rho_{mm}+\kappa(m+1)\rho_{m+1,m+1}-\kappa
m\rho_{mm}.$ (54)

We found these probability coefficients when studying the dynamics of the
Hamiltonian in the degenerate subspace of fixed pseudo-spin. Plugging in the
results there, we have 222To check a limiting case, we set the detunings are
zero. In that case, $|U_{n}|=\sqrt{n}\epsilon$, and we have
$\dot{\rho}_{nn}=r\sin^{2}(\epsilon
T\sqrt{n})\rho_{n-1,n-1}-\left(r\sin^{2}(\epsilon T\sqrt{n+1})+\kappa
n\right)\rho_{nn}+\kappa(n+1)\rho_{n+1,n+1}$ This coincides exactly with the
equation of motion of the so-called micromaser, which describes the
interaction of injected two-level atoms interacting with a cavity (in the
perturbative coupling regime $g\ll\omega$). This is quite interesting as the
micromaser equations assume $g\ll 1$, while here, we are starting from the
limit $g\gg 1$. Moreover, by averaging over decay times as we do in the next
subsection, we will find exactly the standard Scully-Lamb master equation for
a conventional laser. What’s happening here is that in the weak-coupling
regime, assuming the emitter is resonant, the detunings also approximately
vanish between the nearly degenerate levels. And so we get a similar equation,
except that it the conventional case, it is in the photon basis, and here it
is in the DSC photon basis.

$\dot{\rho}_{nn}=\frac{rn\epsilon^{2}}{\Delta_{n}^{2}+n\epsilon^{2}}\sin^{2}(|U_{n}|T)\rho_{n-1,n-1}-\left(\frac{r(n+1)\epsilon^{2}}{\Delta_{n+1}^{2}+(n+1)\epsilon^{2}}\sin^{2}(|U_{n+1}|T)+\kappa
n\right)\rho_{nn}+\kappa(n+1)\rho_{n+1,n+1}.$ (55)

To proceed, we must the emitter loss ($T_{1}$ and $T_{2}$ decay) into account.
Assuming that the emitter loss manifests as exponential decay with rate
$\Gamma$, the effect is to average the probability coefficients over $T$ with
probability distribution $P(T)=\Gamma e^{-\Gamma T}$. Noting that

$\Gamma\int dT~{}e^{-\Gamma T}\sin^{2}(\alpha
T)=\frac{2\alpha^{2}}{\Gamma^{2}+4\alpha^{2}},$ (56)

we have

$\displaystyle\dot{\rho}_{nn}$
$\displaystyle=r\frac{n\epsilon^{2}}{\Delta_{n}^{2}+n\epsilon^{2}}\frac{2U_{n}^{2}}{\Gamma^{2}+4U_{n}^{2}}\rho_{n-1,n-1}$
$\displaystyle-\left(r\frac{(n+1)\epsilon^{2}}{\Delta_{n+1}^{2}+(n+1)\epsilon^{2}}\frac{2U_{n+1}^{2}}{\Gamma^{2}+4U_{n+1}^{2}}+\kappa
n\right)\rho_{nn}$ $\displaystyle+\kappa(n+1)\rho_{n+1,n+1}.$ (57)

Noting that $U_{n}^{2}=\Delta_{n}^{2}+n\epsilon^{2}$, we have

$\displaystyle\dot{\rho}_{nn}$
$\displaystyle=r\frac{2n\epsilon^{2}}{\Gamma^{2}+4(\Delta_{n}^{2}+n\epsilon^{2})}\rho_{n-1,n-1}$
$\displaystyle-\left(r\frac{2(n+1)\epsilon^{2}}{\Gamma^{2}+4(\Delta_{n+1}^{2}+(n+1)\epsilon^{2})}+\kappa
n\right)\rho_{nn}$ $\displaystyle+\kappa(n+1)\rho_{n+1,n+1}.$ (58)

Assuming resonance between the emitter and the light-matter system, we have
finally

$\dot{\rho}_{nn}=\frac{2rn\epsilon^{2}}{\Gamma^{2}+F(n)}\rho_{n-1,n-1}-\left(\frac{2r(n+1)\epsilon^{2}}{\Gamma^{2}+F(n+1)}+\kappa
n\right)\rho_{nn}+\kappa(n+1)\rho_{n+1,n+1},$ (59)

with the nonlinearity, $F(n)$ defined as

$F(n)=4n\epsilon^{2}+\frac{1}{4}\omega^{2}e^{-4g^{2}}(L_{n}(4g^{2})-L_{n-1}(4g^{2}))^{2}.$
(60)

Here, we have used
$\Delta_{n}=\frac{\omega}{2}(\delta-\frac{e^{-2g^{2}}}{2}(L_{n-1}(4g^{2})-L_{n}(4g^{2})))$
with $\delta=0$.

###  B.2 Steady-state dynamics

Perhaps one of the most important results is the steady-state dynamics of the
system. Thus we want to solve $\dot{\rho}_{n}=0$ with the constraint
$\sum\limits_{n}\rho_{n}=1$ (introducing the shorthand $\rho_{n}=\rho_{nn}$).
Writing the steady-state equation as
$0=A_{n}\rho_{n-1}+B_{n}\rho_{n}+C_{n}\rho_{n+1}$, we have the recursion
relation: $\rho_{n+1}=-\frac{B_{n}\rho_{n}+A_{n}\rho_{n-1}}{C_{n}},$ with
$\rho_{0}=1$ and $\rho_{-1}=0$. Since any scale multiple of $\rho$ also solves
this equation, we can normalize the solution at the end to satisfy the
normalization constraint.

This equation can be simplified by noting that $B_{n}=-(A_{n+1}+C_{n-1})$. We
thus have $A_{n}\rho_{n-1}-A_{n+1}\rho_{n}-C_{n-1}\rho_{n}+C_{n}\rho_{n+1}=0$
or alternatively

$A_{n}\rho_{n-1}-C_{n-1}\rho_{n}=A_{n+1}\rho_{n}-C_{n}\rho_{n+1}.$ (61)

Defining the difference $S_{n}=A_{n}\rho_{n-1}-C_{n-1}\rho_{n}$, we see that
$S_{n}=S_{n+1}$. Since $S_{0}=A_{0}\rho_{-1}-C_{-1}\rho_{0}=0$, we have that
$S_{n}=0$ for all $n$, and thus the simpler recursion relation:

$\rho_{n+1}=\frac{A_{n+1}}{C_{n}}\rho_{n}\implies\rho_{n}=\left(\prod\limits_{m=1}^{n}\frac{A_{m}}{C_{m-1}}\right)\rho_{0}.$
(62)

The initial $\rho_{0}$ is taken as 1 understanding that we must normalize the
probability distribution at the end of the calculation. Plugging in the forms
of the $A$ and $C$ coefficients, we have

$\rho_{n}=\frac{1}{Z}\prod\limits_{m=1}^{n}\frac{2r\epsilon^{2}/\kappa}{\Gamma^{2}+F(m)}\equiv\frac{1}{Z}\prod\limits_{m=1}^{n}\frac{\alpha}{1+G(m)}=\frac{\alpha^{n}}{Z}\prod\limits_{m=1}^{n}\frac{1}{1+G(m)},$
(63)

where we have introduced $\alpha=\frac{2r\epsilon^{2}}{\kappa\Gamma^{2}}$,
$G(m)=F(m)/\Gamma^{2}$, and
$Z=1+\sum\limits_{n=1}^{\infty}\left(\prod\limits_{m=1}^{n}\frac{A_{m}}{C_{m-1}}\right)$,
the normalization constant. We note that the factor $\alpha(1+G(n))^{-1}$
essentially “propagates” the probability distribution from $n$ to $n+1$. These
results underlie the results of Fig. 2 and Fig. 3 of the main text.

In Fig. S2, we expand upon Fig. 2 of the main text by showing the statistics
as a function of pump for different coupling parameters, to give the reader of
a clearer sense of the transition from thermal to coherent to Fock statistics.

Figure 5: Evolution of photon statistics with pumping: thermal, coherent,
anti-thermal, and Fock states. Photon probability distributions as a function
for different pump strengths and different coupling strengths. For weak
coupling, the statistics evolve from thermal to coherent with increasing pump.
For the largest couplings considered, the state evolves from thermal (for low
pump) to coherent (for intermediate pump) to a thermal state of _negative
temperature_ for higher pump. As the pump increases, the negative temperature
state converges effectively to a Fock state. Note that the bottom left panel
overlaps with Fig. 2 of the main text.

###  B.3 Summary of changes with $\lambda\neq 0$

In most of the discussion, we have assumed that $\lambda=0$. However, in many
superconducting qubit systems, a nonzero $\sigma_{x}$ term is present due to
an applied bias field which can tune the system. Our main results of
nonperbative nonlinearity, as well as the Fock lasing action in these systems,
are robust to the addition of this term. We summarize the main changes here.

We assume that the generalized Rabi Hamiltonian now takes the full form

$H=\frac{1}{2}(\omega\sigma_{z}+\lambda\sigma_{x})+\omega
a^{\dagger}a+\tilde{g}\sigma_{x}(a+a^{\dagger}).$ (64)

The spectrum is now approximately given by:

$E_{n,\sigma}=n\omega+\frac{\sigma}{2}\sqrt{(\omega D_{n})^{2}+\lambda^{2}},$
(65)

which corresponds to the eigenstates

$\displaystyle\ket{n,+}$
$\displaystyle=\cos(\theta/2)D^{\dagger}(g)\ket{+x,n}+\sin(\theta/2)D(g)\ket{-x,n}$
(66) $\displaystyle\ket{n,-}$
$\displaystyle=\sin(\theta/2)D^{\dagger}(g)\ket{+x,n}-\cos(\theta/2)D(g)\ket{-x,n},$
(67)

where the mixing angle $\theta$ is defined by $\tan(\theta)=\omega
D_{n}/\lambda$.

In principles, the modifications to the analysis of the laser action should
follow through new additions to the matrix elements which couple these
eigenstates. Specifically, we have

$\displaystyle\braket{n^{\prime},+}{(a+a^{\dagger})}{n,+}$
$\displaystyle=(\sqrt{n}\delta_{n^{\prime},n-1}+\sqrt{n+1}\delta_{n^{\prime},n+1})-2g\cos\theta\delta_{nn^{\prime}}$
(68) $\displaystyle\braket{n^{\prime},-}{(a+a^{\dagger})}{n,-}$
$\displaystyle=(\sqrt{n}\delta_{n^{\prime},n-1}+\sqrt{n+1}\delta_{n^{\prime},n+1})+2g\cos\theta\delta_{nn^{\prime}}$
(69) $\displaystyle\braket{n^{\prime},-}{(a+a^{\dagger})}{n,+}$
$\displaystyle=-2g\sin\theta\delta_{nn^{\prime}}.$ (70)

However, we see that the only new terms are only nonzero when the photon
number stays the same. Thus, the only modifications to the equations of motion
come from the eigen-energies. This means that the equations of motion derived
previously still hold valid, but with a new nonlinearity:

$F(n)=4n\epsilon^{2}+\frac{\omega^{2}}{4}\left(\sqrt{D_{n}^{2}+\lambda^{2}}-\sqrt{D_{n-1}^{2}+\lambda^{2}}\right)^{2}.$
(71)

###  B.4 Direct method for evolving the density matrix

In the previous section, we treated the emitter-field interaction as if
excited two-level systems were being injected into the system at rate $r$. We
also treated the interaction with the emitters as sequential: as if one
emitter interacts with the field at any given time, with probability
coefficients averaged over the emitter’s exponential decay probability. In
this section, we provide an alternative treatment of the problem in which we
consider the direct evolution of the density matrix in the presence of
coherent emitter-field interaction, emitter pumping, emitter decay, and field
leakage. This approach, besides being in principle more rigorous, and besides
providing further corroboration of our results above, also allows us to
consider multi-level emitter systems, such as three- and four-level systems,
which are more practical from the standpoint of lasers. This method has been
applied to describe conventional lasers (see scully1999quantum ), but due to
its generality, can be used to describe the Fock laser discussed in this
paper.

The equation of motion for the density matrix is

$\dot{\rho}=-i[H/\hbar,\rho]+\sum\limits_{i}\frac{\gamma_{i}}{2}\left(2J_{i}\rho
J_{i}^{\dagger}-J_{i}^{\dagger}J_{i}\rho-\rho
J_{i}^{\dagger}J_{i}\right)\equiv-i[H,\rho]+\sum\limits_{i}\mathcal{L}_{i}[\rho],$
(72)

where

$H=H_{\text{Rabi}}+\sum\limits_{i=1}^{N}H_{\text{em,i}}+\hbar\epsilon_{i}(\ket{a_{i}}\bra{b_{i}}+\ket{b_{i}}\bra{a_{i}})(b+b^{\dagger}),$
(73)

is the Hamiltonian describing $N$ multi-level emitters (with Hamiltonian
$H_{\text{em,i}}$) coupled to the electromagnetic field associated with matter
strongly coupled to a single electromagnetic field mode with coupling constant
$\epsilon_{i}$. The levels $a$ and $b$ of the $i$th emitter are coupled to the
field and comprise respectively the excited and ground levels of the lasing
transition. We have changed $e,g\rightarrow a,b$ as in what follows, we will
introduce more levels to incorporate realistic decay channels. The terms on
the right of the density matrix equation of motion are Lindblad terms with
decay rates $\gamma_{i}$ and jump operators $J_{i}$. The index $i$ enumerates
over the possible decay mechanisms, as well as all of the emitters. In what
follows, we review several (standard) simplifications of this equation that
render a readily solvable problem.

For simplicity, we will consider the case (as before) where all emitters have
the same energy levels (and thus the same $H_{\text{em,i}}$) and coupling
constant $\epsilon_{i}=\epsilon$ (which is the average coupling dictated by
the emitter distribution and cavity mode profile). Similarly, the decay
constants of all atomic levels are taken as the same. These simplifications do
not negate the effects reported here. Beyond these simplifications, a key
simplification arises because the emitter-field coupling couples all of the
emitters to a single quantum oscillator. In this case, we can consider the
problem as effectively a one-emitter problem where

$H=H_{\text{Rabi}}+\hbar\epsilon(\ket{e}\bra{g}+\ket{g}\bra{e})(a+a^{\dagger})\equiv
H_{0}+V,$ (74)

and the $i$ in the Lindblad terms enumerates only over decay channels. We have
defined for simplicity $H_{0}=H_{\text{Rabi}}+H_{\text{em}}$ and
$V/\hbar=\epsilon(\ket{e}\bra{g}+\ket{g}\bra{e})(a+a^{\dagger})$.

To start, we will consider decay channels for the emitter only, and not the
field, and include the field decay channels at the end of the calculation. In
what follows, we consider an emitter system consisting of lasing levels $a,b$,
ground level $g$, and “bath levels” $c,d$ for which $a$ and $b$ respectively
decay to. The pumping from $g\rightarrow a$ occurs with rate $r$, while the
$a\rightarrow c$ decay occurs with rate $\gamma_{a}$, the $b\rightarrow d$
decay occurs with rate $\gamma_{b}$, the $c\rightarrow g$ decay occurs with
rate $\gamma_{c}$, and the $d\rightarrow g$ decay occurs with rate
$\gamma_{d}$. Thus, the density matrix equation of motion may be written as

$\dot{\rho}=-i[H_{0}/\hbar,\rho]-i[V/\hbar,\rho]+\sum\limits_{i=g,a,b,c,d}\mathcal{L}_{i}[\rho],$
(75)

where the jump operators for $g,a,b,c,d$ are respectively
$\ket{a}\bra{g},\ket{c}\bra{a},\ket{d}\bra{b},\ket{g}\bra{c},\ket{g}\bra{d}$
with corresponding rates $r,\gamma_{a},\gamma_{b},\gamma_{c},\gamma_{d}$.

Let us now write an equation of motion for the matrix elements of the density
matrix, $\rho_{\beta n^{\prime},\alpha n}$, where $\alpha,\beta$ enumerate
over emitter states $g,a-d$ and the $n,n^{\prime}$ enumerate over the
eigenstates of the Rabi Hamiltonian (e.g., the Fock states of DSC photons). We
are considering the Hamiltonian only in one spin projection, as in the
previous treatment, since the spins decouple, both in the conventional Rabi
model, and the generalized one (with $\lambda\neq 0$). To proceed, we will
need the following matrix elements

$\displaystyle\bra{\beta,n^{\prime}}[H_{0}/\hbar,\rho]\ket{\alpha,n}$
$\displaystyle=(\omega_{\beta n^{\prime}}-\omega_{\alpha n})\rho_{\beta
n^{\prime},\alpha n}$ $\displaystyle\bra{b,n^{\prime}}V\rho\ket{\alpha,n}$
$\displaystyle=V_{bn^{\prime},an^{\prime}-1}\rho_{an^{\prime}-1,\alpha n}$
$\displaystyle\bra{a,n^{\prime}}V\rho\ket{\alpha,n}$
$\displaystyle=V_{an^{\prime},bn^{\prime}+1}\rho_{bn^{\prime}+1,\alpha n}$
$\displaystyle\bra{\beta,n^{\prime}}\rho V\ket{b,n}$
$\displaystyle=\rho_{\beta n^{\prime},an-1}V_{an-1,bn}$
$\displaystyle\bra{\beta,n^{\prime}}\rho V\ket{a,n}$
$\displaystyle=\rho_{\beta n^{\prime},bn+1}V_{bn+1,an}.$ (76)

For matrix elements of $V\rho$ and $\rho V$, we have used the structure of the
matrix elements in the section ”Time-evolution of the coupled system”, where
we showed that the effect of the coupling is to change the emitter state, and
to change the number of field quanta by one.

We also need the matrix elements of the Lindblad terms. Let us consider a
generic Lindblad term of the form

$\bra{\beta,n^{\prime}}\mathcal{L}_{i}[\rho]\ket{\alpha,n}=\frac{\gamma_{i}}{2}\bra{\beta,n^{\prime}}2T_{ji}\rho
T^{\dagger}_{ji}-T^{\dagger}_{ji}T_{ji}\rho-\rho
T^{\dagger}_{ji}T_{ji}\ket{\alpha,n},$ (77)

where $T_{ji}=\ket{j}\bra{i}$. $T_{ij}$ is simply $J_{i}$ with the final-state
index $j$ included for clarity. The matrix element follows as

$\bra{\beta,n^{\prime}}\mathcal{L}_{i}[\rho]\ket{\alpha,n}=\frac{\gamma_{i}}{2}(2\delta_{j\beta}\delta_{j\alpha}\rho_{in^{\prime},in}-\delta_{i\alpha}\rho_{\beta
n^{\prime},in}-\delta_{i\beta}\rho_{in^{\prime},\alpha n}).$ (78)

With these matrix elements tabulated, we may write the following set of
equations for the matrix elements of the density matrix:

$\displaystyle\dot{\rho}_{an^{\prime},an}$
$\displaystyle=-i\omega_{n^{\prime}n}\rho_{an^{\prime},an}-\gamma_{a}\rho_{an^{\prime},an}+r\rho_{gn^{\prime},gn}-i(V_{an^{\prime},bn^{\prime}+1}\rho_{bn^{\prime}+1,an}-\rho_{an^{\prime},bn+1}V_{bn+1,an})$
$\displaystyle\dot{\rho}_{bn^{\prime}+1,an}$
$\displaystyle=\left[-i(\omega_{bn^{\prime}+1}-\omega_{an})-\frac{\gamma_{a}+\gamma_{b}}{2}\right]\rho_{bn^{\prime}+1,an}-i(V_{bn^{\prime}+1,an^{\prime}}\rho_{an^{\prime},an}-\rho_{bn^{\prime}+1,bn+1}V_{bn+1,an})$
$\displaystyle\dot{\rho}_{an^{\prime},bn+1}$
$\displaystyle=\left[-i(\omega_{an^{\prime}}-\omega_{bn+1})-\frac{\gamma_{a}+\gamma_{b}}{2}\right]\rho_{an^{\prime},bn+1}-i(V_{an^{\prime},bn^{\prime}+1}\rho_{bn^{\prime}+1,bn+1}-\rho_{an^{\prime},an}V_{an,bn+1})$
$\displaystyle\dot{\rho}_{bn^{\prime}+1,bn+1}$
$\displaystyle=-i\omega_{n^{\prime}+1,n+1}\rho_{bn^{\prime}+1,bn+1}-\gamma_{b}\rho_{bn^{\prime}+1,bn+1}-i(V_{bn^{\prime}+1,an^{\prime}}\rho_{an^{\prime},bn+1}-\rho_{bn^{\prime}+1,an}V_{an,bn+1})$
$\displaystyle\dot{\rho}_{cn^{\prime},cn}$
$\displaystyle=(-i\omega_{n^{\prime}n}-\gamma_{c})\rho_{cn^{\prime},cn}+\gamma_{a}\rho_{an^{\prime},an}$
$\displaystyle\dot{\rho}_{dn^{\prime},dn}$
$\displaystyle=(-i\omega_{n^{\prime}n}-\gamma_{d})\rho_{dn^{\prime},dn}+\gamma_{b}\rho_{bn^{\prime},bn}$
$\displaystyle\dot{\rho}_{gn^{\prime},gn}$
$\displaystyle=(-i\omega_{n^{\prime}n}-r)\rho_{gn^{\prime},gn}+\gamma_{c}\rho_{cn^{\prime},cn}+\gamma_{d}\rho_{dn^{\prime},dn}.$
(79)

While these equations can be generally solved, we focus as in the previous
treatment on the steady state dynamics. As expected from conventional lasers,
the steady state density matrix is diagonal due to decoherence. Numerically,
for this laser system, based on deep strong light-matter coupling, we also
found that the steady-state (found by the null eigenvector of the Liouvillian
($S$ such that $\dot{\rho}=S\rho$)) is diagonal. Let us thus focus on the
steady-state equations for the “photon diagonals” ($n=n^{\prime}$), which are
simply

$\displaystyle 0$
$\displaystyle=r\rho_{gn,gn}-\gamma_{a}\rho_{bn,an}-i(V^{*}\rho_{bn+1,an}-\rho_{an,bn+1}V)$
$\displaystyle 0$
$\displaystyle=\left[i\Delta_{n+1}-\frac{\gamma_{a}+\gamma_{b}}{2}\right]\rho_{bn+1,an}-i(V\rho_{an,an}-\rho_{bn+1,bn+1}V)$
$\displaystyle 0$
$\displaystyle=\left[-i\Delta_{n+1}-\frac{\gamma_{a}+\gamma_{b}}{2}\right]\rho_{an,bn+1}-i(V^{*}\rho_{bn+1,bn+1}-\rho_{an,an}V^{*})$
$\displaystyle 0$
$\displaystyle=-\gamma_{b}\rho_{bn+1,bn+1}-i(V\rho_{an,bn+1}-\rho_{bn+1,an}V^{*})$
$\displaystyle 0$
$\displaystyle=-\gamma_{c}\rho_{cn,cn}+\gamma_{a}\rho_{an,an}$ $\displaystyle
0$ $\displaystyle=-\gamma_{d}\rho_{dn,dn}+\gamma_{b}\rho_{bn,bn}$
$\displaystyle 0$
$\displaystyle=-r\rho_{gn,gn}+\gamma_{c}\rho_{cn,cn}+\gamma_{d}\rho_{dn,dn},$
(80)

where we have defined $\Delta_{n+1}=\omega_{an}-\omega_{bn+1}$.

Immediately, we have $\gamma_{a}\rho_{an,an}=\gamma_{c}\rho_{cn,cn}$ and
$\gamma_{b}\rho_{bn,bn}=\gamma_{d}\rho_{dn,dn}$. The equation for
$\rho_{gn,gn}$ then can be written as

$r\rho_{gn}=\gamma_{a}\rho_{an,an}+\gamma_{b}\rho_{bn,bn}.$ (81)

For simplicity, let us take $\gamma_{a}=\gamma_{b}=\Gamma$, so that

$r\rho_{gn}=\Gamma(\rho_{an,an}+\rho_{bn,bn})=\Gamma(\rho_{nn}-\rho_{cn,cn}-\rho_{dn,dn}-\rho_{gn,gn}),$
(82)

where we have defined the photon populations
$\rho_{nn}=(\text{tr}_{\text{em}}\rho)_{nn}$ in order to express everything in
terms of these populations and arrive at a coarse-grained density matrix for
the field. Let us now consider the case where $\gamma_{c}\gg\gamma_{a}$ and
$\gamma_{d}\gg\gamma_{b}$. In this case, we immediately see that
$\rho_{cn,cn}\approx 0$ and $\rho_{dn,dn}\approx 0$. This is to say that these
levels are depleted immediately after they are populated by the lasing levels.
In this case, $\rho_{gn,gn}=\frac{\Gamma}{(r+\Gamma)}\rho_{nn}$. The steady-
state equations then reduce to the simple inhomogeneous equation:

$\begin{bmatrix}\Gamma&-iV&iV^{*}&0\\\
-iV^{*}&\Gamma+i\Delta_{n+1}&0&iV^{*}\\\ iV&0&\Gamma-i\Delta_{n+1}&-iV\\\
0&iV&-iV^{*}&\Gamma\\\ \end{bmatrix}\begin{bmatrix}\rho_{an,an}\\\
\rho_{an,bn+1}\\\ \rho_{bn+1,an}\\\ \rho_{bn+1,bn+1}\\\
\end{bmatrix}=\frac{r\Gamma\rho_{nn}}{r+\Gamma}\begin{bmatrix}1\\\ 0\\\ 0\\\
0\\\ \end{bmatrix}\equiv r_{a}\rho_{nn}\begin{bmatrix}1\\\ 0\\\ 0\\\ 0\\\
\end{bmatrix},$ (83)

whose solution yields
$\rho_{an,an},\rho_{an,bn+1},\rho_{bn+1,an},\rho_{bn+1,bn+1}$ in terms of
$\rho_{nn}$. To proceed most efficiently, we now connect these density matrix
elements to the equation of motion for the reduced density matrix of the
field. This equation of motion is

$\dot{\rho}_{nn}=-i(\rho_{an-1,bn}V_{bn,an-1}+V_{an,bn+1}\rho_{bn+1,an}-\rho_{bn,an-1}V_{an-1,bn}-V_{bn+1,an}\rho_{an,bn+1}).$
(84)

The coherences can be found from the matrix equation above, which for brevity,
we denote as $M_{n}P_{n}=r_{a}\rho_{nn}e_{1}$ so that
$P_{n}=r_{a}\rho_{nn}(M_{n}^{-1}e_{1})$. We may write this equation in a form
similar to that of the equation for the coarse-grained density matrix of the
previous treatment, i.e., as

$0=A_{n}\rho_{n-1,n-1}-A_{n+1}\rho_{n,n},$ (85)

where
$A_{n}=-ir_{a}((e^{T}_{2}M_{n-1}^{-1}e_{1})V_{bn,an-1}-(e^{T}_{3}M_{n-1}^{-1}e_{1})V_{an-1,bn})$,
and we have looked at the steady-state limit. At this stage, we now add the
Lindblad terms corresponding to the cavity leakage. As per the discussion in
the section “Lindblad terms”, the resulting equation of motion for the field
density matrix is

$\dot{\rho}_{nn}=A_{n}\rho_{n-1,n-1}-(A_{n+1}+\kappa|\langle
n-1|a+a^{\dagger}|n\rangle|^{2})\rho_{n,n}+\kappa|\langle
n|a+a^{\dagger}|n+1\rangle|^{2}\rho_{n+1,n+1}.$ (86)

Defining $B_{n}=-(A_{n+1}+\kappa|\langle n-1|a+a^{\dagger}|n\rangle|^{2})$ and
$C_{n}=\kappa|\langle n|a+a^{\dagger}|n+1\rangle|^{2}$, we have as before
$A_{n+1}+B_{n}+C_{n-1}=0$, enabling us to immediately write

$\rho_{n+1}=\frac{A_{n+1}}{C_{n}}\rho_{n}\implies\rho_{n}=\left(\prod\limits_{m=1}^{n}\frac{A_{m}}{C_{m-1}}\right)\rho_{0}.$
(87)

with the initial $\rho_{0}$ is taken as one with the understanding that we
must normalize the probability distribution at the end of the calculation.
Upon inversion of the matrix $M_{n}$, we immediately find

$A_{n}=\frac{2r_{a}|V_{bn,an-1}|^{2}}{\Gamma^{2}+4|V_{bn,an-1}|^{2}+\Delta_{n}^{2}}.$
(88)

Noting that $|V_{bn,an-1}|^{2}=\epsilon^{2}|\langle
n-1|a+a^{\dagger}|n\rangle|^{2}$, we may write the overall equation as

$\rho_{n}=\frac{1}{Z}\left(\prod\limits_{m=1}^{n}\frac{2r_{a}\epsilon^{2}/\kappa\Gamma^{2}}{1+(4\epsilon^{2}|\langle
m-1|a+a^{\dagger}|m\rangle|^{2}+\Delta_{m}^{2})/\Gamma^{2}}\right)\equiv\frac{\alpha^{n}}{Z}\left(\prod\limits_{m=1}^{n}\frac{1}{1+G(m)}\right),$
(89)

where $\alpha=2r_{a}\epsilon^{2}/\kappa\gamma^{2}$ and
$Z=1+\sum\limits_{n=1}^{\infty}\left(\prod\limits_{m=1}^{n}\frac{A_{m}}{C_{m-1}}\right)$.
Immediately, we see that if we take $|\langle
n-1|a+a^{\dagger}|n\rangle|^{2}=n$ and
$\Delta_{n}^{2}=\frac{1}{4}\omega^{2}e^{-4g^{2}}(L_{n}(4g^{2})-L_{n-1}(4g^{2}))^{2}$
(assuming $\omega_{0}=\omega$) that we recover the results of the previous
treatment. And it may also be easily seen that this agreement persists if we
take the matrix elements and splitting to be governed by the generalized Rabi
model (with $\lambda\neq 0$).

##  Appendix C Numerical calculation of the Fock laser steady state

In this section, we numerically validate the analytical developments of the
previous sections. Since the analytical calculations make use of many
approximations and assumptions, it is important to validate them in terms of a
method which is independent of these assumptions. In what follows, we will use
a method inspired by the observation that the equation of motion for the laser
density matrix effectively describes the interaction of a single gain atom
with a cavity, even when the gain medium is composed of many atoms. This is
because the atoms only couple to each other through the cavity field, as noted
in scully1999quantum . Thus, it follows that laser steady states can be
understood through the steady state of the Liouvillian operator describing a
damped oscillator coupled to a gain atom. By taking the partial trace of the
null eigenvector of the Liouvillian, one finds the steady state probability
distribution of DSC photons. Applied to conventional lasers, one correctly
finds the transition from thermal to coherent state statistics above the laser
threshold.

The Hamiltonian part of the Liouvillian is simply the Hamiltonian of Eq. (2).
Here, we also consider two different possible interaction terms:
$a+a^{\dagger}$ or $b+b^{\dagger}$. The steady-states are also insensitive to
this. In these calculations, we include a reservoir to describe pumping of the
gain medium, as well as its $T_{1}$\- and $T_{2}$-relaxation (here,
$T_{2}=2T_{1}$). We also include a reservoir to describe the decay of the DSC
photon. We take as the jump operator $J^{(+)}$, where $J=a+a^{\dagger}$ or
$b+b^{\dagger}$ and the $+$ superscript means “positive-frequency”, meaning
that we project out all negative frequency components. The steady-state is
insensitive to whether we use $a$ or $b$ (it changes only slightly),
indicating the relative unimportance of the dipole-dipole interaction from the
standpoint of the photon probabilities. The overall Liouvillian is then the
sum of the Hamiltonian (commutator) part and three dissipators:
$\mathcal{D}[\sigma_{\text{em}}^{(-)}],\mathcal{D}[\sigma_{\text{em}}^{(+)}]$
and $\mathcal{D}[J^{(+)}]$ with respective rates $r$, $\Gamma$, and $\kappa$,
to describe gain pumping, gain decoherence, and DSC photon decay. Note that
$\mathcal{D}[O]\rho\equiv O^{\dagger}O\rho+\rho O^{\dagger}O-2O\rho
O^{\dagger}$.

In Fig. S3, we show the probability distribution of DSC photons resulting from
one of these steady-state calculations. Because we have largely neglected spin
in our analytical discussions, we plot the “unpolarized” photon probability
distribution, defined such that $P(n)=P(n,-1)+P(n,1)$. As can be seen, above
threshold, the state has very low number fluctuations, in this case, $\delta
n$ = 1, yielding a state very close to a Fock state. Below threshold, a
thermal state is found. Around threshold, the quasi-uniform state of Fig. 3 of
the main text is found.

These results are insensitive to the presence of a $\lambda$-term, as shown in
Fig. S4. The presence of a $\lambda$ term, all else equal, slightly increases
the photon noise. This is because the presence of a $\lambda$ softens the
anharmonicity (which can be understood from the term
$\sqrt{\lambda^{2}+\omega^{2}e^{-4g^{2}}L^{2}_{n}(4g^{2})}$ in Eq. (60)).

The results are also insensitive to the exact form of interaction and
dissipator (provided that the dissipator doesn’t create spurious excitations).
In Fig. S5, we show the steady-state, computed using interaction terms based
on $a$ or $b$, as well as dissipators based on $a$ or $b$.

Figure 6: Steady-state of the Fock laser, calculated numerically, by finding
the steady-state of the Liouvillian operator. Plot shows the unpolarized
probability distribution for $\lambda=0$. For $\epsilon=10^{-5}\omega$,
$\kappa=10^{-8}\omega$, and $r=10\Gamma$ (such that the population inversion
of the gain is about 90%, the resulting state is nearly a Fock state of 100
DSC photons, with a residual uncertainty of 1. This state has noise 99% below
the shot noise level. Moreover, this calculation shows that the Hamiltonian of
Eq. (3), coupled to damping, supports Fock states as its steady state, from
first principles.  Figure 7: Steady-state of the Liouvillian of the Fock laser
with and without the $\lambda$-term of Eq. (1) of the main text. Here,
$\lambda=0.1\omega$ is sufficient to keep the two spin ladders from
interchanging, and is not found to alter the steady-state appreciably.
$\epsilon,\Gamma,\kappa$ are the same as in the above figure.  Figure 8:
Steady-state of the Liouvillian of the Fock laser with different choices for
the interaction term and dissipator, showing robustness to the exact magnitude
of the dipole-dipole interaction between the emitter and the qubit.

## References

  * [1] Anton Frisk Kockum, Adam Miranowicz, Simone De Liberato, Salvatore Savasta, and Franco Nori.  Ultrastrong coupling between light and matter.  Nature Reviews Physics, 1(1):19–40, 2019.
  * [2] James A Hutchison, Tal Schwartz, Cyriaque Genet, Eloïse Devaux, and Thomas W Ebbesen.  Modifying chemical landscapes by coupling to vacuum fields.  Angewandte Chemie International Edition, 51(7):1592–1596, 2012\.
  * [3] Johannes Flick, Nicholas Rivera, and Prineha Narang.  Strong light-matter coupling in quantum chemistry and quantum photonics.  Nanophotonics, 7(9):1479–1501, 2018.
  * [4] Michael Ruggenthaler, Nicolas Tancogne-Dejean, Johannes Flick, Heiko Appel, and Angel Rubio.  From a quantum-electrodynamical light–matter description to novel spectroscopies.  Nature Reviews Chemistry, 2(3):1–16, 2018.
  * [5] Anoop Thomas, Lucas Lethuillier-Karl, Kalaivanan Nagarajan, Robrecht MA Vergauwe, Jino George, Thibault Chervy, Atef Shalabney, Eloïse Devaux, Cyriaque Genet, Joseph Moran, et al.  Tilting a ground-state reactivity landscape by vibrational strong coupling.  Science, 363(6427):615–619, 2019.
  * [6] E Orgiu, J George, JA Hutchison, E Devaux, JF Dayen, B Doudin, F Stellacci, C Genet, J Schachenmayer, Claudiu Genes, et al.  Conductivity in organic semiconductors hybridized with the vacuum field.  Nature Materials, 14(11):1123–1129, 2015.
  * [7] Gian L Paravicini-Bagliani, Felice Appugliese, Eli Richter, Federico Valmorra, Janine Keller, Mattias Beck, Nicola Bartolo, Clemens Rössler, Thomas Ihn, Klaus Ensslin, et al.  Magneto-transport controlled by landau polariton states.  Nature Physics, 15(2):186–190, 2019.
  * [8] P Forn-Díaz, L Lamata, E Rico, J Kono, and E Solano.  Ultrastrong coupling regimes of light-matter interaction.  Reviews of Modern Physics, 91(2):025005, 2019.
  * [9] Fumiki Yoshihara, Tomoko Fuse, Sahel Ashhab, Kosuke Kakuyanagi, Shiro Saito, and Kouichi Semba.  Superconducting qubit–oscillator circuit beyond the ultrastrong-coupling regime.  Nature Physics, 13(1):44–47, 2017.
  * [10] P Forn-Díaz, J José García-Ripoll, Borja Peropadre, J-L Orgiazzi, MA Yurtalan, R Belyansky, Christopher M Wilson, and A Lupascu.  Ultrastrong coupling of a single artificial atom to an electromagnetic continuum in the nonperturbative regime.  Nature Physics, 13(1):39–43, 2017.
  * [11] Simone De Liberato.  Light-matter decoupling in the deep strong coupling regime: The breakdown of the purcell effect.  Physical review letters, 112(1):016401, 2014.
  * [12] Nicholas Rivera, Johannes Flick, and Prineha Narang.  Variational theory of nonrelativistic quantum electrodynamics.  Physical review letters, 122(19):193603, 2019.
  * [13] Jorge Casanova, Guillermo Romero, Ion Lizuain, Juan José García-Ripoll, and Enrique Solano.  Deep strong coupling regime of the jaynes-cummings model.  Physical review letters, 105(26):263603, 2010.
  * [14] Fumiki Yoshihara, Tomoko Fuse, Ziqiao Ao, Sahel Ashhab, Kosuke Kakuyanagi, Shiro Saito, Takao Aoki, Kazuki Koshino, and Kouichi Semba.  Inversion of qubit energy levels in qubit-oscillator circuits in the deep-strong-coupling regime.  Physical review letters, 120(18):183601, 2018.
  * [15] EK Irish.  Generalized rotating-wave approximation for arbitrarily large coupling.  Physical review letters, 99(17):173601, 2007.
  * [16] 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\.
  * [17] Daniel F Walls and Gerard J Milburn.  Quantum optics.  Springer Science & Business Media, 2007.
  * [18] Nicholas Rivera and Ido Kaminer.  Light–matter interactions with photonic quasiparticles.  Nature Reviews Physics, 2(10):538–561, 2020.
  * [19] Alessandro Ridolfo, Martin Leib, Salvatore Savasta, and Michael J Hartmann.  Photon blockade in the ultrastrong coupling regime.  Physical review letters, 109(19):193602, 2012.
  * [20] Alexandre Le Boité, Myung-Joong Hwang, Hyunchul Nha, and Martin B Plenio.  Fate of photon blockade in the deep strong-coupling regime.  Physical Review A, 94(3):033827, 2016.
  * [21] Marlan O Scully and M Suhail Zubairy.  Quantum optics, 1999.
  * [22] Jason McKeever, Andreea Boca, A David Boozer, Joseph R Buck, and H Jeff Kimble.  Experimental realization of a one-atom laser in the regime of strong coupling.  Nature, 425(6955):268–271, 2003.
  * [23] Kyungwon An, James J Childs, Ramachandra R Dasari, and Michael S Feld.  Microlaser: A laser with one atom in an optical resonator.  Physical review letters, 73(25):3375, 1994.
  * [24] Y-Y Liu, J Stehlik, Christopher Eichler, MJ Gullans, Jacob M Taylor, and JR Petta.  Semiconductor double quantum dot micromaser.  Science, 347(6219):285–287, 2015.
  * [25] J Stehlik, Y-Y Liu, C Eichler, TR Hartke, X Mi, MJ Gullans, JM Taylor, and Jason R Petta.  Double quantum dot floquet gain medium.  Physical Review X, 6(4):041027, 2016.
  * [26] Guido Burkard, Michael J Gullans, Xiao Mi, and Jason R Petta.  Superconductor–semiconductor hybrid-circuit quantum electrodynamics.  Nature Reviews Physics, 2(3):129–140, 2020.
  * [27] MH Devoret, Steven Girvin, and Robert Schoelkopf.  Circuit-qed: How strong can the coupling between a josephson junction atom and a transmission line resonator be?  Annalen der Physik, 16(10-11):767–779, 2007.
  * [28] Nicholas Rivera, Jamison Sloan, Yannick Salamin, and Marin Soljacic.  Macroscopic condensation of photon noise in sharply nonlinear dissipative systems.  arXiv preprint arXiv:2111.03099, 2021.
  * [29] Félix Beaudoin, Jay M Gambetta, and A Blais.  Dissipation and ultrastrong coupling in circuit qed.  Physical Review A, 84(4):043832, 2011.
  * [30] Tuomas Jaako, Ze-Liang Xiang, Juan José Garcia-Ripoll, and Peter Rabl.  Ultrastrong-coupling phenomena beyond the dicke model.  Physical Review A, 94(3):033850, 2016.
  * [31] Daniele De Bernardis, Philipp Pilar, Tuomas Jaako, Simone De Liberato, and Peter Rabl.  Breakdown of gauge invariance in ultrastrong-coupling cavity qed.  Physical Review A, 98(5):053819, 2018.
  * [32] Daniele De Bernardis, Tuomas Jaako, and Peter Rabl.  Cavity quantum electrodynamics in the nonperturbative regime.  Physical Review A, 97(4):043820, 2018.
  * [33] Omar Di Stefano, Alessio Settineri, Vincenzo Macrì, Luigi Garziano, Roberto Stassi, Salvatore Savasta, and Franco Nori.  Resolution of gauge ambiguities in ultrastrong-coupling cavity quantum electrodynamics.  Nature Physics, 15(8):803–808, 2019.
  * [34] Alessio Settineri, Omar Di Stefano, David Zueco, Stephen Hughes, Salvatore Savasta, and Franco Nori.  Gauge freedom, quantum measurements, and time-dependent interactions in cavity and circuit qed.  arXiv preprint arXiv:1912.08548, 2019.