Single-photon emission mediated by single-electron tunneling in plasmonic nanojunctions

Recent scanning tunnelling microscopy (STM) experiments reported single-molecule fluorescence induced by tunneling currents in the nanoplasmonic cavity formed by the STM tip and the substrate.The electric field of the cavity mode couples with the current-induced charge fluctuations of the molecule, allowing the excitation of the mode. We investigate theoretically this system for the experimentally relevant limit of large damping rate $\kappa$ for the cavity mode and arbitrary coupling strength to a single-electronic level. We find that for bias voltages close to the first inelastic threshold of photon emission, the emitted light displays anti-bunching behavior with vanishing second-order photon correlation function. At the same time, the current and the intensity of emitted light display Franck--Condon steps at multiples of the cavity frequency $\omega_c$ with a width controlled by $\kappa$ rather than the temperature $T$. For large bias voltages, we predict strong photon bunching of the order of the $\kappa/\Gamma$ where $\Gamma$ is the electronic tunneling rate. Our theory thus predicts that strong coupling to a single level allows current-driven non-classical light emission.

Electronic transport coupled to the field of an electromagnetic cavity can be realized in a wealth of different systems. This includes in the microwave range carbonnanotubes [1][2][3][4][5], quantum-dots [6][7][8][9], and Josephsonjunctions [10][11][12], or in the optical range, molecules in plasmonic nanocavities, formed by an STM tip with a substrate [13][14][15][16][17][18][19][20][21][22][23][24] and organic microcavities [25][26][27], or with waveguide quantum electrodynamic systems [28][29][30]. The reduction of the cavity volume V allows to increase the zero-point quantum fluctuations of the electric field E zpm ∼ V −1/2 . This motivated optical studies of molecular two-level systems strongly coupled to the cavity field by the dipolar interaction Λ d ∼ pE zpm , (with p the molecule dipole moment). One of the goals of this effort is to reach Λ d larger than κ, which has been and remains challenging, despite recent achievements [31]. On the other side, the coupling of a cavity mode to the current-induced charge fluctuations of a single-electronic level is given by a monopolar coupling constant Λ m ∼ eLE zpm as derived in Ref. [32] (see also [33]), with L the typical extension of the transport region and e the electronic charge. Since typically in a given system eL p, the monopolar coupling constant is much larger than the dipolar one [32]. This probably contributed to the observation of values of Λ m larger than κ in microwave cavities coupled to electronic transport [4,6,8] and even approaching the cavity resonating frequency ω c ( = 1) [11,34,35]. Recent results in plasmonic cavities coupled to electronic transport [17,21,22] open thus the possibility to explore transport through a single electronic level in these structures. This is expected to reach much larger coupling constants than those currently observed for purely dipolar coupling, requiring further theoretical investigations.
The system presents strong analogies with electron-transport coupled to molecular vibrations. This has been investigated in different regimes, leading to the striking prediction of Franck-Condon blockade [36][37][38] and its observation [39,40]. However, there are important differences: The first is the low quality factor of plasmonic cavities, which is typically of the order of 10 [31]. The second, and more interesting, is that the state of the optical or microwave cavity can be directly measured by detecting the emitted photons. It is thus important to investigate how transport through a molecule is linked to the property of the emitted radiation [41][42][43].
In this paper we consider electronic transport through a single-level quantum dot, where the charge on the dot is coupled to the electric field of an electromagnetic cavity. We propose a theoretical model to obtain the cur- Schematic of two metallic electrodes forming a plasmonic nanocavity characterized by a resonating frequency ωc/2π and damping rate κ. A single electronic level ε0 of a molecule in the nanogap couples to the electromagnetic radiation with coupling constant Λm. Electrons can tunnel to and from the dot with tunneling rates Γα. Voltage drops, Vα, with respect to ε0 are indicated. rent through the quantum dot taking into account the cavity dissipation κ, and arbitrary values of the coupling strength in the incoherent transport regime Γ k B T , with Γ the electron tunneling rate, T the temperature, and k B the Boltzmann constant. Similarly to the Franck-Condon case, we find current steps at the inelastic thresholds for photon emission, but with a width controlled by κ rather than T . We also derive the photon distribution and the second-order photon correlation function g (2) (t), where t is the emission time. Its behavior for t = 0 clearly shows that close to the first threshold for photon emission, for Γ/κ 1, and Λ m ≈ ω c , photons anti-bunch: The junction becomes a single-photon source based on single-electron tunneling. This mechanism is different from the one assumed to be responsible for recently observed anti-bunching of emitted light in STM plasmonic nanojunctions involving multiple electronic levels [20]. For large bias voltages we find instead strong photon bunching with g (2) (0) ≈ κ/Γ 1. Model. Figure 1 shows a schematic of the system at hand. The Hamiltonian is written with d † the creation operator for the electron on the dot single level of energyε 0 , a † the creation operator for the photon field, and Λ m = λω c the coupling constant. We follow Ref. [32] for the derivation of the interaction term [33]. We neglect direct coupling of the cavity field to the electrons in the leads, since this effect is analogous to the coupling of the cavity field to the photon bath. We treat the electrons in the leads and the propagating electromagnetic modes as a bath: where c † αk and b † q are the creation operators for the electrons on the leads α = L, R with energy ε αk and for the propagating photons of energy ω q , respectively. The (linear) coupling to the bath is given by with t αk and l q the tunneling amplitudes. We first perform a standard Lang-Firsov unitary transformation on the HamiltonianH = U HU † , with U = e λd † d(a−a † ) . This removes explicitly the electron-photon coupling term in H S , shifts the dot-level energy ε 0 =ε 0 −Λ 2 m /ω c and modifies the d operator in Master-equation. Let us define the reduced density matrix ρ(t) for the d and a degrees of freedom after tracing out the bath. We assume that the molecule is sufficiently isolated from the substrate, as is reasonable for STM experiments performed on thin insulating films [14,[19][20][21][22]. We consider then the relevant regime Γ k B T ω c where the dynamics of ρ(t) can be described by the Born-Markov master equation: The first term contributing to the Liouvilian operator L gives the coherent evolution of ρ(t). The sec-  ond one describes the damping of the cavity mode [44,45]: Bose distribution of photons in the bath at the cavity frequency. The last term describes incoherent electron tunneling L e ρ(t) = D − ρ(t) − ρ(t)D + , D † + h.c. [46], where , and Γ α = 2π k |t αk | 2 δ(ω − ε αk ) is the tunneling rate from the lead α, that in the usual wide-band approximation becomes ω-independent. Finally, D I (−t) is the D operator in the interaction representation with respect toH S . We introduced the short-hand notation  Electronic current. Using the previous results, we derive the expression for the average electronic dc-current evaluated at lead α where we introduced the Fourier transform of f ± α , S DD † and S D † D . Equation (3) enables to calculate the current in presence of strong damping rates κ, that for plasmonic cavities reaches low quality factors ω c /κ ≈ 10 [31]. It allows to include the damping of the electromagnetic field during the tunneling process. This expression and ρ st can be evaluated numerically by projecting on the charge and harmonic oscillator basis. Figure 2(a) reports the electronic current I for strong coupling, λ = 1.4, as a function of the relative voltage drops eV α = µ α − ε 0 between the chemical potential of lead α and the dot energy level [49]. Specific current-voltage characteristics, corresponding to symmetric (V L = −V R ) and asymmetric (eV R = −0.2ω c ) voltage drops, are shown respectively as full and dashed lines in Fig. 2(b), for weak (λ = 0.1), moderate (λ = 0.6), and strong (λ = 1.4) coupling strengths. These exhibit similar features of the Franck-Condon blockade regime [36][37][38], with inelastic steps observed each time the voltage drop eV L = nω c matches a multiple of the cavity-photon frequency. This is the threshold for oneelectron tunneling while emitting n photons in the cavity. The step heights are given by the Poisson distribution P n (λ) = e −λ 2 λ 2n /n! and the width of the inelastic steps by the cavity-losses κ, which exceed the temper-ature broadening. This is analogous to the broadening of phonon sidebands by frictional damping [36], but our treatment is not bound to thermal equilibrium.
Emitted light. We consider now the emitted light power ∼ κω c a † a , by plotting the average population of the cavity mode in Fig. 2(c) as a function of V L . We find that the photon population also increases with bias voltage in a step-like manner [50], correlated to the evolution of the electronic current [46], thus confirming that single-electron tunneling is at the origin of light-emission inside the cavity. From Eq. (2), performing a secular approximation, we derive a rate equation for the photon population P n (t) [33]. A relevant experimental regime in plasmonic cavities is κ Γ and ω c k B T . Since the time between two tunneling events is much longer than the damping time of the cavity, typically the circulating photon leaks out before a new photon is emitted in the cavity. In this limit P 0 ≈ 1, and we find for the other populations P n = Γ 0n /nκ 1, with Γ 0n = α Γ α P n (λ) n F (nω c − µ α ) the cavity 0 to nphotons transition rate induced by a single tunneling event. The expression for a † a = n Γ 0n /κ describes then accurately the emitted power.
Correlation function g (2) . In order to characterize the statistics of the emitted light, we compute the secondorder correlation function g (2) [44,[51][52][53]. Let us begin with the t = 0 case, g (2) (0) = ( n n(n−1)P n )/( n nP n ) 2 . One can readily verify that in thermal equilibrium g (2) (0) = 2. On Fig. 3 we show g (2) (0) as a function of V R and V L , for three different values of λ. As expected, for e|V L − V R | k B T , one always finds the value of 2, corresponding to thermal equilibrium (pink regions on the diagonal V L = V R ). Out of equilibrium we find either bunching g (2) (0) > 1 or antibunching g (2) (0) < 1. The anti-bunching appears for sufficiently strong coupling and is indicated by the blue regions with dashed border (where g (2) (0) = 1). Treating both the electron tunneling and the thermal excitations as weak perturbation to the distribution P n = δ n,0 we obtain an analytical expression for g (2) (0) (see [33] The prediction of this expression agrees very well with the numerical calculations for k B T ≥ 0.1ω c , (cf. thick lines in Fig. 4(a) and (b)). At lower temperature we thus show only the result obtained with the analytical expression, that does not suffer from numerical instability. In the strong coupling regime, we predict a smooth crossover from the equilibrium value g (2) (0) = 2 at low voltage, to an anti-bunching g (2) (0) < 1 regime that appears close to the first inelastic threshold (|eV L | ≈ ω c or |eV R | ≈ ω c ) [33,54]. As temperature decreases the region of antibunching expands, eventually including almost the full bias range [0, 2ω c ], cf. Fig. 4(b). Figure 4(c) shows the minimum value taken by g (2) (0) when minimized on the plane (V L , V R ) for a given value of the coupling λ. One finds that anti-bunching can be observed for λ > 0.17 for an asymmetric bias configuration, thus being at reach of present experiments. The case of symmetric bias is also shown, with anti-bunching beginning at λ > 0.76. Let us discuss now the anti-bunching mechanism for symmetric bias. When P n+1 P n , at lowest order g (2) (0) = 2P 2 /P 2 1 . This expression gives 2 for thermal distribution P n ∼ e −nωc/k B T . Anti-bunching is thus achieved for 2P 2 P 2 1 . At very low temperature the only way to populate the state 2 is either a 2-photons transition (for eV L > 2ω c ), or an electron-tunneling assisted transition from the state 1 to the state 2, controlled by Γ 12 = e −λ 2 λ 2 2 − λ 2 2 α Γ α n F (ω c − µ α )/2 (for eV L > ω c ). As a confirmation one finds that P 1 ≈ P 0 Γ 01 /κ and P 2 ≈ P 1 Γ 12 /2κ, testifying that the result is just due to a balance between the light leaked out of the cavity and the photons emitted in the cavity by the tunneling electrons. From the explicit expression of Γ 12 at low temperature (k B T ω c ) one then finds that the minimum value of g (2) (0) is approximated by Fig. 4(c)], tracing the λ-dependence of Γ 12 . Its vanishing for λ = √ 2 is thus at the origin of the anti-bunching. A similar effect has also been recently reported in dc-biased Josephson junction coupled to microwave resonators [11,12].
Conclusions. Charge fluctuations induced by electronic transport in a molecular single-electronic level are expected to couple strongly to the plasmonic mode formed by an STM tip and the substrate. We derived an expression for the current taking into account the strong damping of these cavities and obtained the current, emitted light intensity, and the correlation function g (2) (t). We showed that when the coupling strength is of the order λ ∼ 1, Franck-Condon steps appear in both the current and the light intensity. Non-classical light can be emitted for a coupling strength in the range 0.17 < λ < 1.8 for bias voltage near the one photon emission threshold. This prediction can be relevant for a series of experiments on STM cavities [16][17][18][19][20][21][22][23][24]. The importance of single-level descriptions has further been reinforced experimentally after our initial submission [56]. The fast evolution of the cavity (κ > k B T ) might question the validity of the Markov approximation. We esti-mate that the non-Markovian contributions are negligible as far as k B T ω c , but investigating low-frequency response of plasmonic cavities could unravel such effects. Another interesting direction is the study of current response in presence of light irradiation of the plasmonic junction. Finally, we note that a local quantum emitter can be used to sense the environment of the molecule with the minimal quantum of energy.
We provide supplementary information about the microscopic Hamiltonian describing the coupling between the plasmon mode of a STM cavity and the electronic single-level of a molecule embedded between the STM tip and the substrate. The dependence of the tunneling current across the junction with the plasmon-molecule coupling strength is analyzed, and shows analogous features to the Franck-Condon blockade regime. The second-order correlation function of the cavity plasmon field is investigated further both analytically and numerically, in the full bias-voltage range and in a regime where dissipation κ of the cavity is larger than the tunneling rate Γ.

Microscopic Hamiltonian
In this section, we provide more details about the microscopic Hamiltonian describing the coupling between a STM local plasmon mode, and a single molecule embedded between the STM tip and the substrate. The form of the interaction Hamiltonian in Coulomb gauge is derived in Eq. (28) of Ref. [1], for the case of a microwave cavity mode coupled to a quantum dot connected to a single electronic reservoir. Similarly to this case, the interaction Hamiltonian between the electronic STM nanocircuit (including two electronic leads and the electronic single-level of the molecule) and the cavity plasmon mode is given by [1] The first term in Eq. (2) stands for the direct (monopolar) coupling between the plasmon mode and charge fluctuations of the dot-level. The corresponding coupling strength is given by the matrix element with φ d ( r) the dot electronic wave function, and V ⊥ ( r) = −iω c C( r0, r) A( l)·d l the photonic pseudo-potential [1] obtained as the work of the cavity electrical field −iω c A( l) ( A is the vector potential) on a path C ( r 0 , r) connecting any reference point r 0 to the point r. The second term in Eq.
(2) describes modulation of the electronic tunneling from the leads to the dot induced by the photonic potential. The corresponding matrix element for this coupling is with φ αk ( r) the Bloch function of the metallic lead α = L, R with quasi-momentum k. Finally, the last term of Eq. (2) stands for the direct coupling between the electrons in the leads and the cavity plasmon mode

Orders of magnitude
The direct (monopolar) term proportional to Λ m (see Eq. (2)) contributes to shifting the dot-level energy by an amount δε 0 = ε 0 −ε 0 = −Λ 2 m / ω c . A simple estimation of the coupling strength Λ m is obtained when neglecting the variation of the photonic pseudo-potential on the scale of the typical extension of the molecular orbitals Λ m ∼ −eV ⊥ ( r m ) ≈ eLE zpm , with r m the location of the molecule, L the typical extension of the transport region, and E zpm zero-point vacuum fluctuations of the cavity electric field.
The third and last term in Eq. (2) stands for processes (allowed by charge conservation) by which electrons in the leads directly decay to the cavity plasmon mode. We arXiv:1907.11269v2 [cond-mat.mes-hall] 30 Oct 2019 have not included explicitly this term in the Hamiltonian written in the paper, since it is equivalent to an additional dissipation channel for the plasmon mode, already taken into account by our phenomenological damping term ∝ κ.
We thus have shown that the monopolar term provides an important (even dominant) contribution to the interaction Hamiltonian in Eq. (2) describing the coupling of a single-electronic level to the cavity plasmon mode. We present in Fig. S1 the electronic current I as a function of the voltage drop eV L in a symmetrical junction (V L = −V R ) for various values of the coupling strength λ. We show that I(V L ) presents steps each times the voltage drop reaches a n-photon emission threshold V L ≈ nω c . The height of those steps increases with λ and after reaching a maximum, decreases again. Furthermore, Fig. S1 shows that as λ increases, the current gets blocked at low voltages, a phenomenon sharing strong analogies to the Franck-Condon blockade regime observed in molecular junctions with phonons [2,3].
We present in Fig. S2 the evolution of the electronic current I at low-voltage eV L = 0.5ω c with coupling strength λ. Using Franck-Condon blockade theory at low temperature (T ω c ) and in the equilibrated regime for the photon population [2], we expect the current to scale as Clearly, while Eq. (6) properly predicts the appearance of inelastic thresholds in the current I(V L ), when V α ≈ nω c . It also predicts a width of the steps that is proportional to temperature T , in contrast to Fig. S1 which exhibits steps with a width given by the dissipation rate of the cavity κ. This is clearly due to the fact that Eq. (6) does not properly include a description of cavity losses.
However, some qualitative features present in Fig. S1 are preserved by Eq. (6). For instance, for voltages below the first inelastic threshold V L < ω c , Eq. (6) predicts an exponential suppression of the low-bias current: ln (I/Γ) ≈ −λ 2 , which is indeed obtained in Fig. S2. Also, Eq. (6) predicts that the typical size of the first in- with η V a voltage of order κ, scales as ∆I ≈ Γλ 2 e −λ 2 , thus reaching a maximum value for λ ≈ 1, which is consistent with the curves of Fig. S1.
To show further the impact of λ on the current suppression, we plot in Fig. S3 the voltage drop V 90 needed to recover 90% of the maximum current as a function of λ. Fig. S3 shows that V 90 increases with λ in a parabolic manner: the higher λ, the stronger the current suppression at low-voltage, and the higher the needed voltage to reach saturation. It is interesting to notice that the classical analogue to Franck-Condon blockade theory [4] predicts that a critical voltage 2ω c λ 2 is needed to unblock the junction, which scales accordingly to V 90 in Fig. S3.

Expression of the dc-tunneling current
The current operator across lead α = L, R is obtained from the continuity equation for tunneling charges as with e the elementary charge. The average current at time t is computed using the density matrix ρ(t) of the full system, including the molecular junction with its electronic reservoirs, and the cavity mode with the external photon bath The full density matrix in Eq. (8) evolves in time with respect to the full HamiltonianH (after the Lang-Firsov transformation has been performed) aṡ The initial condition for the density matrix in Eq. (10) is chosen to be the product of density matrices of the electronic (ρ α=L,R ) and bosonic (ρ b ) environments each supposed to be at equilibrium, times the initial reduced density matrix ρ(0) of system (dot plus cavity-mode). We further define the "bare" Hamiltonian H 0 =H S +H B = H −H I containing the system and baths Hamiltonians. The "perturbation Hamiltonian" V =H I contains the tunneling Hamiltonian (dressed by photons) and the dissipative coupling between the cavity-mode and the photon bath. Integrating Eq. (9) with initial condition Eq. (10), we get an exact expression for the average current with all operators written in interaction picture with respect to H 0 .
We further perform the Born approximation ρ(t ) ≈ ρ L ⊗ ρ R ⊗ ρ b ⊗ ρ(t ) and truncate consistently Eq. (11) at lowest non-vanishing order in the tunneling rates Γ α . We note that the first term in Eq. (11) vanishes due to our choice of initial conditions. We further expand the second term after injecting the expression of the interaction Hamiltonian V written in the paper. After coming back to Schrödinger representation, a subsequent Markov approximation enables to neglect the time dependence of ρ(t ) ≈ ρ(t) in the time integral, and to take the limit t 0 → −∞. We obtain the final expression for the average dc-current I α = lim t→+∞ I α (t) with A(τ )B 0 = Tr e iH0τ Ae −iH0τ Bρ st the average value of the system operators A(τ ) and B(0) written in interaction representation, using as initial condition the stationary density matrix ρ(0) ≡ ρ st .
In order to incorporate the contribution of cavity dissipation to the time evolution of the D and D † operators in Eq. (12), we further replace each term of the form A(τ )B 0 by the correlation function S AB (τ ) = Tr e iHτ Ae −iHτ Bρ st ≡ Tr S Ae Lτ Bρ st . This ansatz is strictly speaking beyond Born approximation, being all orders in the damping rate κ. Equation (3) in the main text of our paper follows after Fourier transformation of the correlation functions S DD † (τ ) and S D † D (τ ).  Figure S3.
Value of V for which 90% of the maximum current is achieved as a function of the coupling strength λ in the case of a fully symmetric junction (VL = −VR = V , ΓL = ΓR = 0.001ωc and kBT = 0.01ωc). In dashed black, expected classical value V90 = 2ωcλ 2 at which the junction gets unblocked.

Rate equation for Pn(t)
We derive in this section the rate equation for the photon population P n (t) = q=0,1 Π qn (t), with Π qn the probability of occupation of the state with q = 0, 1 charge on the dot and n ∈ N photons of the cavity mode. We obtain the following rate equation for Π qn (t) withq = 1, 0 when q = 0, 1, κ ↑ = κn B (ω c ), and The transition rates Γ (qn)(qm) involving photon emission or absorption mediated by electron tunneling events are given by Fermi golden rule and |M mn | 2 is the Franck-Condon matrix element [2] for the considered transition. The stationary solution P st n = lim t→+∞ P n (t) of Eq. (13) enables to calculate the second order correlation function of the cavity mode at equal time with n k ≡ +∞ n=0 P st n n k . In the following, we focus on the case of symmetric bias voltage V L = −V R = V /2, such that Γ (0n)(1m) = Γ (1n)(0m) ≡ Γ nm . A simpler and closed rate equation can be derived in this case for the photon population P n (t), after integrating out the charge degree of freedom of the doṫ Large bias voltage limit VL ≥ 2ωc From now, we focus on the regime κ Γ, for which the cavity damping rate dominates over the dissipation rate induced by tunneling electrons. We also consider low-temperatures T = 0.1ω c , such that κ ↑ κ ↓ ≈ κ. In this regime, the cavity is thus always close to its quantum ground state, namely P 0 ≈ 1 P n≥1 . For V L ≥ 2ω c , we keep in Eq. (17) the leading orders in the photon populations and rates, such thaṫ Eq. (18) can be solved exactly leading to the stationary population Incorporating Eq. (19) into Eq. (16) enables to derive an analytical expression for g (2) (0) that is valid at finite voltage, down to the two-photon emission threshold V L ≥ 2ω c . At the second inelastic threshold V L ≈ 2ω c , emission of photon pairs at equal time is strongly enhanced, thus resulting into a strong bunching of the g (2) (0) versus V curve ( see black dashed curve in Fig.S4 ).
The infinite bias voltage limit V ω c simplifies further in Eq. (19) to such that We thus predict a saturation of g (2) (0) at large voltages to a value g (2) (0) 1, corresponding to strong bunching of the emitted light. For the parameters of Fig.3 in the paper, Γ ≈ 10 −3 ω c and κ ≈ 0.1ω c , such that g (2) (0) ≈ 50, as observed in Fig.S4.

Low-voltage limit VL ≤ 2ωc
In the range of voltages V L ≤ 2ω c , we have to be more careful about how to keep the leading terms in Eq. (17), in particular in order that P st n converges to thermal equilibrium at V L = 0. We found that below the two-photon inelastic threshold (V L ≤ 2ω c ), it is sufficient to keep only states involving at most 2 photons, such that Eq. (17) can be truncated and simplified aṡ P 2 ≈ Γ 02 P 0 + (2κ ↑ + Γ 12 ) P 1 − 2κ ↓ P 2 (22) P 0 ≈ κ ↓ P 1 − (κ ↑ + Γ 01 + Γ 02 ) P 0 (23) The stationary solution of this system of equations provides the rate for the transition from 0 to n photons in the cavity, the corresponding rate for the transition from 1 to 2 photons in the cavity, and ∆ = ( Close to the one-photon inelastic threshold (V L ≈ ω c ), we find that Eq. (25) predicts a minimum value for g (2) (0) given by min g (2) (0) ≈ Γ 12 /Γ 01 ≈ 2 − λ 2 2 /2. Eq. (25) compares well with the numerics as shown in Fig.S4 (see black dotted curve). TIME DEPENDENCE OF g (2) (t) In this section, we are interested in the full time dependence of the second order correlation function g (2) (t) of the plasmonic field g (2) (t) is then computed numerically using the second part of the equation, that is an expression of quantum regression theorem [5]. The curves g (2) (t) are shown in Fig. S5 for various values of λ, at a fixed value of bias voltage close to the first photon-emission threshold V L = −V R = 1.2ω c . In the limit of vanishing plasmonmolecule coupling λ → 0, the plasmon mode gets to thermal equilibrium with the external electromagnetic bath. For a damped electromagnetic cavity at thermal equilibrium, we obtain g (2) (t) = a †2 a 2 e −κt + n B (ω c ) a † a [1 − e −κt ] a † a 2 , with an exponential decay given by κ the damping rate of the cavity. The curve g (2) (t) decays from g (2) (0) = 2 its value at thermal equilibrium, to a value g (2) (t → +∞) = 1 characterizing uncorrelated photon emission events at large times. Note that in this case, since g (2) (t) < g (2) (0) close to t = 0, the probability of emission of two photons at the same initial time is higher than the one of emitting two successive photons, thus characterizing bunching of the emitted light. The fact that g (2) (t) > 1 for all times means that the statistics of emitted light is super-Poissonian. As expected, the obtained blue curve (λ = 0.1) in Fig. S5 is close to the one predicted by Eq. (29).
Upon increasing λ, however, we notice a crossover towards a different regime (see magenta curve for λ = 1.4 in Fig. S5) for which, g (2) (t) > g (2) (0) close to t = 0, namely the probability of emitting two photons at successive times is higher than the one of emitting two photons at the same time. This unambiguously characterizes photon antibunching of the emitted light [6]. Moreover, upon decreasing temperature, antibunching gets almost perfect g (2) (0) ≈ 0 < 1. The fact that g (2) (t) < 1 at any times is a signature that the emitted light has a sub-Poissonian statistics. In the intermediate coupling regime λ ≈ 0.8 (red curve in Fig. S5), g (2) (t) oscillates weakly in time around 1. The statistics of emitted light is thus close to the Poisson distribution obtained for a coherent classic field. We note that after checking the behavior of the full time dependence encoded in g (2) (t) (see Fig.S5), the knowledge of g (2) (0) is sufficient in Fig.S5 to characterize the existence of a crossover from a regime of