Transport of pseudothermal photons through an anharmonic cavity

Under nonequilibrium conditions, quantum optical systems reveal unusual properties that might be distinct from those in condensed matter. The fundamental reason is that photonic eigenstates can have arbitrary occupation numbers, whereas in electronic systems these are limited by the Pauli principle. Here, we address the steady-state transport of pseudothermal photons between two waveguides connected through a cavity with Bose–Hubbard interaction between photons. One of the waveguides is subjected to a broadband incoherent pumping. We predict a continuous transition between the regimes of Lorentzian and Gaussian chaotic light emitted by the cavity. The rich variety of nonequilibrium transport regimes is revealed by the zero-frequency noise. There are three limiting cases, in which the noise-current relation is characterized by a power-law, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$S\propto J^\gamma$$\end{document}S∝Jγ. The Lorentzian light corresponds to Breit-Wigner-like transmission and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\gamma =2$$\end{document}γ=2. The Gaussian regime corresponds to many-body transport with the shot noise (\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\gamma =1$$\end{document}γ=1) at large currents; at low currents, however, we find an unconventional exponent \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\gamma =3/2$$\end{document}γ=3/2 indicating a nontrivial interplay between multi-photon transitions and incoherent pumping. The nonperturbative solution for photon dephasing is obtained in the framework of the Keldysh field theory and Caldeira-Leggett effective action. These findings might be relevant for experiments on photon blockade in superconducting qubits, thermal states transfer, and photon statistics probing.

Experimental and theoretical studies of out of equilibrium cavity and circuit QED have shown remarkable progress during the last decade [1][2][3][4][5] . This research area covers a diverse class of driven-dissipative phenomena and quantum phase transitions [6][7][8][9] . There, observation of photon-photon correlations and quantum state transfer becomes possible in hybrid systems [10][11][12][13][14] where transmission lines are coupled to nonlinear quantum oscillators, such as superconducting transmon qubits or anharmonic cavities. In particular, the photon-photon interaction can be probed in the photon blockade effect as suggested in Ref. 15 , an optical counterpart of the Coulomb blockade in electronic devices. This is an intriguing phenomenon where a driven quantum anharmonic oscillator emits anticorrelated photon 'trains' indicating their sub-Poissonian statistics [16][17][18] .
In this work, we explore the transport of incoherent photons through the cavity with anharmonicity (Kerr energy) smaller than the excitation frequency and bandwidth of the input signal. This is a bosonic counterpart of thermal or voltage biased electronic level with Coulomb interaction 19,20 . Our findings are motivated mostly by experiments on photons statistics [21][22][23] and thermal states propagation 24 that are relevant for various applications such as a microwave non-classical light emission 25 , thermometry 26 , and quantum states transfer 27 .
The moderate anharmonicity does not lead to a well-developed photon blockade, i.e., photon number in the cavity is proportional to the input drive and there is no saturation of the photon current. Nevertheless, it is known that even weak anharmonicity results in the antibunching of photons and their negative correlations 28,29 . In our studies, we find that wideband incoherent pumping induces bunching of photons. We predict a transition from Lorentzian to Gaussian pseudothermal light, as follows from second-order intensity correlator g (2) . However, the emitted light possesses a partial coherence. This results in intriguing behavior in the integral characteristics of the emitted (transmitted) photons such as their zero-frequency noise, written S. Contrary to g (2) that resolves short timescales, S is provided by photon counting during long time intervals. The unusual noise-current relations derived here represent an interplay between photon-photon interaction and strong incoherent drive which brings the system far from equilibrium.
The aim of this work is twofold. First, we analyze nonequilibrium noise, Fano factor, decoherence, and transmission spectrum of the cavity photons. We consider a wide range of behaviors ranging from the single-particle transfer with Breit-Wigner transmission to a many-body regime with unconventional shot noise. The second purpose is methodological. We apply the Keldysh path integral technique 30,31 , and adopt the concept of dissipative Caldeira-Leggett action 32,33 to a driven oscillator with Bose-Hubbard interaction. Although these methods found their application in condensed matter theory a long time ago, in the cavity and circuit QED, they have started to gain attention much later. Nowadays, applications of Keldysh formalism to quantum optics is an active area of research 4,[34][35][36][37] . In our methodology, a nonperturbative solution for quasi-classical fluctuations, which takes into account many-body effects and decoherence of transmitted photons is proposed.
The sketch of the hybrid system under consideration is shown in Fig. 1(a); it is implemented as two open waveguides coupled through an anharmonic cavity 15 . Alternatively, this can be a circuit QED setup, similar to that reported in Ref. 38 , where a multilevel transmon qubit is coupled to superconducting transmission lines.
In Fig. 1(a), the left waveguide is connected to a source of thermal photons, while outgoing light is measured by a photon detector in the right waveguide. The distribution function of incident states in the left waveguide is assumed flat N L,ω = F in the frequency range ω ∈ [ω 0 − �; ω 0 + �] and zero otherwise, as shown in Fig. 1; ω 0 is the cavity mode frequency and is the source linewidth. The outgoing states in the right waveguide can have Lorentzian or Gaussian nonequilibrium distribution functions. The Gaussian distribution is determined by a width κ depending nontrivially on F.
The number of photons F in a particular incident mode can be less or greater than unity. We assume the following ordering of the relevant frequencies where ε is the anharmonicity, Ŵ L and Ŵ R are relaxation rates due to the coupling with the left and the right waveguides, and their sum Ŵ = Ŵ L + Ŵ R determines the bare relaxation rate. (We set Planck constant = 1 hereafter.) The condition (1) is motivated by typical parameters of nonlinear optical cavities 16 where the anharmonicity is induced by a coupling of the cavity mode with a trapped atom. Also, this condition is accessible in superconducting systems based on transmon qubits coupled to transmission lines 14,17,18 . According to (1), antiresonant processes (counter-rotating wave) are neglected and the quantum dynamics is determined by the following U(1) Hamiltonian with continuous symmetry, The first term in (2) is Ĥ ac = ω 0â †â + εâ †â (â †â − 1) ; it governs cavity field dynamics. â † and â are, respectively, boson creation and annihilation operators acting in a basis of Fock states |n� with different photon numbers, and the nonlinearity ε defines Bose-Hubbard interaction strength. Ĥ L = k E L,kb † kb k and Ĥ R = p E R,pĉ † pĉ p describe photon dynamics in the left and right waveguides, respectively; b † k , b k and ĉ † p , ĉ p are creation and H =Ĥ ac +Ĥ L +Ĥ R +Ĥ tL +Ĥ tR . The difference in filling color transparencies stands for a difference in light intensities. The standing wave between blue mirrors is the cavity mode. The current J is measured in the output detector. (b) The spectrum of the incident and outgoing photon states. Incident photons have flat spectrum N L,ω = F of the bandwidth (gray color), nonequilibrium outgoing distribution function have Gaussian shape of width κ . Fluctuating potential φ(t) is induced by the Bose-Hubbard interaction. † kâ + (H.c.) and Ĥ tR = t R pĉ † pâ + (H.c.) , where t L and t R are coupling amplitudes. Broad spectra of E L,k and E R,p determine densities of states, ν L and ν R , at the cavity mode ω 0 and relaxation rates, The measured photon current is defined as the average number of photons n t 0 that have passed during the measurement time t 0 , J = 1 t 0 �n t 0 � . Brackets denote an average for a quantum mechanical operator Ô with a nonequilibrium density matrix ρ , �Ô� = Tr[Ôρ] . The trace is calculated with the help of the nonequilibrium Keldysh technique and path integrals over complex boson fields. For the limit (1), the wideband drive induces the current J = 2F Ŵ L Ŵ R Ŵ that does not depend on the interaction ε . The average photon number in the cavity, �â †â � = F Ŵ L Ŵ , does not saturate as F increases. This value of �â †â � might be associated with effective temperature T eff ∼ F Ŵ L Ŵ ω 0 , however, the incident photons are not distributed by the Bose-Einstein law and the cavity has non-Gibbs density matrix.
Zero frequency noise of the photon current is defined as the second cumulant of counted photons during t 0 interval, S = 1 t 0 �n 2 t 0 � − �n t 0 � 2 . Below we show that the noise-current relation, S(J), and transmission probability function T ω , which appears in a photonic counterpart of the Landauer formula, depend nontrivially on ε . These quantities, as well as g (2) -correlator and Fano factor, demonstrate the richness of Lorentzian-to-Gaussian crossover in a pseudothermal light.

Results
Non-equilibrium effective action. We reduce the description to the Keldysh effective action S eff [φ] for the Hubbard-Stratonovich real field φ(t) . It decouples the Bose-Hubbard interaction term. The dynamics of φ(t) determines the decoherence of transmitted photons. A particular configuration of φ(t) has a contribution ∼ e iS eff [φ] in the partition function. The dynamics of φ(t) is influenced by different multi-photon transitions in the cavity shown in Fig. 2. The particle number operator n =â †â does commute with the cavity Hamiltonian, [Ĥ ac ,n] = 0 . Hence, its Fock states |n� are the same as those of the harmonic oscillator, i.e., they can be represented through the canonical coordinate and Hermite polynomials. The wavefunctions in canonical coordinates should not be confused with multi-mode field configurations in the real space. Contrary to harmonic oscillator, the eigenvalues spectrum of Ĥ ac is a non-equidistant, E n = ω 0 n + εn(n − 1) . It is shown that a description of multi-photon transitions in the nonlinear single-mode cavity is reduced to an effective theory for φ(t) . It has a stochastic behavior, which resembles Kubo model 39 of an oscillator with random modulated frequency. For the sake of compactness, the results of this Section are presented for the symmetric setup with Ŵ L = Ŵ R = Ŵ/2 . We find that the stationary part of φ(t) is a nonequilibrium saddle point of the action. Its value defines a shift of the cavity mode frequency by the occupation number F and the anharmonicity, ω ac = ω 0 + 1 2 (F − 2)ε . A nonperturbative solution for the decoherence follows from Gaussian theory for fluctuations near the saddle point. In our approach, we distinguish the fields φ + (t) and φ − (t) residing on the upward and backward branches of the Keldysh contour K . Hence, the stochastic and quantum phases are introduced as )dt ′ , respectively. Their dynamics is governed by the action of Caldeira-Leggett type 32,33 , Figure 2. Sketch of multi-photon transitions in the cavity with a negative anharmonicity ε < 0 . Eigenstates |n� have non-equidistant spectrum E n = ω 0 n + εn(n − 1) where n is the respective photon number. A multi-photon transition from |0� to |n� is accompanied by the absorption of n photons of the frequency ω n = ω 0 + ε(n − 1) . The decrease of ω n with n for ε < 0 is illustrated by a change of color from blue to magenta. where Fourier transformed phases read as � ω = �(t)e iωt dt and ϕ ω = ϕ(t)e iωt dt . The first term is the 'kinetic' part, whereas the second term is Caldeira-Leggett dissipative part. It is written as a matrix with Keldysh causality structure.
In the equilibrium situation, the many-body density matrix is ρ = is bosonic occupation number. In our situation with flat distribution functions, we find that dissipative terms vanish in the limit of large , α R ω = (α A ω ) * = 4i π� 3 εŴF 2 . Oppositely, the fluctuational Keldysh term does not vanish and depends non-linearly on the occupation number F, α K ω = 4 ŴF(F+2) 4Ŵ 2 +ω 2 . The distribution function in the anharmonic cavity is found, demonstrates a break of the fluctuation-dissipation theorem in our nonequilibrium regime.
Decoherence. Caldeira-Leggett action is reduced to a Langevin equation 40 d dt �(t) = εξ(t) for the stochastic phase of transmitted photons. The random force in the r.h.s. has the following correlator, �ξ(t)ξ(0)� = 1 4 α K (t) , where α K (t) = F(F + 2)e −2Ŵ|t| is the time-resolved Keldysh part from (3). The non-Gaussian effects of the noise can be sensitive to the dynamics of ϕ ; this is beyond the scope of our consideration. The non-Gaussian effects might be important in the limit of large ε > ω 0 when the anharmonic cavity becomes a two-level system.
The decoherence of photons is determined by the envelope The same structure of a correlator appears in the Kubo model. There is Gaussian decay at short timescale with the rate κ = 1 2 ε √ F(F + 2) , which increases with the anharmonicity and is nonanalytic by F due to the square root. The behavior of z(t) plays a central role in our findings because it determines all characteristics of the outgoing light. We note that the Gaussian decay of the envelope appears in the spin-boson model with the sub-Ohmic dissipative environment with 1/ω spectrum 41 . The effective sub-Ohmic spectrum of phase fluctuations in our problem is due to multi-photon transitions.
Intensity correlator. The second-order correlator of outgoing photons, which can be probed in Hanbury-

Brown and Twiss interferometry 42 , is
Thus, the emitted signal is chaotic pseudothermal light that can be Lorentzian or Gaussian depending on the anharmonicity ε and mean photon current. Namely, the Lorentzian light with g (2) (t) ≈ 1 + e −2Ŵt occurs at J ≪ J * ε , whereas the Gaussian, g (2) (t) = 1 + e −2κ 2 t 2 , at J ≫ J * ε is found. The value of J * ε is related to the ratio between the anharmonicity and relaxation as J * ε = Ŵ 2 1 + 4 Ŵ 2 ε 2 − 1 . The correlator at zero time g (2) (0) > 1 is related to the bunching and g (2) (0) < 1 to the antibunching of emitted photons 43 . In our case, g (2) = 2 that is consistent with that pseudothermal photons are bunched. This occurs because the incoherent drive with a flat spectrum resembles the high-temperature limit of Bose-Einstein distribution where photon number fluctuations in a single mode are proportional to mean photon number squared.
However, the presence of coherence in the emitted light suppresses positive correlations. As shown below this is reflected in the emergence of the shot noise and in the suppression of the Fano factor. This is also indicated by the fast decay of g (2) (t) to the unity, which is associated with a coherent light. The decay occurs at the timescale t ∼ 1/κ , which becomes very short at large F and ε.

Transmission spectrum. It is found that the single-photon Breit-Wigner transmission function
Here, τ (t) is time-resolved τ ω . Similarly to g (2) , the Lorentzian transmission spectrum τ ω is changed to Gaussian at J ≫ J * ε which reads This is a nonperturbative result with respect to the interaction and incoherent drive. This transmission describes nonequilibrium and many-body photon transfer. The crossover between Lorentzian and Gaussian spectra in T ω , when F increases, is shown in Fig. 3; results are obtained after numerical integration. It should be noted that the photon-photon interaction redistributes the spectral weight of T ω around ω ac , keeping its integral value constant. The latter explains that fact that the anharmonicity does not change the current in the wideband input drive regime.
Noise. Correlations in the output field are determined by g (2) (t) . A relevant information is given by a short timescale of t 0 ≪ 1/ Ŵ . Distinctly, the zero-frequency noise 44 contains information about the long timescale, Scientific Reports | (2021) 11:8328 | https://doi.org/10.1038/s41598-021-87536-w www.nature.com/scientificreports/ t 0 ≫ 1/ Ŵ . The indicated change of the asymptotic behaviors in z(t) or T ω is reflected in the richness of a lowfrequency noise-current relation S(J). The Lorentzian emitted light shows the quadratic noise-current relation, This scaling resembles a situation of equilibrium radiation in a cavity at a very high temperature where the variation, ��n 2 �� = �n 2 � − �n� 2 , and average photon number are related as ��n 2 �� ∼ �n� 2 . The upper boundary for the current J * ε can be large or small compared to Ŵ depending on the ratio ε/ Ŵ . As shown in Fig. 4(a), the region of small ε corresponds to a Lorentzian light with the exponent γ ≈ 2 in the noise-current relation represented as S ∝ J γ . We note that the quadratic dependence and the Lorentzian spectrum of S therm at finite ω (see Section "Generating functional method" and Eq. (39) for details), is analogous to a modulation noise in quantum point contacts 45 .
In the many-body regime with Gaussian light the following expression is found The absence of the quadratic scaling in S(J) indicates the shot noise behavior. Also, this means that the Gaussian pseudothermal light has a partial coherence. In the limit of strong anharmonicity, ε ≫ Ŵ , the scale J * ε ∼ Ŵ 3 ε 2 is small compared to Ŵ . Consequently, the ratio Ŵ/J in the square root of (8) can be both large or small compared to unity. This shows that two different asymptotical regimes of the Gaussian noise do exist, they are associated with Ŵ/J → 0 and Ŵ/J → ∞ limits. These limits have a continuous crossover between each other at Ŵ/J ∼ 1.
If the drive is sufficiently strong, such that J ≫ Ŵ holds, then we arrive at the linear shot noise-current relation We note that the scalings similar to S shot ∝ J and S therm ∝ J 2 were discussed in Refs. 46,47 . Namely, the photon current emitted by a coherent source and propagated through a random media shows a linear noisecurrent relation if the scattering region is absorbing, and quadratic if it is amplifying.
At low currents, Ŵ ≫ J ≫ J * ε , the ratio Ŵ/J in (8) dominates over the unity and we arrive at unconventional scaling of the noise, S ′ shot ∝ J 3/2 . This is one of the remarkable findings of this work. The change of the scaling exponent in S ∝ J γ from γ = 2 to γ = 3/2 as ε increases can be understood as emerging of a partial coherence in the propagated light. The difference between these nonequilibrium fixed points is illustrated by blue and orange curves in Figs. 5(a) and 5(b). In the limit of weak anharmonicity, ε ≪ Ŵ , the asymptotic behavior of S ∝ J 3/2 does not appear. Here, we obtain a simple crossover between the thermal-like behavior and the shot noise. It is demonstrated as blue curves, which start at γ = 2 and saturate smoothly at γ = 1 . Parameter domains corresponding to different asymptotic limits are collected in Table 1.
Fano factor. It is instructive to analyze the Fano factor, the ratio of the low-frequency noise and current, F = S/J . In mesoscopic physics, F defines an effective charge that is transmitted coherently. This applies, e.g., to single-electron transistors 48 , to helical electrons scattering 49,50 and to multiple Andreev reflection in super- Contours of constant γ show a crossover between Lorentzian and Gaussian pseudothermal light. The Lorentzian light corresponds to singleparticle transport with Breit-Wigner transmission and γ = 2 , and the Gaussian is to many-body transport and the shot noise. The fractional γ = 3/2 asymptotic appears at low currents and large anharmonicity ε/ Ŵ > 1 .     51 . In our optical system, it can be associated with the emergence of positive ( F > 1 ) and negative ( F < 1 ) correlations between photons 52,53 .
In the regime of conventional shot noise we find universal result F shot = π 8 Ŵ ε ≪ 1 that does not depend on the drive amplitude. The unconventional shot noise S ′ shot ∝ J 3/2 corresponds to Fano factor F ′ shot = π 8 ŴJ ε 2 which is also smaller than unity. The regime of Lorentzian light yields F therm = J 2Ŵ which can be less or greater than unity. Thus, the Fano factor represents the complexity of nonequilibrium properties of the emitted light. Its color plot, the contours with F = 1 (black curve), and J * ε (white dashed curve) are presented in Fig. 4(b).

Discussion and outlook
Studies of driven-dissipative optical systems with Bose-Hubbard interaction is a large research area. In particular, incoherent drive effects were explored in experiments on the photon blockade 17,18 and thermal states propagation 24 . The photon blockade is characterized by the antibunching of photons and their sub-Poissonian statistics. The thermal states propagation is an alternative regime, which is characterized by a suppressed photon blockade. In that experiment 24 , the bunching effect and super-Poissonian statistics of emitted photons were observed. In this work, we investigated how the incoherent pump influences the photon transport in a wideband limit and arbitrary excitations number. Namely, the central question of this work is how does a photonic counterpart of thermal transport through an electronic quantum dot with a non-equidistant spectrum look like? This question seems actual in view of increasing interest to novel phenomena and phases in nonequilibrium cavity and circuit QED where the fluctuation-dissipation relation and detailed balance condition are violated. The bandwidth of the drive signal is assumed to be greater than characteristic relaxation rates. In this limit, we found that the increase of the Bose-Hubbard interaction and incoherent pump induce an intriguing transition from effectively thermal fluctuations to partially coherent nonequilibrium noise. A complexity of this transition is demonstrated in Fig. 5(b) where cross-sections of the phase diagram from Fig. 4(a) are shown for particular values of the interaction.
First, the Lorentzian transmission function and noise spectrum are transformed into Gaussian. The latter is a signature of many-body interaction in the anharmonic cavity and the consequence of transmitted photons decoherence. The time-resolved intensity correlator g (2) and photons dephasing reveal transitions from exponential to Gaussian decay laws demonstrating a certain similarity with Kubo model of an oscillator with stochastic frequency and also with the spin-boson model with sub-Ohmic 1/ω spectrum.
Second, the low-frequency noise S and the Fano factor are nonlinear functions of the drive amplitude, or, the average transmitted current. The scaling exponent in the noise-current ratio, S ∝ J γ , can have three universal values. They indicate (i) single-photon transfer and effectively thermal state noise with γ = 2 , and (ii) partial coherence in strongly nonequilibrium limit with the shot noise ( γ = 1 ) and suppression of the Fano factor that becomes smaller than unity, and (iii) unconventional noise-current relation with the fractional exponent γ = 3/2 . In the regime (iii), that is found for low currents, the fluctuations of outgoing states are contributed by two competing components, the coherent emission and pseudothermal dissipative dynamics in a large Fock space of the cavity mode. This regime is supposed to be relevant for experiments reported in the Ref. 24 , where the difference between fluctuations in thermal radiation and the Poissonian regime was observed, and also in Refs. [21][22][23] , where the statistics of incident photons was measured by a continuous wave mixing.
In our solution, we found a single maximum in the transmission spectrum at the frequency ω ac . This is in contrast to the case of the coherent drive at a particular frequency, where multi-photon supersplittings 54 are clearly observed in experiments 14 . Their absence in our solution is explained by the averaging over a large amount of multi-photon transitions induced by the wideband and noisy drive.
The effective Caldeira-Leggett action proposed here describes fluctuations in the Gaussian (quadratic) order approximation. This approach provides a quasi-classical theory for noise represented as symmetrized correlators due to the Keldysh time ordering. The non-Gaussian fluctuations effects are the intriguing direction of recent studies related to non-linear quantum environment 35 and dynamical Coulomb blockade 55 ; in particular, the possibility of different time orderings on the Keldysh contour becomes important in the description of quantum noise and nonsymmetrized correlators.

Summary
To conclude, we investigated steady-state photon transport in the anharmonic single-mode oscillator. It can be transmon qubit or nonlinear optical cavity incoherently driven and coupled to the input and output waveguides. We predict a rich behavior of the noise and intensity correlators in the nonequilibrium. The nonequilibrium conditions reveal intriguing effects due to the photon-photon interaction, they include pseudothermal fluctuations with partial coherence, Gaussian transmission spectrum with effectively sub-Ohmic dephasing, and unconventional shot noise with fractional power-law scaling.
The assumption of the wideband noisy drive field allows for a compact analytical solution. The decoherence due to the photon-photon interaction is described within a nonperturbative approach reduced to Caldeira-Leggett effective action. In this many-body problem, a stochastic phase of transmitted photons depends on various multi-photon transitions. We expect that the nontrivial behavior of the noise might not have a direct analogy in condensed matter systems.
The methodology is based on the Keldysh field theory. The solution for a system far from equilibrium is obtained. Our results give an insight into photon transport in the cavity and circuit QED, where the interaction and nonequilibrium dynamics appear on an equal footing. They can be generalized, in particular, for drivendissipative Bose-Hubbard lattices 56,57 . It is suggested that these findings are experimentally accessible by methods based on photon counting 6,17 , homodyne detection of intensity correlations 16

Methods
In this part of the paper, we present the methodology of calculations. We start from a formulation of nonequilibrium field theory that corresponds to a microscopic Hamiltonian (2). It is based on the Keldysh Green functions approach and Hubbard-Stratonovich transformation decoupling the non-Gaussian interaction term. In particular, the Landauer formula for photon current is derived within the Keldysh formalism. Then, the Caldeira-Leggett action is formulated and applied to calculations of decoherence, transmission function, noise spectrum, and g (2) -correlator. The last term in the exponent corresponds to fluctuating potential φ(t) that randomly modulates oscillator frequency. As a result of this transformation, the cavity action becomes the following:

Effective functional.
With the use of this representation, the total action is formulated through matrix Green functions. Here, a transformation to the usual time axis t ∈ [−∞, ∞] doubles the number of variables, , where these fields reside on upward or backward parts of the Keldysh contour labeled by ± indices. This is associated with Keldysh τ-space parametrized by Pauli matrices τ x,y,z and identity matrix τ 0 . After the standard Keldysh rotation to classical and quantum components for complex fields, , we arrive at the following representation of the total action: Each element of the matrix is a block in τ-space that possesses causality structure. Let us consider the left waveguide's block: The retarded (advanced) Green functions are G R(A) = 1/(ω − E L,k ± i0) ), they differ by the sign of the infinitesimal frequency shift i0 along the imaginary axis. The Keldysh component 2i0f L,k encodes a distribution function N L,k of incident photons as f L,k = 2N L,k + 1 . The similar definition applies for Ǧ −1 R . The local in time inverted Green function of the cavity is nonstationary because it involves φ c (t) and φ q (t) . The nondiagonal elements t L,Rτx in (12) mix the cavity and waveguides' fields corresponding to k and p modes. After the integration over b k and c p , we obtain the effective action for the anharmonic cavity, where the inverted Green function becomes a nonlocal in time due to the self-energy � (t) = dω 2π e −iωt� ω . It has the causality structure, , the retarded and advanced components are � R(A) ω = ±iŴ where the total relaxation Ŵ = Ŵ L + Ŵ R is a sum of the rates determined by the couplings to left and right waveguides,  + 1) , brings information on the nonequilibrium distribution function in the cavity mode N ac,ω = Ŵ L Ŵ N L,ω + Ŵ R Ŵ N R,ω . The integration over the cavity fields ā, a provides the effective action for the fluctuating potential The trace and logarithm here are determined over discretized time axis. As long as this action is non-Gaussian due to the matrix logarithm, we apply Gaussian expansion near a saddle point.

Saddle point equation.
A configuration of φ c (t) that determines a saddle point of S eff is given by the zero variation with respect to φ q , δ δφ q S eff [φ c , φ q ] = 0 . Applying this to (17), we find φ c dt − iεTrǦ ac [φ c ] = 0 . Here, we use the standard notation for trace Tr(f (t − t ′ )) = f (0) dt = dω 2π f ω dt , where the time integral stands for a large constant that cancels out in the saddle point equation. In our system, the saddle point configuration is associated with the static part of the field, φ 0 . An explicit form of this equation is In our case of flat N L,ω = F in the range ω ∈ [ω 0 − �; ω 0 + �] with � ≫ Ŵ, ε , and the absence of incident modes in the right waveguide, i.e., N R,ω = 0 , the integral in (18) gives a constant that does not depend on φ 0 . The equation (18) becomes linear and we find φ 0 = ε + Fε Ŵ L Ŵ . Thus, the photon Green function at the saddle point is found after the inversion of Ǧ −1 ac [φ] from (16) with φ = φ 0 : Caldeira-Leggett approach for fluctuations. Non-stationary configurations of the field, δφ c (t) = φ c (t) − φ 0 , determine many-body transitions which induce decoherence. Let us go back to the effective action S eff . It is non-Gaussian because of the matrix logarithm. We apply quasi-classical second-order expansion of the logarithm by stochastic, δφ c , and quantum, φ q , fluctuations: This approximation for S eff provides Caldeira-Leggett action (3) which rules fluctuational and dissipative dynamics of φ(t) . For nonequilibrium distributions indicated above, we find from (20) that the retarded and advanced parts in (3) read They are not universal, i.e., they depend on the bandwidth and vanish as 1/� 3 . The Keldysh part, however, is non-zero This fact demonstrates a break of fluctuation-dissipation relation out of the equilibrium, i.e., We note that the action S CL determines stochastic Langevin equation 1 ε + 1 2 α R ω δφ c,ω = ξ ω where the stochastic force ξ has the correlator given by �ξ −ω ξ ω � = 1 4 α K ω . Addressing the photons decoherence in Keldysh formalism, we introduce stochastic �(t) = t −∞ δφ c (t ′ )dt ′ and quantum ϕ(t) = t −∞ φ q (t ′ )dt ′ phases. In our case with vanishing α R , the Langevin equation for the transmitted photon phase is reduced to that corresponds to a frictionless drift.
Calculations of transmission functions. Non-interacting limit. Lorentzian spectrum. Analogy with Landauer formula. We address a photon transport in a spirit of the Landauer formula determined by the transmission probability T ω . We start from a single-particle problem at ε = 0 . This limit is exactly solvable. After that, we study the case of ε = 0 taking into account many-body interaction utilizing Caldeira-Leggett action (3).
Let us define the photon current as time derivative of the photon number operator n R = pĉ † pĉ p in the right waveguide as J = d dt �n R � = i �[Ĥ,n R ]� . We note, that tunnelling Hamiltonians can be represented as (20) www.nature.com/scientificreports/ H tL = t LB †â + t * Lâ †B and Ĥ tR = t RĈ †â + t * Râ †Ĉ where B = kb k and Ĉ = pĉ p are 'local' operators of waveguide fields. In this representation, the photon current reads as We employ Keldysh Green functions technique to calculate these averages. In the non-interacting case, the action reads The inverted Green functions of the cavity mode is Ǧ −1 0,ω = (ω − ω 0 )τ x . Green functions of left and right waveguides are, respectively, ǧ L,ω = − kǦ L,ω,k and ǧ R,ω = − pǦ R,ω,p . They act in space of variables ψ = [B ω , a ω , C ω ] . The summations over k and p are found after replacements k,p → ν L,R dE.
Interacting case: Gauge transformation. In this part, we present a generalization on the many-body problem where ε = 0 addressing the situation of flat distribution functions N L,ω = F and N R,ω = 0 . The action with stochastic potential becomes local in time due to flat distributions In contrast to the action in the non-interacting limit S (0) , the Green function Ǧ −1 0 in the interacting case is replaced with Ǧ −1 )τ x − δφ q (t)τ 0 . The locality in time gives an advantage in analytic calculations. Namely, it is possible to get rid of time-dependent δφ c (t) and φ q (t) by the gauge transformation. This transformation from ψ = [B, a, C] T to new fields, ψ ′ = [B ′ , a ′ , C ′ ] T , is distinct on upward and backward parts of K and reads It provides a nonperturbative solution for photon decoherence. A dynamics of these new fields ψ ′ is ruled by the action S (0) found for ε = 0 . Hence, the solution for cavity-to-right waveguide propagator ǧ acR at ε = 0 is obtained from non-interacting result ǧ 0R from (26) when the gauge inversed to (29) is applied: This function is not stationary due to the fluctuating potential. After the averaging ǧ acR (t, t ′ ) with respect to S CL , we arrive at the stationary Green function that depends on the time difference The relevant modification appears in the envelope related to stochastic fluctuations z(t) = �e −i(�(t)−�(0)) � . The average with Gaussian action S CL yields z(t) = exp(−D(t)) where the symmetrized autocorrelation function is D(t) = 1 2 (��(t)�(0)� + ��(0)�(t)�) − �� 2 (0)� . With the use of Langevin equation (23), we find that found after the Fourier transform of (22), gives  (31) are due to the mixing of fields on different parts of the contour K ; they read ��(t)ϕ(0)� = iεtθ(t) and �ϕ(t)�(0)� = −iεtθ(−t) . In our Gaussian theory, classical-quantum terms yield simply a complex phase e −iε(t−t ′ ) that results in a shift of the cavity mode frequency a quantum counterpart of the nonlinear frequency shift in a classical anharmonic oscillator.
Gaussian spectrum. The result (31) provides the expression for the transmission function in the time domain Here, τ (t) is Fourier transformed Lorentzian transmission with the maximum at ω ac . The Fourier transform of (34) reads T ω = z ω−ω ′ τ ω ′ dω ′ 2π . Analytic calculation of this integral is challenging. Nevertheless, we find asymptotic expressions for the strongly nonequilibrium regime with the use of the saddle point method. In this solution, the correlator in the exponent of z(t) is approximated as D(t) ≈ κ 2 t 2 . Here, new relaxation rate It shows that the coherence loss rate grows non-linearly as F increases. The Fourier transformation gives the transmission spectrum of Gaussian shape: . This result applies to a regime of strong driving or strong interaction, such that the condition κ ≫ Ŵ is satisfied. Let us analyze this condition for the symmetric case of Ŵ L = Ŵ R = Ŵ/2 . The decay rate can be equivalently expressed through the current as κ = ε Ŵ JŴ + J 2 where J = FŴ/2 does not depend on anharmonicity in our problem with the wide spectrum of the incoherent drive. The condition κ ≫ Ŵ is resolved as J ≫ J * ε : The scale J * ε corresponds to the lower boundary for the current at a given ε when it drives the system out of the equilibrium and induces Gaussian chaotic light emission. If the current is small, J ≪ J * ε , then the system is in a fixed point related to thermal behavior, as shown in the phase diagram in Fig. 4(a). We learn from (35) that there is two asymptotics in J * ε that distinguish weak, ε ≪ Ŵ , and strong, ε ≫ Ŵ , interaction limits. Both of them overlap with the Lorentzian and Gaussian sectors.
Photon current and cavity photon number. As it follows from (27), the photon current at N L,ω = F and N R,ω = 0 reads The average photon number in this limit, is obtained as �â †â � = iG < ac,ω dω 2π where the 'greater' Green function is G < ac,ω = 1 2 (G K ac,ω − G R ac,ω + G A ac,ω ) . We note that the frequency shift due to the anharmonicity can be represented via (37) as ω ac = ω 0 − ε + �â †â �ε . As it should be, the same value can be found from a mean-field approach for Bose-Hubbard interaction where â †ââ †â →â †â �â †â �.
The current given by (36) is not modified by the interaction in the case of flat distribution functions. As follows from the representation of the current through g < acR (t, t) at coincident times, we find J = Fz(0) dω 2π τ ω from (34). As long as z(0) = 1 for any interaction, J does not depend on ε . The same logic applies to �â †â � which also does not change as ε increases.
Calculations of the noise and intensity fluctuations. Generating functional method. We start calculations of noise from the non-interacting limit and then generalize for the interacting case. Cumulants of transmitted photons can be calculated through variations of the generating functional logarithm, ln Z[η] . The generating functional is given by the path integral (32) (24) through fields a and C on opposite branches (labeled by τ z = ±1 ) of the contour K as J =ψJψ . The integration over ψ , ψ in (38) results in ln Z[η] = Tr ln(1 + GηJ) . Photon current J is found as first variation of ln Z[η] by η , where it is set to zero η → 0 . This gives J = dω 2π itr G ω J that is reduced to Landauer formula (27); here, " tr " stands for trace over τ , L, R, and oscillator indices. In order to find the spectrum of the symmetrized noise, S (0) ω = 1 2 (�δJ −ω δJ ω � + �δJ ω δJ −ω �) , we need to find second variation of ln Z[η] by η −ω and η ω . In our particular case of flat N L,ω , the spectrum is Lorentzian Asymptotic behavior of the noise-current relation. Hereafter, we express F through the current J according to (36) and suppose that Ŵ L and Ŵ R can be unequal. This allows analyzing zero-frequency noise-current ratio, S(J), where both S and J are measured by the output detector. The zero-frequency noise found in (39) has quadratic scaling, S therm = J 2 2Ŵ at J ≪ J * ε , that corresponds to Lorentzian light emitted into the output waveguide. The decoherence is important at J ≫ J * ε when pseudothermal output light becomes Gaussian and the thermal noise is changed to shot noise. Technically, the inclusion of the stochastic �(t) in the noise calculation is performed with the gauge transformation, similarly to that applied for the transmission function T(t). Symmetrized correlator of the noise, S(t) = 1 2 (�δJ(t)δJ(0)� + �δJ(0)δJ(t)�) , reads S(t) = z 2 (t)S (0) (t) where S (0) (t) = 1 2 J 2 e −2Ŵ|t| is time-resolved Lorentzian correlator found in the non-interacting limit (39). The shot noise expression follows from the Fourier transform of S(t) where the correlator in the exponent of z(t) is Gaussian D(t) ≈ κ 2 t 2 . The low-frequency result, S ≡ S ω=0 , is given by S = J 2 ∞ 0 e −2Ŵ|t|−2κ 2 t 2 . Assuming κ ≫ Ŵ , that is equivalent to J ≫ J * ε , we calculate this integral and arrive at one of central results Let us analyze it in details. If the interaction is weak, ε ≪ Ŵ , then we always have Ŵ R /J ≪ 1 according to J * ε ∼ Ŵ 2 ε from (35). Thus, in the nonequilibrium regime at weak anharmonicity, we find conventional shot noise from (40) with a linear noise-current ratio S shot = π 2 Ŵ R ε J at J ≫ Ŵ 2 ε . This result is considered as a nonequilibrium fixed point with the exponent γ = 1.
For strong interaction, ε ≫ Ŵ L,R , the thermal noise sector in Fig. 4(a) shrinks because of the vanishing J * ε ∼ Ŵ 3 ε 2 . This means that we can go to low current domain where J ≪ Ŵ R . Here, the noise-current relation (40) has the following asymptotic This non-analytical dependence occurs at small J ∼ J * ε , i.e., it is parametrically close to zero current line J = 0 in Fig. 4(a) at large ε . This feature is also demonstrated in Fig. 5(b) for ε > Ŵ where orange and red curves drop rapidly from γ = 2 to γ = 3/2 . Then, curves saturate to conventional shot noise with γ = 1 as the current increases.
There are smooth crossovers between S therm , S ′ shot and S shot at intermediate ε ∼ Ŵ as a consequence of fluctuations in zero-dimensional cavity. At weak anharmonicity, there is no S ′ shot behavior. Instead, there is a crossover from S therm to S shot at J ∼ Ŵ 2 ε . It is shown in Fig. 5(b) where blue curve decays smoothly from γ = 2 to γ = 1.
Intensity correlators. Time-resolved intensity-intensity correlator in the right waveguide, is an indicator for bunching or antibunching of photons in the output field. In our solution we assumed that anomalous term �Ĉ(t)Ĉ(0)� is zero, and this four-point correlator is defined through the two-point one, g (2) (t) = 1 + |g (1) (t)| 2 , according to the Wick theorem applicable in our quasi-classical solution. We find g (1) (t) = z(t)ǧ RR (t)/ǧ RR (0) where ǧ RR (t) ∝ e −Ŵt is the Lorentzian propagator while z(t) involves ε and F.
Time dependence of the correlator shows exponential decay in the non-interacting limit g (2) (t) = 1 + e −2Ŵt at t > 0 . In the interacting case, we find   www.nature.com/scientificreports/ There is Gaussian law of the correlations decay g (2) (t) = 1 + e −2κ 2 t 2 for t ≪ 1/ Ŵ . For large timescale, t ≫ 1/ Ŵ , there is exponential decay g (2) (t) = 1 + e −2 κ 2 Ŵ t . As follows from the condition κ ≫ Ŵ on the nonequilibrium, the decay rate in the latter case exceeds that from the non-interacting limit as κ 2 Ŵ ≫ Ŵ.