Probing dynamical symmetry breaking using quantum-entangled photons

The interaction between light and matter lies at the heart of all photophysics and spectroscopy. Typically, one treats the interaction within a semi-classical approximation, treating light as an oscillating classical electro-magnetic wave as given by Maxwell's equations. It is well recognized that light has a quantum mechanical discreteness (photons) and one can prepare entangled interacting photon states. The pioneering work by Hanbury Brown and Twiss in the 1950's, who measured intensity correlations in light originating from thermal sources, set the stage for what has become quantum optics [1–7]. Quantum photons play a central role in a number of advanced technologies including quantum cryptography [8], quantum communications [9], and quantum computation [10, 11]. Only recently has it been proposed that entangled photons can be exploited as a useful spectroscopic probe of atomic and molecular processes [11–16].

The spectral and temporal nature of entangled photons offer a unique means for interrogating the dynamics and interactions between molecular states. The crucial consideration is that when entangled photons are created, typically by spontaneous parametric down-conversion (SPDC), there is a precise relation between the frequency and wavevectors of the entangled pair. For example if we create two entangled photons from a common laser source, energy conservation dictates that ${\omega }_{{\rm{laser}}}={\omega }_{1}+{\omega }_{2}$. Hence measuring the frequency of either photon will collapse the quantum entanglement and the frequency of the other photon will be precisely defined. Moreover, in the case of multi-photon absorption, entangled 2-photon absorption is greatly enhanced relative to classical 2-photon absorption since the cross-section scales linearly rather than quadratically with intensity. Recent work by Schlawin et al indicate that entangled photon pairs may be useful in controlling and manipulating population on the 2-exciton manifold of a model biological energy transport system [17]. The non-classical features of entangled photons have also been used as a highly sensitive detector of ultra-fast emission from organic materials [18]. Beyond the potential practical applications of quantum light in high-fidelity communication and quantum encryption, by probing systems undergoing spontaneous symmetry breaking with quantum photons one can draw analogies between bench-top laboratory based experiments and experimentally inaccessible systems such as black holes, the early Universe, and cosmological strings [19, 20].

We begin with a brief overview of the photon coincidence experiment and the preparation of two-photon entangled states, termed 'Bell-states'. We then use the input/output approach of Gardner and Collett [21] to develop a means for computing the transmission function for a material system placed in one of the arms of the Hong–Ou–Mandel (HOM) apparatus sketched in figure 1. We provide a precise connection between a material sample, described in terms of a model Hamiltonian, and the resulting signal in the context of the interference experiment described by Kalashnikov et al in [16] using the Dicke model for an ensemble of two-level atoms as input [22]. The Dicke model is an important test-case for a variety of quantum effects occurring in photonic and plasmonic cavities [23]. It allows for a non-trivial steady state identified as a non-equilibrium phase transition to a super-radiant regime. The model is especially important since an analytical solution for this regime exists. For this reason, we have adopted the Dicke model in our calculations.

Zoom In Zoom Out Reset image size

Our goal is to distinguish between the elementary excitations (i.e. the fluctuations) above the super-radiant cooperative state and their signatures in the photon coincidence measurements. This is in contrast to conventional techniques probing elementary excitations above the ground states of an ensemble of quantum emitters. In a broader context, one can also model Frenkel excitons in molecular aggregates as a series of two-level systems and placed into a homogeneous electric field as a modified Dicke model thereby providing connection to realistic material systems [24, 25]. Lastly, in realistic experimental situations some deviations are expected such as field inhomogeneity, energetic disorder, appearance of higher-lying excited states (e.g., multi-level systems), cascades, and complex many-body interactions. Such higher order interactions can be identified by comparison with the predictions based upon the Dicke model.

We consider the interferometric scheme implemented by Kalashnikovet al [16]. A CW laser beam is incident on a nonlinear crystal, creating an entangled photon pair state by SPDC, which we shall denote as a Bell state,

where ${ \mathcal F }({\omega }_{1},{\omega }_{2})$ is the bi-photon field amplitude and ${B}_{S,I}^{\dagger }({\omega }_{i})$ creates a photon with frequency ${\omega }_{i}$ in either the signal or idler branch. The ket $| 0\rangle$ is the vacuum state and $| {\omega }_{1}{\omega }_{2}\rangle$ denotes a two photon state. In general, energy conservation requires that the entangled photons generated by SPDC obey ${\omega }_{L}={\omega }_{1}+{\omega }_{2}$. Similarly, conservation of photon momentum requires ${{\bf{k}}}_{L}={{\bf{k}}}_{1}+{{\bf{k}}}_{2}$. By manipulating the SPDC crystal, one can generate entangled photon pairs with different frequencies. As a result, the bi-photon field is strongly anti-correlated in frequency with

where ${\omega }_{L}$ is the central frequency of the bi-photon field. This aspect was recently exploited in [15], which used a visible photon in the idler branch and an infrared (IR) photon in the signal branch, interacting with the sample.

As sketched in figure 1, both signal and idler are reflected back towards a beam-splitter (BS) by mirrors M1 and M2. M1 introduces an optical delay with transmission function ${\rm{\Phi }}(\omega )$ which we will take to be of modulo 1. In the other arm, we introduce a resonant medium at S with transmission function ${ \mathcal S }(\omega )$. Not shown in our sketch is an optional pumping laser for creating a steady state exciton density in S. For the case of an isolated line with resonance frequency at ${\omega }_{o}$, the transmission function ${ \mathcal S }(\omega )$ can be written as

$$\begin{eqnarray}&&{ \mathcal S }(\omega )=\exp \left[-{\rm{i}}\displaystyle \frac{b}{(\omega -{\omega }_{o})+{\rm{i}}\gamma }\right],\end{eqnarray} \tag{ 3 }$$

where $b=\alpha L/2{T}_{2}$, $\alpha L$ is the optical thickness (α is a Bouger coefficient), and ${\gamma }^{-1}={T}_{2}$ is the dephasing time of the medium. The motivation for this paper is to provide a framework for a general system with more complex dynamics.

Upon interacting with both the delay element and the medium, the Bell-state can be rewritten as

Finally, the two beams are re-joined by a BS producing the mapping [26]

$$\begin{eqnarray}&&{B}_{I}^{\dagger }({\omega }_{1}){B}_{S}^{\dagger }({\omega }_{2})\mapsto \displaystyle \frac{1}{2}[{A}_{1}^{\dagger }({\omega }_{1})+{{\rm{i}}{A}}_{2}^{\dagger }({\omega }_{1})][{A}_{2}^{\dagger }({\omega }_{2})+{{\rm{i}}{A}}_{1}^{\dagger }({\omega }_{2})],\end{eqnarray} \tag{ 5 }$$

whereby ${A}_{i}^{\dagger }({\omega }_{j})$ creates a photon with frequency ${\omega }_{j}$ in the $i{\rm{t}}{\rm{h}}$ exit channel. This yields a final Bell state:

The first two terms in this state correspond to the two possible outcomes where one photon is transmitted into each of the outgoing channels. In other words, the signal photon is transmitted to detector C1 and the idler photon is transmitted to detector C2 or vice versa. The other two terms correspond to the cases where both signal and idler branch photons are transmitted to either C1 or C2. The coincidence detection discriminates the former from the latter and we write the coincident term as

We can then take the counting rate as proportional to the probability

If we assume that the delay stage is dispersionless with ${\rm{\Phi }}(\omega )={{\rm{e}}}^{{\rm{i}}\omega t}$ and take the symmetric case of ${\omega }_{1}={\omega }_{L}+z$ and ${\omega }_{2}={\omega }_{L}-z$,

$$\begin{eqnarray}{P}_{c}({t}_{{\rm{delay}}}) & = & \displaystyle \frac{1}{4}{\displaystyle \int }_{-\infty }^{+\infty }{\rm{d}}{z}| { \mathcal F }(z){| }^{2}\{| { \mathcal S }({\omega }_{L}-z){| }^{2}+| { \mathcal S }({\omega }_{L}+z){| }^{2}\\ & & -2\mathrm{Re}[{{ \mathcal S }}^{* }({\omega }_{L}-z){ \mathcal S }({\omega }_{L}+z){{\rm{e}}}^{-2{{\rm{i}}{z}{t}}_{{\rm{delay}}}}]\},\end{eqnarray} \tag{ 9 }$$

where ${t}_{\mathrm{delay}}$ is the time lag between entangled photons traversing the upper and lower arms of the HOM apparatus. In the absence of a sample, the coincidence count is exactly equal to zero when z = 0. It is important to note that two photons traversing the idler and signal branches need to be of different frequencies in order to have any observable effect. In fact, taking the limit that the bi-photon amplitude is extremely narrow about ${\omega }_{L}$ gives ${P}_{c}({t}_{\mathrm{delay}})=0$, which is the Hong–Ou–Mandel effect [27].

A crucial component of our approach is the action of the sample at S which introduces a transmission function $S(\omega )$ into the final Bell state. We wish to connect this function to the dynamics and molecular interactions within the sample. To accomplish this, we use the input/output formulation of quantum optics and apply this to an ensemble of identical 2-level states coupled to a common photon mode [28]. The technical details of our approach are presented in the appendices of this paper. In short, we begin with a description of the material system described by N two-level spin states coupled to common set of photon cavity modes.

$$\begin{eqnarray}{\hat{H}}_{{\rm{sys}}} & = & \displaystyle \sum _{j}\displaystyle \frac{{\rm{\hslash }}{\omega }_{o}}{2}{\hat{\sigma }}_{z,j}+\displaystyle \sum _{k}{\rm{\hslash }}({\omega }_{k}-{\rm{i}}\kappa ){\hat{\psi }}_{k}^{\dagger }{\hat{\psi }}_{k}\\ & & +\displaystyle \sum _{k,j}\displaystyle \frac{{\rm{\hslash }}{\lambda }_{{kj}}}{\sqrt{N}}({\hat{\psi }}_{k}^{\dagger }+{\hat{\psi }}_{k})({\hat{\sigma }}_{j}^{+}+{\hat{\sigma }}_{j}^{-}),\end{eqnarray} \tag{ 10 }$$

where $\{{\hat{\sigma }}_{z,j},{\hat{\sigma }}_{j}^{\pm }\}$ are local spin-1/2 operators for site j, ${\rm{\hslash }}{\omega }_{o}$ is the local excitation energy, and ${\lambda }_{{kj}}$ is the coupling between the kth photon mode and the jth site, which we will take to be uniform over all sites. We introduce κ as the decay rate of a cavity photon. We allow the photons in the cavity (S) to exchange quanta with photons in the HOM apparatus and derive the Heisenberg equations of motion corresponding to input and output photon fields within a steady-state assumption. This allows us to compute (appendix A) the scattering matrix connecting an incoming photon with frequency ν from the field to an outgoing photon with frequency ν returned to the field viz.

$$\begin{eqnarray}&&{{\rm{\Psi }}}_{{\rm{out}}}(\nu )=-{\hat{{\rm{\Omega }}}}_{{\rm{out}}}^{(-)\dagger }{\hat{{\rm{\Omega }}}}_{{\rm{in}}}^{(+)}{{\rm{\Psi }}}_{{\rm{in}}}(\nu ),\end{eqnarray} \tag{ 11 }$$

where ${\hat{{\rm{\Omega }}}}_{{\rm{in}},{\rm{out}}}^{(\pm )}$ are Møller operators that propagate an incoming (or outgoing) state from $t\to -\infty$ to t = 0 where it interacts with the sample or from t = 0 to an outgoing (or incoming) state at $t\to +\infty$ and give the S-matrix in the form of a response function

$$\begin{eqnarray}&&{ \mathcal S }(\nu )=\langle \delta {{\rm{\Psi }}}_{{\rm{out}}}^{\dagger }(\nu )\delta {{\rm{\Psi }}}_{{\rm{out}}}(\nu ^{\prime} )\rangle \delta (\nu -\nu ^{\prime} ),\end{eqnarray} \tag{ 12 }$$

where the $\delta {{\rm{\Psi }}}_{{\rm{out}}}(\nu )$ are fluctuations in the output photon field about a steady-state solution. The derivation of $S(\nu )$ for the Dicke model and its incorporation into equation (9) is a central result of this work and is presented in appendix B of this paper. In general, $S(\nu )$ is a complex function with a series of poles displaced above the real ν axis and we employ a sync-transformation method to integrate equation (9). The approach can be applied to any model Hamiltonian system and provides the necessary connection between a microscopic model and its predicted photon coincidence.

In treating this as a scattering problem, we assume that the individual single photons impinging on the sample are uncorrelated with previous and subsequent photons. That is to say that the intensity of laser is small such that a single photon as left the cavity before the next entangled photon from the signal-arm of the HOM apparatus interacts with the sample. Under this approximation, we can treat the matter/photon interaction within a linear response approximation. A nonlinear theory must also include cross-correlation terms between the in- and out-going photon components and the material.

Before discussing the results of our calculations, it is important to recapitulate a number of aspects of the Dicke model and how these features are manifest in the photon coincidence counting rates. As stated already, we assume that the sample is in a steady state by exchanging the photons in the HOM apparatus with photons within the sample cavity and that ${ \mathcal S }(\nu )$ can be described within a linear-response theory. Because the cavity photons become entangled with the material excitations, the excitation frequencies are split into lower photonic (${\omega }_{-}$) and upper excitonic (${\omega }_{+}$) branches. Figure 2 shows the typical branching structure for the model for upper and lower branch polariton modes.

Zoom In Zoom Out Reset image size

At very low values of ${\lambda }_{k}$, the imaginary parts of the eigenvalues are equal and ${\omega }_{+}={\omega }_{-}$. In this over-damped regime, photons leak from the cavity before the photon/exciton state has undergone a single Rabi oscillation. At ${\lambda }_{k}=\kappa /2$ the system becomes critically damped and for ${\lambda }_{k}\lt \kappa /2$ and the degeneracy between the upper and lower polariton branches is lifted. As ${\lambda }_{k}$ increases above a critical value given by

$$\begin{eqnarray}&&{\lambda }_{c}=\sqrt{\displaystyle \frac{{\omega }_{k}{\omega }_{o}}{4}\left(1+\displaystyle \frac{{\kappa }^{2}}{{\omega }_{k}^{2}}\right)},\end{eqnarray} \tag{ 13 }$$

the system undergoes a quantum phase transition when ${\omega }_{-}=0$. Above this regime, excitations from the non-equilibrium steady state become collective and super-radiant. For our numerical results, unless otherwise noted we use dimensionless quantities, taking ${\omega }_{o}={\omega }_{k}=1.5$ for both the exciton frequency and cavity mode frequency, $\kappa =0.1$ for the cavity decay. These give a critical value of ${\lambda }_{c}=0.7517$. Also, within our model calculations, unless otherwise stated, we assume the central frequency of the signal and idler laser modes to be resonant with the cavity ${\omega }_{L}={\omega }_{k}={\omega }_{o}$.

We first consider the photon coincidence in the normal regime. Figures 3(a), (b) show the variation of the photon coincidence count when the laser frequency is resonant with the excitons (${\omega }_{k}={\omega }_{o}$). For low values of ${\lambda }_{k}$, the system is in the over-damped regime and the resulting coincidence scan reveals a slow decay for positive values of the time-delay. This is the perturbative regime or over-damped in which the scattering photon is dephased by the interaction with the sample, but there is insufficient time for the photon to become entangled with the sample. For ${\lambda }_{k}\gt \kappa /2$, the scattering photon is increasingly entangled with the material and further oscillatory structure begins to emerge in the coincidence scan.

Zoom In Zoom Out Reset image size

In the strong coupling regime, ${P}_{c}(t)$ becomes increasingly oscillatory with contributions from multiple frequency components. The origin of the structure is further revealed upon taking the Fourier cosine transform of ${P}_{c}(t)$ (equation (9)) taking the bandwidth of the bi-photon amplitude to be broad enough to span the full spectral range. The first two terms in the integral of equation (9) are independent of time and simply give a background count and can be ignored for purpose of analysis. The third term depends upon the time delay and is Fourier-cosine transform of the bi-photon amplitude times the scattering amplitudes,

$$\begin{eqnarray}&&{{ \mathcal P }}_{c}(\omega )=| { \mathcal F }(\omega ){| }^{2}{{ \mathcal S }}^{* }({\omega }_{L}-\omega ){ \mathcal S }({\omega }_{L}+\omega ).\end{eqnarray} \tag{ 14 }$$

As we show in the appendix B, ${ \mathcal S }(\omega )$ has a series of poles on the complex plane that correspond to the frequency spectrum of fluctuations about the matter-radiation steady state as given by the eigenvalues of ${{ \mathcal M }}_{s}$ in equation (B8). In figures 4(a), (c), (d) we show the evolution of pole-structure of ${{ \mathcal S }}^{* }({\omega }_{L}-z){ \mathcal S }({\omega }_{L}+z)$ (solid black) superimposed over the Fourier-cosine transform of the coincidence counts (${{ \mathcal P }}_{c}(\omega )$) with increasing coupling ${\lambda }_{k}$ for cases where the central laser frequency ${\omega }_{L}$ is resonant or off-resonant with the cavity (S). The location of poles can be readily understood in terms of eigenfrequencies of the polariton branches (see figure 2) with origins at $\pm ({\omega }_{L}+{\omega }_{o})$ and $\pm ({\omega }_{L}-{\omega }_{o})$.

Zoom In Zoom Out Reset image size

Both upper and lower polariton branches contribute to the scattering function and hence to the overall response indicating the entanglement between excitonic and photonic modes within the sample. However, the lion's share of the contribution to the response revealing that both excitonic $({\omega }_{+})$ and photonic $({\omega }_{-})$ branches contribute to the overall photon coincidence counting rates.

A closer examination of the pole structure in the vicinity of the phase transition reveals that two of the ${\omega }_{k}^{(-)}$ modes become degenerate over a small range of ${\lambda }_{k}$ but with different imaginary components indicating that the two modes decay at different rates. This is manifest in figure 4(b) by the rapid variation and divergence in the ${{ \mathcal S }}^{* }({\omega }_{L}-z){ \mathcal S }({\omega }_{L}+z)$ about ${\lambda }_{c}$ [29]. While the parametric width of this regime is small, it depends entirely upon the rate of photon exchange between the cavity and the laser field (κ).

We present here a formalism and method for connecting the photon coincidence signals for a sample placed in a HOM apparatus to the optical response of the coupled photon/material system. Our formalism reveals that by taking the Fourier transform of the ${P}_{c}(t)$ coincidence signal reveals the underlying pole structure of the entangled material/photon system. Our idea hinges upon an assumption that the interaction with the material preserves the initial entanglement between the two photons and that sample on the entanglement introduces an additional phase lag to one of the photons which we formally introduce in the form of a scattering response function ${ \mathcal S }$. The pole-structure in the output comes about from the further quantum entanglement of the signal photon with the sample.

Encoded in the time-delay signals is important information concerning the inner-workings of a quantum phase transition. Hence, we conclude that entangled photons with interferometric detection techniques provide a viable and tractable means to extract precise information concerning light–matter interactions. In particular, the approach reveals that at the onset of the symmetry-breaking transition between normal and super-radiant phases, two of the eigenmodes of the light–matter state exhibit distinctly different lifetimes. This signature of an intrinsic aspect of light–matter entanglement may be observed in a relatively simple experimental geometry with what amounts to a linear light-scattering/interferometry set up.

At first glance, it would appear that using quantum photons would not offer a clear advantage over more standard spectroscopies based upon a semi-classical description of the radiation field. However, the entanglement variable adds an additional dimension to the experiment allowing one to preform what would ordinarily be a nonlinear experiment using classical light as a linear experiment using quantized light. The recent works by Kalashnikov et al that inspired this work are perhaps the proverbial tip of the iceberg [15, 16].

The work at the University of Houston was funded in part by the National Science Foundation (CHE-1664971, MRI-1531814) and the Robert A Welch Foundation (E-1337). AP acknowledges the support provided by Los Alamos National Laboratory Directed Research and Development (LDRD) Funds. CS acknowledges support from the School of Chemistry & Biochemistry and the College of Science of Georgia Tech. JJ acknowledges the support of the Army Research Office (W911NF-13-1-0162). ARSK acknowledges funding from EU Horizon 2020 via Marie Sklodowska Curie Fellowship (Global) (Project No. 705874).

ERB, CS, and AP conceived and developed the core ideas. ERB, HL, and AP developed and implemented the formalism. ERB and JJ developed the contour integration method. All authors participated in the discussion of the methods and results. All authors participated in the drafting and final editing of this manuscript.

The authors declare no competing financial interests.

Our theoretical approach is to treat S as a material system interacting with a bath of quantum photons. We shall denote our 'system' as those degrees of freedom describing the material and the photons directly interacting with the sample, described by H_sys and assume that the photons within sample cavity are exchanged with external photons in the bi-photon field under the rotating-wave approximation,

$$\begin{eqnarray}&&{H}_{r}+{H}_{{rs}}={\rm{\hslash }}{\int }_{-\infty }^{\infty }\{{{zB}}_{k}^{\dagger }(z){B}_{k}(z)-{\rm{i}}\kappa (z)({\psi }_{k}^{\dagger }{B}_{k}(z)-{B}_{k}^{\dagger }(z){\psi }_{k})\}{\rm{d}}{z},\end{eqnarray} \tag{ A1 }$$

where $[{B}_{k}(z),{B}_{k^{\prime} }^{\dagger }(z^{\prime} )]={\delta }_{{kk}^{\prime} }\delta (z-z^{\prime} )$ are boson operators for photons in the laser field, and $[{\psi }_{k},{\psi }_{k}^{\dagger }]={\delta }_{{kk}^{\prime} }$ are boson operators for cavity photons in the sample that directly interact with the material component of the system.

The Heisenberg equations of motion for the reservoir and system photon modes are given by

$$\begin{eqnarray}&&{\partial }_{t}{B}_{k}(z)=-{{\rm{i}}{z}{B}}_{k}(z)+\kappa (z){\psi }_{k}\end{eqnarray} \tag{ A2 }$$

and

$$\begin{eqnarray}&&{\partial }_{t}{\psi }_{k}=-\displaystyle \frac{{\rm{i}}}{{\rm{\hslash }}}[{\psi }_{k},{H}_{{\rm{sys}}}]-\int \kappa (z){B}_{k}(z;t){\rm{d}}{z},\end{eqnarray} \tag{ A3 }$$

where the integration range is over all z. We can integrate formally the equations for the reservoir given either the initial or final states of the reservoir field

$$\begin{eqnarray}{B}_{k}(z;t)=\left\{\begin{array}{ll}{{\rm{e}}}^{-{\rm{i}}{z}(t-{t}_{i})}{B}_{k}(z;{t}_{i})+\kappa (z){\displaystyle \int }_{{t}_{i}}^{t}{{\rm{d}}{s}{\rm{e}}}^{-{\rm{i}}{z}(t-s)}{\psi }_{k}(s) & {\rm{for}}\,t\gt {t}_{i}\\ {{\rm{e}}}^{-{\rm{i}}{z}(t-{t}_{f})}{B}_{k}(z;{t}_{f})-\kappa (z){\displaystyle \int }_{t}^{{t}_{f}}{{\rm{d}}{s}{\rm{e}}}^{-{\rm{i}}{z}(t-s)}{\psi }_{k}(s) & {\rm{for}}\,t\lt {t}_{f}.\end{array}\right.\end{eqnarray} \tag{ A4 }$$

We shall eventually take ${t}_{i}\to -\infty$ and ${t}_{f}\to +\infty$ and require that the forward-time propagated and reverse-time propagated solutions are the same at some intermediate time t. If we assume that the coupling is constant over the frequency range of interest, we can write

$$\begin{eqnarray}&&\kappa (z)=\sqrt{\gamma /2\pi },\end{eqnarray} \tag{ A5 }$$

where γ is the rate that energy is exchanged between the reservoir and the system. This is the (first) Markov approximation.

Using these identities, one can find the Heisenberg equations for the cavity modes as

$$\begin{eqnarray}&&{\partial }_{t}{\psi }_{k}=-\displaystyle \frac{{\rm{i}}}{{\rm{\hslash }}}[{\psi }_{k},{H}_{{\rm{sys}}}]-\sqrt{\displaystyle \frac{\gamma }{2\pi }}{\int }_{-\infty }^{+\infty }{B}_{k}(z;t){\rm{d}}{z},\end{eqnarray} \tag{ A6 }$$

where H_sys is the Hamiltonian for the isolated system. We can now cast the external field in this equation in terms of its initial condition:

$$\begin{eqnarray}{\partial }_{t}{\psi }_{k} & = & -\displaystyle \frac{{\rm{i}}}{{\rm{\hslash }}}[{\psi }_{k},{H}_{{\rm{sys}}}]-\sqrt{\displaystyle \frac{\gamma }{2\pi }}{\displaystyle \int }_{-\infty }^{+\infty }{{\rm{e}}}^{{\rm{i}}{z}(t-{t}_{{\rm{i}}})}{B}_{k0}(z){\rm{d}}{z}\\ & & -\displaystyle \frac{\gamma }{2\pi }{\displaystyle \int }_{-\infty }^{+\infty }{\rm{d}}{z}{\displaystyle \int }_{{t}_{{\rm{i}}}}^{t}{{\rm{e}}}^{{\rm{i}}{z}(t-t^{\prime} )}{\psi }_{k}(t^{\prime} ){\rm{d}}{t}^{\prime} .\end{eqnarray} \tag{ A7 }$$

Let us define an input field in terms of the Fourier transform of the reservoir operators:

$$\begin{eqnarray}&&{\psi }_{k,{\rm{in}}}(t)=-\displaystyle \frac{1}{\sqrt{2\pi }}{\int }_{-\infty }^{+\infty }{{\rm{d}}{z}{\rm{e}}}^{-{\rm{i}}{z}(t-{t}_{i})}{B}_{{ko}}(z).\end{eqnarray} \tag{ A8 }$$

Since these depend upon the initial state of the reservoir, they are essentially a source of stochastic noise for the system. In our case, we shall use these as a formal means to connect the fields inside the sample to the fields in the laser cavity.

For the term involving $\gamma /2\pi$, the integral over frequency gives a delta-function:

$$\begin{eqnarray}&&{\int }_{-\infty }^{\infty }{{\rm{d}}{z}{\rm{e}}}^{{\rm{i}}{z}(t-t^{\prime} )}=2\pi \delta (t-t^{\prime} )\end{eqnarray} \tag{ A9 }$$

then

$$\begin{eqnarray}&&{\int }_{{t}_{o}}^{t}{\rm{d}}{t}^{\prime} \delta (t-t^{\prime} ){\psi }_{k}(t^{\prime} )=\displaystyle \frac{\psi (t)}{2}\,\mathrm{for}\,({t}_{o}\lt t\lt {t}_{f}).\end{eqnarray} \tag{ A10 }$$

This gives the forward equation of motion.

$$\begin{eqnarray}&&{\partial }_{t}{\psi }_{k}=-\displaystyle \frac{{\rm{i}}}{{\rm{\hslash }}}[{\psi }_{k},{H}_{{\rm{sys}}}]+\sqrt{\gamma }{\psi }_{k,{\rm{in}}}(t)-\displaystyle \frac{\gamma }{2}{\psi }_{k}(t).\end{eqnarray} \tag{ A11 }$$

We can also define an output field by integrating the reservoir backwards from time t_f to time t given a final state of the bath, B_kf.

$$\begin{eqnarray}&&{B}_{k}={{\rm{e}}}^{-{\rm{i}}{z}(t-{t}_{f})}{B}_{{kf}}-\sqrt{\displaystyle \frac{\gamma }{2\pi }}{\int }_{t}^{{t}_{f}}{{\rm{e}}}^{-{\rm{i}}{z}(t-t^{\prime} )}{\psi }_{k}(t^{\prime} ){\rm{d}}{t}^{\prime} .\end{eqnarray} \tag{ A12 }$$

This produces a similar equation of motion for the output field

$$\begin{eqnarray}&&{\partial }_{t}{\psi }_{k}=-\displaystyle \frac{{\rm{i}}}{{\rm{\hslash }}}[{\psi }_{k},{H}_{{\rm{sys}}}]-\sqrt{\gamma }{\psi }_{k,{\rm{out}}}(t)+\displaystyle \frac{\gamma }{2}{\psi }_{k}(t).\end{eqnarray} \tag{ A13 }$$

Upon integration:

$$\begin{eqnarray}&&{\psi }_{k,{\rm{out}}}(t)=\displaystyle \frac{1}{\sqrt{2\pi }}{\int }_{-\infty }^{+\infty }{\rm{d}}\nu {{\rm{e}}}^{-{\rm{i}}\nu (t-{t}_{f})}{B}_{{kf}}(\nu ).\end{eqnarray} \tag{ A14 }$$

At the time t, both equations must be the same, so we can subtract one from the other

$$\begin{eqnarray}&&{\psi }_{k,{\rm{in}}}+{\psi }_{k,{\rm{out}}}=\sqrt{\gamma }{\psi }_{k}\end{eqnarray} \tag{ A15 }$$

to produce a relation between the incoming and outgoing components. This eliminates the nonlinearity and explicit reference to the bath modes.

We now write ${\rm{\Psi }}=\{{\psi }_{k},{\psi }_{k}^{\dagger },{S}_{1},{S}_{2},\,\cdots \}$ as a vector of Heisenberg variables for the material system $\{{S}_{1},{S}_{2},\,\cdots \}$ and cavity modes $\{{\psi }_{k},{\psi }_{k}^{\dagger }\}$. Taking the equations of motion for the all fields to be linear and of the form

$$\begin{eqnarray}&&{\partial }_{t}{\rm{\Psi }}={{ \mathcal M }}_{{\rm{in}}}\cdot {\rm{\Psi }}+\sqrt{\gamma }{{\rm{\Psi }}}_{{\rm{in}}},\end{eqnarray} \tag{ A16 }$$

where ${{ \mathcal M }}_{{\rm{in}}}$ is a matrix of coefficients which are independent of time. The input vector ${{\rm{\Psi }}}_{{\rm{in}}}$ is non-zero for only the terms involving the input modes. We can also write a similar equation in terms of the output field; however, we have to account for the change in sign of the dissipation terms, so we denote the coefficient matrix as ${{ \mathcal M }}_{{\rm{out}}}$. In this linearized form, the forward and reverse equations of motion can be solved formally using the Laplace transform, giving

$$\begin{eqnarray}&&({{ \mathcal M }}_{{\rm{in}}}-{\rm{i}}{z}){\rm{\Psi }}(z)=-\sqrt{\gamma }{{\rm{\Psi }}}_{{\rm{in}}}(z),\end{eqnarray} \tag{ A17 }$$

$$\begin{eqnarray}&&({{ \mathcal M }}_{{\rm{out}}}-{\rm{i}}{z}){\rm{\Psi }}(z)=+\sqrt{\gamma }{{\rm{\Psi }}}_{{\rm{out}}}(z).\end{eqnarray} \tag{ A18 }$$

These and the relation ${{\rm{\Psi }}}_{{\rm{in}}}+{{\rm{\Psi }}}_{{\rm{out}}}=\sqrt{\gamma }{\rm{\Psi }}$ allows one to eliminate the external variables entirely:

$$\begin{eqnarray}&&{{\rm{\Psi }}}_{{\rm{out}}}(z)=-({{ \mathcal M }}_{{\rm{out}}}-{\rm{i}}{z}){({{ \mathcal M }}_{{\rm{in}}}-{\rm{i}}{z})}^{-1}{{\rm{\Psi }}}_{{\rm{in}}}(z).\end{eqnarray} \tag{ A19 }$$

This gives a precise connection between the input and output fields. More over, the final expression does not depend upon the assumed exchange rate between the internal ${\psi }_{k}$ and external ${B}_{k}(z)$ photon fields. The procedure is very much akin to the use of Møller operators in scattering theory. To explore this connection, define ${\hat{{\rm{\Omega }}}}_{{\rm{in}}}^{(\pm )}$ as an operator which propagates an incoming solution at $\mp \infty$ to the interaction at time t = 0 and its reverse ${\hat{{\rm{\Omega }}}}_{{\rm{out}}}^{(\pm )}$ which propagates an out-going solution at $\mp \infty$ back to the interaction at time t = 0.

$$\begin{eqnarray}&&{\hat{{\rm{\Omega }}}}_{{\rm{in}},{\rm{out}}}^{(\pm )\dagger }{\hat{{\rm{\Omega }}}}_{{\rm{in}},{\rm{out}}}^{(\pm )}=I\end{eqnarray} \tag{ A20 }$$

and

$$\begin{eqnarray}&&{\hat{{\rm{\Omega }}}}_{{\rm{in}}}^{(\pm )}={({{ \mathcal M }}_{{\rm{in}}}\mp {\rm{i}}{z})}^{-1},\end{eqnarray} \tag{ A21 }$$

$$\begin{eqnarray}&&{\hat{{\rm{\Omega }}}}_{{\rm{out}}}^{(\pm )}={({{ \mathcal M }}_{{\rm{out}}}\pm {\rm{i}}{z})}^{-1}.\end{eqnarray} \tag{ A22 }$$

Thus, we can write equation (A19) as

$$\begin{eqnarray}&&{{\rm{\Psi }}}_{{\rm{out}}}(z)=-{\hat{{\rm{\Omega }}}}_{{\rm{out}}}^{(-)\dagger }{\hat{{\rm{\Omega }}}}_{{\rm{in}}}^{(+)}{{\rm{\Psi }}}_{{\rm{in}}}(z).\end{eqnarray} \tag{ A23 }$$

To compute the response function, we consider fluctuations and excitations from a steady-state solution:

$$\begin{eqnarray}&&{\rm{\Psi }}(t)={{\rm{\Psi }}}_{{ss}}+\delta {\rm{\Psi }}(t).\end{eqnarray} \tag{ A24 }$$

The resulting linearized equations of motion read

$$\begin{eqnarray}&&\displaystyle \frac{{\rm{d}}}{{\rm{d}}{t}}\delta {\rm{\Psi }}(t)={{ \mathcal M }}_{s}\delta {\rm{\Psi }}(t)\end{eqnarray} \tag{ A25 }$$

implying a formal solution of

$$\begin{eqnarray}&&\delta {\rm{\Psi }}(t)={{\rm{e}}}^{{{ \mathcal M }}_{s}t}\delta {\rm{\Psi }}(0).\end{eqnarray} \tag{ A26 }$$

From this we deduce that the eigenvalues and eigenvectors of ${{ \mathcal M }}_{s}$ give the fluctuations in terms of the normal excitations about the stationary solution. Using the input/output formalism, we can write the outgoing state (in terms of the Heisenberg variables) in terms of their input values:

$$\begin{eqnarray}&&\delta {{\rm{\Psi }}}_{{\rm{out}}}(z)=-({{ \mathcal M }}_{{\rm{out}},s}-{\rm{i}}{z}{I}){({{ \mathcal M }}_{{\rm{in}},s}-{\rm{i}}{z}{I})}^{-1}\delta {{\rm{\Psi }}}_{{\rm{in}}}(z),\end{eqnarray} \tag{ A27 }$$

where as given above, $\delta {\rm{\Psi }}(z)$ is a vector containing the fluctuations about the stationary values for each of the Heisenberg variables. The ${{ \mathcal M }}_{{\rm{in}},s}$ and ${{ \mathcal M }}_{{\rm{out}},s}$ are the coefficient matrices from the linearisation process. The input field satisfies $\langle \delta {\psi }_{{\rm{in}}}(z)\delta {\psi }_{{\rm{in}}}^{\dagger }(z^{\prime} )\rangle =\delta (z-z^{\prime} )$ and all other terms are zero. Thus, the transmission function is given by

$$\begin{eqnarray}&&\delta (z-z^{\prime} ){ \mathcal S }(z)=\langle \delta {\psi }_{{\rm{out}}}^{\dagger }(z)\delta {\psi }_{{\rm{out}}}(z^{\prime} )\rangle .\end{eqnarray} \tag{ A28 }$$

In other words, the ${ \mathcal S }(z)$ is the response of the system to the input field of the incoming photon state producing an output field for the out-going photon state.

Let us consider an ensemble of N identical two-level systems corresponding to local molecular sites coupled to a set of photon modes described by ${\psi }_{k}$.

$$\begin{eqnarray}\hat{H} & = & \displaystyle \sum _{j}\displaystyle \frac{{\rm{\hslash }}{\omega }_{o}}{2}{\hat{\sigma }}_{z,j}+\displaystyle \sum _{k}{\rm{\hslash }}{\omega }_{k}{\hat{\psi }}_{k}^{\dagger }{\hat{\psi }}_{k}+\displaystyle \sum _{k,j}\displaystyle \frac{{\rm{\hslash }}{\lambda }_{{kj}}}{\sqrt{N}}({\hat{\psi }}_{k}^{\dagger }+{\hat{\psi }}_{k})({\hat{\sigma }}_{j}^{+}+{\hat{\sigma }}_{j}^{-}),\end{eqnarray} \tag{ B1 }$$

where $\{{\hat{\sigma }}_{z,j},{\hat{\sigma }}_{j}^{\pm }\}$ are local spin-1/2 operators for site j, ${\rm{\hslash }}{\omega }_{o}$ is the local excitation energy, and ${\lambda }_{{kj}}$ is the coupling between the kth photon mode and the jth site, which we will take to be uniform over all sites. Defining the total angular momentum operators

$$\begin{eqnarray*}&&{\hat{J}}_{z}=\displaystyle \sum _{j}{\hat{\sigma }}_{z,j}\,\,\mathrm{and}\,\,{\hat{J}}_{\pm }=\displaystyle \sum _{j}{\hat{\sigma }}_{j}^{\pm }\end{eqnarray*}$$

and

$$\begin{eqnarray*}&&{\hat{J}}^{2}={\hat{J}}_{z}^{2}+({\hat{J}}_{+}{\hat{J}}_{-}+{\hat{J}}_{-}{\hat{J}}_{+})/2\end{eqnarray*}$$

as the total angular momentum operator, this Hamiltonian can be cast in the form in equation (B1) by mapping the total state space of N spin 1/2 states onto a single angular momentum state vector $| J,M\rangle$.

$$\begin{eqnarray}\hat{H} & = & {\rm{\hslash }}{\omega }_{o}{\hat{J}}_{z}+\displaystyle \sum _{k}{\rm{\hslash }}{\omega }_{k}{\hat{\psi }}_{k}^{\dagger }{\hat{\psi }}_{k}+\displaystyle \sum _{k}\displaystyle \frac{{\rm{\hslash }}{\lambda }_{k}}{\sqrt{N}}({\hat{\psi }}_{k}^{\dagger }+{\hat{\psi }}_{k})({\hat{J}}_{+}+{\hat{J}}_{-}).\end{eqnarray} \tag{ B2 }$$

Note that the ground state of the system corresponds to $| J,-J\rangle$ in which each molecule is in its electronic ground state. Excitations from this state create up to N excitons within the system corresponding to the state $| J,+J\rangle$. Intermediate to this are multi-exciton states which correspond to various coherent superpositions of local exciton configurations. For each value of the wave vector k one obtains the following Heisenberg equations of motion for the expectation values of the operators

$$\begin{eqnarray}&&\displaystyle \frac{\partial {\hat{\psi }}_{k}}{\partial t}=(-{\rm{i}}{\omega }_{k}-\kappa ){\hat{\psi }}_{k}-{\rm{i}}\displaystyle \frac{{\lambda }_{k}}{\sqrt{N}}({\hat{J}}_{+}+{\hat{J}}_{-}),\end{eqnarray} \tag{ B3 }$$

$$\begin{eqnarray}&&\displaystyle \frac{\partial {\hat{\psi }}_{k}^{\dagger }}{\partial t}=({\rm{i}}{\omega }_{k}-\kappa ){\hat{\psi }}_{k}^{\dagger }+{\rm{i}}\displaystyle \frac{{\lambda }_{k}}{\sqrt{N}}({\hat{J}}_{+}+{\hat{J}}_{-}),\end{eqnarray} \tag{ B4 }$$

$$\begin{eqnarray}&&\displaystyle \frac{\partial {\hat{J}}_{\pm }}{\partial t}=\pm {\rm{i}}{\omega }_{o}{\hat{J}}_{\pm }\mp 2{\rm{i}}{\hat{J}}_{z}\displaystyle \sum _{k}\displaystyle \frac{{\lambda }_{k}}{\sqrt{N}}({\hat{\psi }}_{k}+{\hat{\psi }}_{k}^{\dagger }),\end{eqnarray} \tag{ B5 }$$

$$\begin{eqnarray}&&\displaystyle \frac{\partial {\hat{J}}_{z}}{\partial t}={\rm{i}}\displaystyle \frac{{\lambda }_{k}}{\sqrt{N}}({\hat{J}}_{-}-{\hat{J}}_{+})({\hat{\psi }}_{k}+{\hat{\psi }}_{k}^{\dagger }),\end{eqnarray} \tag{ B6 }$$

where κ is the decay of photon ${\hat{\psi }}_{k}$ into the reservoir. These are nonlinear equations and we shall seek stationary solutions and linearize about them. The linearized equation for the fluctuations reads

$$\begin{eqnarray}&&\displaystyle \frac{{\rm{d}}}{{\rm{d}}{t}}\delta {\rm{\Psi }}(t)={{ \mathcal M }}_{s}\delta {\rm{\Psi }}(t)\end{eqnarray} \tag{ B7 }$$

with

$$\begin{eqnarray}{{ \mathcal M }}_{s}=\left[\begin{array}{ccccc}-(\kappa -{\rm{i}}{\omega }_{k}) & 0 & {\rm{i}}{\lambda }_{k} & {\rm{i}}{\lambda }_{k} & 0\\ 0 & -(\kappa +{\rm{i}}{\omega }_{k}) & -{\rm{i}}{\lambda }_{k} & -{\rm{i}}{\lambda }_{k} & 0\\ 2{\rm{i}}{\lambda }_{k}{\overline{J}}_{z} & 2{\rm{i}}{\lambda }_{k}{\overline{J}}_{z} & -{\rm{i}}{\omega }_{o} & 0 & 2{\rm{i}}{\lambda }_{k}({\overline{\psi }}_{k}+{\overline{\psi }}_{k}^{\dagger })\\ -2{\rm{i}}{\lambda }_{k}{\overline{J}}_{z} & -2{\rm{i}}{\lambda }_{k}{\overline{J}}_{z} & 0 & {\rm{i}}{\omega }_{o} & -2{\rm{i}}{\lambda }_{k}({\overline{\psi }}_{k}+{\overline{\psi }}_{k}^{\dagger })\\ {\rm{i}}{\lambda }_{k}({\overline{J}}_{-}-{\overline{J}}_{+}) & {\rm{i}}{\lambda }_{k}({\overline{J}}_{-}-{\overline{J}}_{+}) & {\rm{i}}{\lambda }_{k}({\overline{\psi }}_{k}+{\overline{\psi }}_{k}) & -{\rm{i}}{\lambda }_{k}({\overline{\psi }}_{k}+{\overline{\psi }}_{k}) & 0\end{array}\right],\end{eqnarray} \tag{ B8 }$$

where ${\overline{J}}_{z,\pm }$, ${\overline{\psi }}_{k}$, and ${\overline{\psi }}_{k}^{\dagger }$ denote the steady-state solutions. Where we have removed N from the equations of motion by simply rescaling the variables. The cavity photon decay $\kappa =\pm \gamma /2$ for the input and output equations of motion, respectively. The model has both trivial and non-trivial stationary solutions corresponding to the normal and super-radiant regimes. For the normal regime,

$$\begin{eqnarray}&&{\overline{\psi }}_{k}={\overline{\psi }}_{k}^{\dagger }={\overline{J}}_{\pm ,s}=0\end{eqnarray} \tag{ B9 }$$

and

$$\begin{eqnarray}&&{\overline{J}}_{z}=\pm \displaystyle \frac{N}{2},\end{eqnarray} \tag{ B10 }$$

which correspond to the case where every spin is excited or in the ground state. Since we are primarily interested in excitations from the electronic ground state, we initially focus our attention to these solutions.

Non-trivial solutions to these equations predict that above a critical value of the coupling $\lambda \gt {\lambda }_{c}$, the system will undergo a quantum phase transition to form a super-radiant state. It should be pointed out that in the original Dicke model, above the critical coupling, the system is no longer gauge invariant leading to a violation of the Thomas–Reiche–Kuhn (TRK) sum rule. Gauge invariance can be restored; however, the system no longer undergoes a quantum phase transition [30]. However, for a driven, non-equilibrium system such as presented here, the TRK sum rule does not apply and the quantum phase transition is a physical effect.

The non-trivial solutions for the critical regime are given by

$$\begin{eqnarray}&&{\overline{\psi }}_{k}=\mp \displaystyle \frac{{\lambda }_{k}}{{\omega }_{k}-{\rm{i}}\kappa }\sqrt{1-{\left(\displaystyle \frac{{\lambda }_{c}}{{\lambda }_{k}}\right)}^{4}},\end{eqnarray} \tag{ B11 }$$

$$\begin{eqnarray}&&{\overline{J}}_{-}=\pm \displaystyle \frac{1}{2}{\left(1-{\left(\displaystyle \frac{{\lambda }_{c}}{{\lambda }_{k}}\right)}^{4}\right)}^{1/2},\end{eqnarray} \tag{ B12 }$$

$$\begin{eqnarray}&&{\overline{J}}_{z}=-\displaystyle \frac{1}{2}{\left(\displaystyle \frac{{\lambda }_{c}}{{\lambda }_{k}}\right)}^{2}.\end{eqnarray} \tag{ B13 }$$

The (imaginary) eigenvalues of ${{ \mathcal M }}_{s}$ gives 4 non-zero and 1 trivial normal mode frequencies (for $\kappa =0$), which we shall denote as

$$\begin{eqnarray}&&{\omega }_{\pm }^{2}=-\displaystyle \frac{1}{2}({\omega }_{k}^{2}+{\omega }_{o}^{2}\pm \sqrt{{({\omega }_{k}^{2}-{\omega }_{o}^{2})}^{2}+16{\lambda }_{k}^{2}{\omega }_{k}{\omega }_{o}}).\end{eqnarray} \tag{ B14 }$$

One also obtains the critical coupling constant

$$\begin{eqnarray}&&{\lambda }_{c}=\sqrt{\displaystyle \frac{{\omega }_{k}{\omega }_{o}}{4}\left(1+\displaystyle \frac{{\kappa }^{2}}{{\omega }_{k}^{2}}\right)}.\end{eqnarray} \tag{ B15 }$$

Figure 2 gives the normal mode spectrum for a resonant system with ${\omega }_{k}={\omega }_{o}=1.5$ and $\kappa =0.1$ (in reduced units).

The integral for the photon coincidence can be problematic to evaluate numerically given the oscillatory nature of the sinc function in ${ \mathcal F }(z)$. To accomplish this, we define a sinc-transformation based upon ${ \mathcal F }(z)$ using the identity

$$\begin{eqnarray}&&\mathrm{sinc}(z)=\displaystyle \frac{\sin (z)}{z}=\displaystyle \frac{1}{2}{\int }_{-1}^{1}{{\rm{e}}}^{{\rm{i}}{k}{z}}{\rm{d}}{k}\end{eqnarray} \tag{ C1 }$$

which yields

$$\begin{eqnarray}&&{ \mathcal F }(z)=\displaystyle \frac{1}{2}{\int }_{-1}^{1}{{\rm{e}}}^{{{\rm{i}}{k}{b}{z}}^{2}}\,{\rm{d}}{k}.\end{eqnarray} \tag{ C2 }$$

From this we can re-write each term in equations (9) in the form

$$\begin{eqnarray}&&I(t)={\int }_{-\infty }^{\infty }{\rm{d}}{z}| { \mathcal F }(z){| }^{2}{{ \mathcal S }}^{* }({\omega }_{L}-z){ \mathcal S }({\omega }_{L}+z){{\rm{e}}}^{2{\rm{i}}{z}{t}},\end{eqnarray} \tag{ C3 }$$

$$\begin{eqnarray}&&=\,\displaystyle \frac{1}{4}{\int }_{-1}^{1}{\rm{d}}{k}{\int }_{-1}^{1}{\rm{d}}{k}^{\prime} { \mathcal G }(k-k^{\prime} ,t),\end{eqnarray} \tag{ C4 }$$

whereby we denote

$$\begin{eqnarray}&&{ \mathcal G }(q,t)={\int }_{-\infty }^{\infty }{\rm{d}}{z}\,{{\rm{e}}}^{-{{\rm{i}}{b}{z}}^{2}q+2{\rm{i}}{t}{z}}{{ \mathcal S }}^{* }({\omega }_{L}-z){ \mathcal S }({\omega }_{L}+z).\end{eqnarray} \tag{ C5 }$$

The integrand is highly oscillatory along the z-axis; however, for non-zero $k-k^{\prime} =q$, ${ \mathcal G }$ becomes a Gaussian integral under coordinate transformation obtained by completing the square:

$$\begin{eqnarray}&&-{{\rm{i}}{b}{q}{z}}^{2}+2{\rm{i}}{t}{z}=-{\rm{i}}{b}{q}\left({z}^{2}-2\displaystyle \frac{t}{{bq}}z\right)\end{eqnarray} \tag{ C6 }$$

$$\begin{eqnarray}&&=-{\rm{i}}{b}{q}\left[{\left(z-\displaystyle \frac{t}{{bq}}\right)}^{2}-{\left(\displaystyle \frac{t}{{bq}}\right)}^{2}\right]\end{eqnarray} \tag{ C7 }$$

$$\begin{eqnarray}&&=\,-{\rm{i}}{b}{q}{\left(z-\displaystyle \frac{t}{{bq}}\right)}^{2}-\displaystyle \frac{{t}^{2}}{{\rm{i}}{b}{q}}.\end{eqnarray} \tag{ C8 }$$

For $q\gt 0$, we take $u=(\sqrt{{\rm{i}}})\left(z-\tfrac{t}{{bq}}\right)$, and for $q\lt 0$, we take $u=(-{\rm{i}}\sqrt{{\rm{i}}})\left(z-\tfrac{t}{{bq}}\right)$. Solving for z yields: $z=(-{\rm{i}}\sqrt{{\rm{i}}})u+\tfrac{t}{{bq}}$ and $z=(\sqrt{{\rm{i}}})u+\tfrac{t}{{bq}}$, respectively. In short, the optimal contour of the Gaussian integral is obtained by rotating by $\pi /4$ from the real-axis in the anticlockwise direction for the case of $q\gt 0$ and by $\pi /4$ in the clockwise direction for the case of $q\lt 0$ as indicated in figure 5.

Zoom In Zoom Out Reset image size

The spectral response ${{ \mathcal S }}^{* }({\omega }_{L}-z){ \mathcal S }({\omega }_{L}+z)$ has a number of poles on the complex z plane above the real-z axis. We now use the residue theorem to evaluate the necessary poles which result as the contour rotates from the real axis to the complex $\pm \pi /4$ axis. The 8 second order poles, $\{{\rho }_{n}\}$, are defined by roots of the denominators $D({\omega }_{L}-z)$ and $D({\omega }_{L}+z)$ and located $\kappa /2$ above the real axis at locations symmetrically placed around the origin. For the counter-clockwise rotation ($q\gt 0$), poles included to the right of the real crossing point will be added and those to the left will be ignored (${{ \mathcal P }}_{R}(t)$) ; whereas for a clockwise rotation ($q\lt 0$), the left-hand poles will be subtracted and the right-hand poles will be ignored (${{ \mathcal P }}_{L}(t)$).

For the unique case q = 0, all poles are summed (${{ \mathcal P }}_{{\rm{all}}}(t)$).

$$\begin{eqnarray}&&{ \mathcal P }(t)=2\pi {\rm{i}}\displaystyle \sum _{n}\mathop{\mathrm{lim}}\limits_{z\to {\rho }_{n}}{}_{}\displaystyle \frac{{\rm{d}}}{{\rm{d}}{z}}[{(z-{\rho }_{n})}^{2}{{ \mathcal S }}^{* }({\omega }_{L}-z){ \mathcal S }({\omega }_{L}+z){{\rm{e}}}^{-{{\rm{i}}{b}{z}}^{2}q+2{\rm{i}}{z}{t}}].\end{eqnarray} \tag{ C9 }$$

Thus, we obtain

$$\begin{eqnarray}{ \mathcal G }(q\gt 0,t) & = & {{\rm{e}}}^{\tfrac{{{\rm{i}}{t}}^{2}}{{bq}}}{\displaystyle \int }_{-\infty }^{\infty }{{\rm{d}}{u}{\rm{e}}}^{-{{bqu}}^{2}}{{ \mathcal S }}^{* }\left({\omega }_{L}-(-{\rm{i}}\sqrt{{\rm{i}}})u-\displaystyle \frac{t}{{bq}}\right){ \mathcal S }\left({\omega }_{L}+(-{\rm{i}}\sqrt{{\rm{i}}})u+\displaystyle \frac{t}{{bq}}\right)\\ & & +{{ \mathcal P }}_{R}(t),\end{eqnarray} \tag{ C10 }$$

$$\begin{eqnarray}&&{ \mathcal G }(q=0,t\gt 0)={{ \mathcal P }}_{{\rm{all}}}(t),\end{eqnarray} \tag{ C11 }$$

$$\begin{eqnarray}&&{ \mathcal G }(q=0,t\lt 0)=0,\end{eqnarray} \tag{ C12 }$$

$$\begin{eqnarray}{ \mathcal G }(q\lt 0,t) & = & {{\rm{e}}}^{\tfrac{{{\rm{i}}{t}}^{2}}{{bq}}}{\displaystyle \int }_{-\infty }^{\infty }{{\rm{d}}{u}{\rm{e}}}^{+{{bqu}}^{2}}{{ \mathcal S }}^{* }\left({\omega }_{L}-(\sqrt{{\rm{i}}})u-\displaystyle \frac{t}{{bq}}\right){ \mathcal S }\left({\omega }_{L}+(\sqrt{{\rm{i}}})u+\displaystyle \frac{t}{{bq}}\right)\\ & & -{{ \mathcal P }}_{L}(t),\end{eqnarray} \tag{ C13 }$$

whereby ${{ \mathcal P }}_{L}$, ${{ \mathcal P }}_{L}$, and ${{ \mathcal P }}_{{\rm{all}}}$ are summed over the left, right, or all poles, respectively. The resulting expressions are analytic (albeit lengthy) and defined by exponentials and low order polynomials. Completion of the k and $k^{\prime}$ integrals yield an exact expression for the response.

In the non-critical regime, we find the response ${ \mathcal S }(z)$ in a closed form as

$$\begin{eqnarray}&&{ \mathcal S }(z)=-\displaystyle \frac{\pi }{D{(z)}^{2}}({S}_{o}+{S}_{2}{z}^{2}+{S}_{4}{z}^{4}+{S}_{6}{z}^{6}+{S}_{8}{z}^{8}),\end{eqnarray} \tag{ D1 }$$

where

$$\begin{eqnarray}&&{S}_{0}={\omega }_{o}^{2}({\kappa }^{4}{\omega }_{o}^{2}+2{\kappa }^{2}(16{\lambda }_{k}^{4}-4{\lambda }_{k}^{2}{\omega }_{k}{\omega }_{o}+{\omega }_{k}^{2}{\omega }_{o}^{2})+{\omega }_{k}^{2}{({\omega }_{k}{\omega }_{o}-4{\lambda }_{k}^{2})}^{2}),\end{eqnarray} \tag{ D2 }$$

$$\begin{eqnarray}&&{S}_{2}=-2{\omega }_{o}({\kappa }^{4}{\omega }_{o}-{\kappa }^{2}(4{\lambda }_{k}^{2}{\omega }_{k}-2{\omega }_{k}^{2}{\omega }_{o}+{\omega }_{o}^{3})+{\omega }_{k}({\omega }_{k}^{2}+{\omega }_{o}^{2})({\omega }_{k}{\omega }_{o}-4{\lambda }_{k}^{2})),\end{eqnarray} \tag{ D3 }$$

$$\begin{eqnarray}&&{S}_{4}={\kappa }^{4}+2{\kappa }^{2}({\omega }^{2}-2{\omega }_{o}^{2})-8{\lambda }^{2}{\omega }_{k}{\omega }_{o}+{\omega }_{k}^{4}+4{\omega }_{k}^{2}{\omega }_{o}^{2}+{\omega }_{o}^{4},\end{eqnarray} \tag{ D4 }$$

$$\begin{eqnarray}&&{S}_{6}=2({\kappa }^{2}-{\omega }_{k}^{2}-{\omega }_{o}^{2}),\end{eqnarray} \tag{ D5 }$$

$$\begin{eqnarray}&&{S}_{8}=1,\end{eqnarray} \tag{ D6 }$$

$$\begin{eqnarray}&&D(z)=[({z}^{2}-{\omega }_{o}^{2})({\omega }_{k}^{2}+{(\kappa +{\rm{i}}{z})}^{2})+4{\lambda }_{k}^{2}{\omega }_{k}{\omega }_{o}].\end{eqnarray} \tag{ D7 }$$

The four complex-valued roots of $D(z)=0$ correspond to the eigenvalues of the ${{ \mathcal M }}_{s}$ matrix in equation (B8). It is possible using Mathematica to obtain a closed-form expression in the super-radiant regime; however, the resulting expression is very lengthy and we do not reproduce it here.