The Circuit Quantum Electrodynamical Josephson Interferometer

Arrays of circuit cavities offer fascinating perspectives for exploring quantum many-body systems in a driven dissipative regime where excitation losses are continuously compensated by coherent input drives. Here we investigate a system consisting of three transmission line resonators, where the two outer ones are driven by coherent input sources and the central resonator interacts with a superconducting qubit. Whereas a low excitation number regime of such a device has been considered previously with a numerical integration, we here specifically address the high excitation density regime. We present analytical approximations to these regimes in the form of two methods. The first method is a Bogoliubov or linear expansion in quantum fluctuations which can be understood as an approximation for weak nonlinearities. As the second method we introduce a combination of mean-field decoupling for the photon tunneling with an exact approach to a driven Kerr nonlinearity which can be understood as an approximation for low tunneling rates. In contrast to the low excitation regime we find that for high excitation numbers the anti-bunching of output photons from the central cavity does not monotonously disappear as the tunnel coupling between the resonators is increased.


Introduction
In recent years, a new direction of research in cavity quantum electrodynamics (cavity-QED) has developed, in which multiple cavities that are coupled via the exchange of photons are considered. Such setups are particularly intriguing if the cavities are connected to form regular arrays and if the strong coupling regime is achieved in each cavity of the array. Such devices would give rise to novel structures where coherent light matter interactions exceed dissipative processes simultaneously in multiple locations of the array [1][2][3][4][5][6].
Whereas it is rather challenging to build mutually resonant high finesse cavities in the optical regime, it is for microwave photons perfectly feasible to engineer large arrays of mutually resonant cavities on one chip in an architecture known as circuit-QED [6,7]. Here multiple superconducting transmission line resonators with virtually identical lengths in the centimeter range can be coupled via capacitors or inductive links [6].
The individual transmission line resonators can feature strong optical nonlinearities by coupling them to superconducting qubits such as transmons [8,9] or phase qubits [7].
Yet in all experiments that involve light-matter interactions, some photons will inevitably be lost from the structure due to imperfect light confinement or emitter relaxation. To compensate for such losses, cavity arrays are thus most naturally studied in a regime where an input drive continuously replaces the dissipated excitations. This mode of operation eventually gives rise to a driven dissipative regime, where the dynamical balance of loading and loss processes leads to the emergence of stationary states [10][11][12][13].
The properties of these stationary states may however vary significantly if one changes system parameters such as the photon tunneling rate between cavities, the light-matter interaction in a cavity or even the relative phase between a pair of coherent input drives [12,14].
A device that is ideally suited for studying the effect of relative phases between input drives is the so called quantum optical Josephson interferometer [14], which consists of two coherently driven linear cavities connected through a central cavity with a single-photon nonlinearity. Here, the interplay between tunneling and interactions in the steady state of the system has been analyzed for regimes of weak input drives where the number of excitations in each cavity is rather small by Gerace et al. [14]. For opposite phases of the input drives one finds a destructive interference which suppresses population of the central cavity, whereas for input drives in phase the central cavity is populated with anti-bunched excitations due to its strong nonlinearity.
For the considered regime of low excitation numbers the approach [14] employed a full numerical analysis relying on an excitation number truncation of the Hilbert spaces. For high input drives however, the theoretical description of the system poses a challenge as conventional numerical methods quickly become computationally infeasible. Hence, such a numerical approach cannot describe regimes with more intense input drives where excitation numbers begin to grow and is therefore unable to explore a possible transition from a quantum regime with anti-bunched output photons to a classical regime with numerical Bogoliubov U/ J/ P-function mean-field Thursday, November 28, 13 Figure 1: Sketch of parameter regions where the discussed approaches provide accurate descriptions. With increasing driving strength Ω, the validity range of the numerical approach shrinks to larger nonlinearities whereas the validity ranges for the Bogoliubov and P-function mean-field approaches grow. For definitions of U , J, Ω and γ see equation (1). The validity boundaries of the various approaches are discussed further in figures 8, 9 and 12. uncorrelated output photons.
Here we present approaches that are capable of describing this transition. For a parameter regime with a weak nonlinearity in the central resonator we expand the intra cavity fields to linear order in the quantum fluctuations around their coherent parts. This Bogoliubov-type expansion provides a good approximation provided the density of quantum fluctuations around the coherent background is small compared to the excitation number in the latter. This regime is realized for weak nonlinearities in the central resonator.
In turn for strong nonlinearities but moderate photon tunneling rates between the resonators we combine an exact solution for a driven Kerr nonlinearity by Drummond and Walls [15] with a mean-field decoupling of the photon tunneling between resonators. As the solution is based on a description of the quantum state in terms of a P-function we will refer to it as P-function mean-field [16]. We find that this approach provides an accurate description in a large parameter range provided the photon tunneling is low enough such that the fields in the outer resonators do not deviate significantly from coherent ones. An overview of the discussed approaches and the parameter regimes where they provide an accurate description is sketched in figure 1.
The remainder of the paper is organized as follows, in section 2 we describe the setup and model we consider and in section 3 we revisit the numerical approach by Gerace et al. [14] to present results for time resolved correlation functions which so far have not been considered. In the next section 4 we discuss the approach based on a Bogoliubov expansion and the results it yields. In section 5 we then introduce our Here only the central conductors of the resonators are drawn. The three-resonator chip allows for a coupling J between two neighboring resonators [17]. A voltage-node in the center is the reason for an off-center implementation of the nonlinearity U that is capacitively coupled to the middle resonator. The two outer transmission line resonators are driven by coherent input sources with amplitudes Ω. Each transmission line resonator has an individual dissipation-rate, which depends on the resonator's coupling capacitors at both ends. We consider the case where the combined decay rate of the central resonator and the transmon is equal to the decay rate of the two outer resonators γ.
mean-field approach based on an exact solution for the central resonator. Finally section 6 presents conclusions and a discussion of the parameter regimes covered by each of the discussed descriptions.

Setup
The circuit QED Josephson interferometer we consider is a multicavity driven-dissipative system consisting of three coupled transmission line resonators where the two outer ones are coherently driven and the central transmission line resonator interacts with a superconducting qubit, see figure 2 for a sketch. We consider here a transmon as the qubit and assume a regime with strong interactions between the central resonator and this qubit so that polaritons, combined light matter excitations, become a useful description of the nonlinear central resonator. For large enough qubit-resonator coupling one can focus the description only on one excitation species in the central resonator [18]. Employing this approach, the Hamiltonian for the circuit QED Josephson interferometer can be written in a frame rotating at the frequency ω in of the input drives aŝ Hereâ k annihilates an excitation in transmission line resonator k. For the two outer resonators the excitations are photons whereas for the central resonators they are excitations of the considered polariton mode. ∆ k = ω k − ω in is the detuning between the input drives and resonance frequency of excitations in the k-th transmission line resonator. For simplicity, we only consider the on-resonance case with ∆ k = 0.
The parameter U describes the effect of the Kerr-nonlinearity due to the transmon in the central resonator and J is the tunneling rate between neighboring resonator sites, Ω is the drive amplitude, assumed to be equal for the two inputs, and ϕ describes the phase difference between them. At this point it is important to note that the nonlinearity of a transmon qubit is attractive (corresponding to U > 0 in our notation).
The circuit QED Josephson interferometer is an open quantum system due to experimental imperfections such as the finite reflectivity of the coupling capacitors at each end of the transmission line resonators and relaxation of the transmon. We consider here the case where excitations in all three transmission line resonators have the same dissipation rates. As excitations in the second resonator can also decay via transmon relaxation, equal dissipation rates can be achieved for example by choosing suitable capacitances at the ends of the transmission lines. Since the couplings to the environments are sufficiently small [19,20] we can make use of the master equation approach [21]. As a simplification we assume the reservoir state to be in vacuum, which is justified by the fact that circuit QED experiments are typically performed at a temperature in the millikelvin range [22]. These considerations lead to the master equatioṅ where the Liouvillian L is given by the dot denotes a time derivative and γ describes the excitation loss rate of an individual resonator site.
Since the excitation losses are continuously compensated by a coherent input drive, a dynamical equilibrium leading to a steady state will be established. We are here interested in properties of this steady state, such as the mean excitation number and its correlation functions, mainly for the central resonator.
In cavity or circuit quantum electrodynamics experiments the output fields of the investigated resonators are typically accessible by measurements. Hence in our device the quantity of interest in measurements is the output field of the central resonator. It is linked to the field in the resonator via the input-output relation [23], where a out and a in are the output and input fields and a res the field in the resonator which leaks out at a rate κ. In cavity and circuit QED the photons in the cavity usually couple to the photon modes in the output channel so that one can identify the cavity field a res with the annihilation operator of photons in the cavity. In the central cavity of our setup however, two polaritonic excitations build up due to the strong coupling to the qubit and the photon field in the central resonator can be expressed as a superposition of them, a res = ηa 2 + η ã 2 with coefficients η and η (Note that η can always be chosen to be positive and that |η| 2 + |η | 2 = 1). In the regime of strong coupling we are interested in, the frequencies of the polariton modes a 2 andã 2 differ strongly enough so that the input drives selectively only excite a 2 and we can neglect theã 2 excitations in the input-output relation since they only contribute vacuum output. Hence for calculating the output fields of the central resonator we can use the modified input-output relation, where γ = ηκ. For vacuum input fields, we thus obtain an out photon flux a † out a out = γ a † 2 a 2 for the central resonator.
The quantumness or non-classicality of light fields can be characterized by the statistics of their excitations. Here a quantity of central interest is the so-called (time-resolved) second order correlation function, here for resonator 2, The significance of the physical interpretation is to have a value for the probability to measure a second particle at time t + τ after a particle has been measured at time t and compare this to a coherent field where g (2) = 1. Therefore, if g (2) < 1 it is less likely to measure a second particle. In this case we speak of an anti-bunched light field that is necessarily non-classical, whereas for g (2) > 1 we have bunched light meaning a higher probability than for a coherent field. We note that the g (2) -function for the output fields is identical to the g (2) -function of intra-cavity fields provided the input fields are in vacuum.
Let us finally stress that although we refer to an implementation in circuit electrodynamics here, our calculations do in most parts not make use of details of this technology and are therefore applicable to other realizations as well. Yet we prefer to discuss them in the context of circuit QED as this technology offers excellent perspectives for realizing a coupling between resonators as we envision it here.

The Numerical Approach in the Weak Driving Regime
In the regime where the laser intensity is rather low, we are able to numerically solve the master equation directly. Here, an entire regime of rather high nonlinearities U/γ > 1, which lessen the probability of high intra-cavity excitations, can be covered for arbitrary values of the tunneling rate J. We introduce a cut-off  For the steady states we consider, g (2) (τ ) as given in equation (6) is independent of the time t and only depends on the time delay τ . In the weak driving regime the results for no delay τ = 0 have already been covered in [14] and we therefore focus for this regime on the time-resolved second order correlation function for τ = 0.
Applying the quantum regression theorem (QRT) [24] we find the results presented in Fig. 3. Here, as expected, the g (2) (τ )-function converges to 1 for a long time delay τ → ∞. The amplitudes of the flucutations around g (2) (τ ) = 1 increase with increasing nonlinearities and increasing tunneling rates J, however at τ = 0 the system is more excitation anti-bunched if the tunneling strength is weaker and the nonlinearity higher.
With the used approach of directly solving the master-equation numerically in a truncated Hilbert-space one can cover a regime for arbitrary values of J if the nonlinearity fulfills U/γ > 1 in the weak driving regime of |Ω|/γ < 2. Nevertheless, this method is unable to show a clear crossover between anti-bunched and uncorrelated excitations because the weak drive and high nonlinearity restrict it to low numbers of excitations, implying anti-bunched statistics. Therefore we will present a method to extend the investigation of this set-up to the strongly driven regime, motivated by the possibility of high intra-resonator excitations.

The Quantum Fluctuations Approach in the Strong Driving Regime
For higher laser intensities, the number of excitations increases rapidly. Since the two outer transmission line resonators do not directly couple to qubits, their only deviation from a coherent state stems from the coupling to the central resonator, which features a nonlinearity due to its coupling to a transmon qubit.
Therefore it is for strong inputs no longer computationally feasible to numerically solve the full master equation. Assuming a weak nonlinearity, we thus use a generalized Bogoliubov approach to non-equilibrium steady states [12], which linearizes the equations for the quantum fluctuations around the classical fields and can be understood as an approximation for low nonlinearities. We expand the creation and annihilation operators around the coherent state by writinĝ where the complex number α j = â j represents the coherent background and the operators δa j and δa † j describe the quantum fluctuations around it. Using equations (7) and (8) where L (n) [ρ] denotes the nth order of quantum fluctuations δa j in the Liouvillian L.

Steady State Coherent Background
To determine the coherent background in terms of the α j , we require that the first order Liouvillian vanishes for all quantum fluctuations, L (1) ! = 0. The result is equivalent to the solution for the steady state coherent background. With the above requirement we get a cubic equation for the steady state mean excitation number n 2 = |α 2 | 2 of the coherent background in the second cavity which has only one physically meaningful solution with n 2 ≥ 0: where η = 3 72γ 4 J 2 U 4 Ω 2 (cos(ϕ) + 1) + √ 3 γ 6 U 6 6912γ 2 J 4 U 2 Ω 4 cos 4 ϕ 2 The values of α 1 and α 3 can be determined from the solution (11) via the additional relations The solution (11) is plotted in Fig. 4. One can clearly see that the coherent background decreases for increasing nonlinearities. This feature appears since the frequencies of the input drives are resonant with the transition between zero-and single-excitation states in the resonators. Hence with increasing nonlinearity it becomes less probable that higher excitation states become populated as well. Moreover the population of the coherent part of the intra-resonator fields shows a maximum for specific values of J ≥ 0.
These maxima can be understood when writing the Hamiltonian (1) in terms of Bloch modes [12], which have frequencies that depend on the tunneling rate J. Thus with increasing J a balance establishes between the higher detuning between Bloch modes and input drives leading to less excitations and an increased tunneling leading to more excitations in the central cavity. Here, it can already be seen that the number of polaritons in the second resonator exceeds the computational feasibility of conventional numerical methods to calculate the 3-site system. but changes the severeness of the interference effect in agreement with the previously visualized results in Fig. 4. We further find (not plotted) that for an increase in driving strength the oscillation due to the phase-difference ϕ is steeper, changes the value of J for which the background occupation number has a maximum and decreases the incline around it. Also the sharpness of the peak of maximal background occupation is less pronounced for higher nonlinearities.

Steady State Quantum Fluctuations
We are mainly interested in the statistics of output photons from the central resonator. These are described by the mean excitation number n 2 and the second-order correlation function g 2 . To calculate the mean excitation number in the steady state, we can write in the expansion approach, where we have dropped the site index as we are only focusing on resonator 2. The second order correlation function has been introduced in Eq. (6). The case with zero time-delay τ = 0 is important in order to determine whether the output light field exhibits Poissonian or sub-Poissonian statistics. We expand the g (2) -function to get Therefore, in order to calculate the important statistical properties of our system we need to combine the previously derived results of the coherent background with the mean correlation values of the quantum fluctuations. We thus solve the master equation (2) to second order in the quantum fluctuations for ρ, using L (1) = 0 and plugging the solution (11) into L (2) . For the considered case of ∆ = 0, the eigenvalues of L (2) all have real parts −γ/2 so that the solutions are stable. Neglecting higher order terms in the master equation, i.e. L (3),(4) ≈ 0, is a good approximation provided that quantum fluctuations are small compared to the coherent background, To analyze the validity of the solutions we find, we thus check whether they fulfill this property.
For the central transmission line resonator, the solution for the mean quantum fluctuation number reads Considering the fact that we neglected the terms of third and fourth order in quantum fluctuations, which contain additional pre-factors of the nonlinearity U , we come to the conclusion that our approach can also be considered an expansion in U . As a consequence, we can deduce the results represent a good approximation in a regime of weak nonlinearities U . Moreover, the approach becomes much better for higher laser-input intensities because the coherent background excitation number in Eq.(11) scales much stronger with the laser input amplitude Ω than the number of quantum fluctuations in Eq. (19).
In contrast to the mean excitation number in the second resonator, the quantum fluctuations do play a significant role for the g (2) -function, see Eq. (16). The second order correlation function of the central resonator in the strong driving regime is plotted in Fig. 6. Depending on the tunneling rate J we observe smooth transitions from a coherent to an anti-bunched excitation field as the nonlinearity is increased beyond specific values. Fig. 7 shows that the oscillations in the excitation number due to interference effects can be observed in the g (2) -function as well.

Limitations of the Bogoliubov Approach
The major limitation of the Bogoliubov approach is obviously the validity of the assumption (18). In order to determine the regime where this is fulfilled, we plot the ratio δa † 1, 2, 3). For parameters where this ratio is sufficiently small, the approximation works fairly well. The ratio δa † 2 δa 2 / |α 2 | 2 is plotted in Fig. 8, whereas Fig. 9 depicts the same ratio for the two outer resonator sites, i.e.      As a conclusion, we can safely assume that our approximation works well in a regime for strong input drives Ω/γ > 5 and low nonlinearities U/γ < 0.2 for arbitrary values of the tunneling rate. For larger drives, such as Ω/γ > 10 the range of validity extends to larger U/γ. Thus we are motivated to find access to the regime of larger nonlinearities U/γ because, as is evident from Fig. 8, the assumption δa † i δa i |α i | 2 breaks down in that regime, primarily in the second resonator. We present an approach to this parameter regime in the next section.

P-Function Mean-Field Approach to the Low Tunnelling Regime
As one can see from figures 8 and 9, the main limitation of a Bogoliubov approximation is the ratio of the mean quantum fluctuation number to the mean excitation number in the second cavity which significantly increases for large nonlinearities U and weak drives Ω. Hence if we are able to take higher quantum fluctuations in resonator 2 into account the results should become more exact and be applicable to a larger parameter range. On the other hand the deviation from a coherent state for the first and third transmission line resonator is extremely small because they do not contain any nonlinearity themselves and the nonlinearity of the center cavity only slightly affects the two outer cavities. Thus, if we are able to find a way to solve the second transmission line resonator's mean values exactly, our results should become better. This will be the ansatz of our mean-field approach where we ignore quantum fluctuations in the two outer resonators and substituteâ 1,3 → α 1,3 . With this substitution, the Hamiltonian for resonator 2 takes the formĤ which can be written in the form of a coherently driven Kerr nonlinearity [15], by introducing the driveΩ = −2J (α 1 + α 3 ). To determine α 1 and α 3 self consistently, we use equations (13) and (14) to obtain,Ω in the ∆ = 0 case. Our approach can thus be understood as a mean-field decoupling of the three cavities.
With the consistency condition (22), the Hamiltonian (21) becomes a single site model which can be solved exactly using a P-function based method introduced by Drummond and Walls [15].
What we have done is basically a trade-off: We allow for a more crude approximation of the quantum tunneling effects in the two outer cavities by completely neglecting their quantum fluctuations. The latter are however expected to be very small, see Fig. 9. In turn we gain an exact solution for the resulting one-site problem, giving a better approximation of the nonlinear effects in the second transmission line resonator. In contrast to the Bogoliubov expansion that has been performed in section 4, we here do not expand in U but in the tunneling rate J [25]. With the approach introduced here we thus cover a regime with arbitrary values for the nonlinearity U but only for moderate values of J. The solution for the normal-ordered mean values of the middle transmission line resonator in the ∆ = 0 case reads where is the generalized Gauss hypergeometric series, defined in terms of the Γ-function. Note thatΩ is a function of â 2 , see Eq. (22). The starting point of our approach is thus to find the solution of the expectation value â 2 from Eq. (23) which, after inserting equation (22), becomes a nonlinear algebraic equation for â 2 , The result for â 2 can be put into Eq. (22) allowing for an expression of the "new" driveΩ. Based on this solution the correlation values for all higher normal-ordered momenta in Eq. (23) can be subsequently derived in a self-consistent manner. Fig. 10 shows the behaviour of the second order correlation function for excitations in the second cavity depending on the external experimental parameters. As expected, for low nonlinearities the state in the middle transmission line resonator is close to a coherent state. An increase of the nonlinearity leads to a state becoming more and more anti-bunched. In contrast to the low excitation regime [14], we here find that g (2) (0) does not monotonously approach unity as the tunneling J is increased. Instead anti-bunching can also become more pronounced. We attribute this non-monotonous behavior to the competition of two effects. On the one hand, an increased tunneling rate means that more photons from the coherent fields in the outer resonators flow into the central resonator and weaken the anti-bunching. On the other hand, an increased tunneling rate detunes the normal modes of the three-cavity system from the laser drives so that the outer cavities are less populated. If this second effect starts to dominate over the previous one, the 2 -function (b) of the central resonator in the strong driving regime, Ω/γ = 5, as a fundtion of J/γ and U/γ. The plots show the results gained from equation (23). Due to the previously explained nature of our approach being a good approximation for lower tunneling strength, we are able to cover an additional regime compared to the Bogoliubov expansion, which broke down for higher values of the nonlinearity. (Compare with Fig. 8) photon flow into the central resonator decreases and its excitations show more tendency to anti-bunch.
This signature is also visible, albeit less pronounced, in Fig. 8.
In Fig. 11 we investigate the influence of the phase-difference ϕ on the excitation statistics of the central transmission line resonator. The mean excitation number fluctuates with the phase-difference ϕ, exhibiting maxima at ϕ = 2nπ and minima at ϕ = 2(n + 1)π, where n is an integer, due to constructive and destructive interference of the input laser drives respectively. An increase in the tunneling rate from J/γ = 0.1 to J/γ = 0.5 enables a higher number of excitations to tunnel into, and therefore excite, the central transmission line resonator, leading to a higher number of excitations altogether. The g (2) -function also shows oscillations, however, their behavior is very different for different tunneling rates. Whereas maxima appear for both tunneling rates at a phase-difference ϕ = 2nπ, the second order correlation function for J/γ = 0.1 shows additional local maxima at ϕ = (2n + 1)π. These may be attributed to the low excitation number at these points which induce a tendency for g (2) -functions to increase.

Comparison to the numerical solution
Finally we quantify the parameter range where the P-function based mean-field approach we introduce here is in good agreement with a full numerical solution for regimes where the latter can be obtained. Fig.   12a compares the mean-field solution (red) and full numerical solutions (dashed lines) for Ω/γ = 1.5 and U/γ = 2 as the tunneling rate J is increased. Whereas the numerical solutions for cut-offs at four (blue)  Figure 11: Dependence of the second resonator's excitation statistics on the phase difference ϕ of the input laser drives for Ω/γ = 5, U/γ = 0.5 and J/γ = 0.1 (blue) or Ω/γ = 10, U/γ = 2 and J/γ = 0.5 (black) as obtained from equation (23). a: Mean excitation number and b: Second order correlation function. The oscillation in the mean excitation number can be attributed to the interference of the two input laser fields. and five (green) excitations agree very well indicating that the cut-offs are high enough, one can see deviations between the numerical and the mean-field solutions as J increases. These should be interpreted as a failure of the mean-field approach for larger J. Fig. 12b in turn compares the mean-field solution (red) and full numerical solutions (dashed lines) for J/γ = 0.2 and U/γ = 0.5 as the drive strength Ω is increased. Here the numerical solutions converge to the mean-field solution as the cut-off number is increased, indicating that the mean-field approach is a very accurate approximation throughout. The agreements and disagreements displayed in Figs. 12a and b are in agreement with the intuition that enhancing the coherent drives at the outer cavities should make the fields in those cavities more coherent and hence the mean-field approximation more accurate whereas an increase in the tunneling rate enhances the influence of the nonlinearity on the outer cavities and hence drives them away from a coherent state which renders the mean-field approach inaccurate.

Conclusion and Discussion of Accessible Regimes
In summary we have complemented the discussion of a quantum optical Josephson interferometer with two methods. One is a linearization or Bogoliubov expansion of the intra-cavity fields around their coherent parts. The second method, which we have introduced here, combines an exact solution for the central resonator with a mean-field decoupling of the tunneling terms to the outer resonators.
With the proposed methods the circuit QED Josephson interferometer can be described for a significant range of all experimentally adjustable parameters, such as the input drive Ω, the nonlinearity U and the tunneling strength J for a given dissipation rate γ. If one is faced with low laser intensities and thus low  Figure 12: Comparison of the results for mean polariton number in the central cavity obtained by the Pfunction mean-field method, red line, and by the conventional numerical method, dashed lines, with cut-offs at three (black), four (blue) and five (green) excitations. a: Deviation due to an increase in the tunneling rate for Ω/γ = 1.5 and U/γ = 2. b: Deviations due to an increase in driving strength for J/γ = 0.2 and U/γ = 0.5.
intra-cavity excitations, the solution can be obtained by solving the master-equation numerically via an excitation cut-off in the matrix representation of the ladder operators as has been shown in section 3. In a strong driving regime with intermediate to high tunneling rates but low nonlinearities a good approximation to the excitation output statistics can be given by a Bogoliubov expansion in the quantum fluctuations around the coherent background parts of the intra-cavity fields, which we presented in section 4. However, for a strong input drive, low tunneling rates and intermediate to high values of the nonlinearity the Bogoliubov expansion breaks down. Here we derived a method in section 5, which we termed P-function mean-field approach, that works extremely well in said regime. The only regime for which a good approximation could not be gained was one where the drive is strong enough to produce high intra-cavity excitations, but the tunneling rate and nonlinearity are still relatively high. As a result they are responsible for non-neglectible quantum effects in comparison to the drive or to each other and therefore neither an expansion in U nor in J would be valid. A summary of the validity ranges for each of the presented methods is sketched in figure 1.

Author's contributions
RJ performed the calculations and analysis, MH conceived and supervised the project, both authors discussed the results and wrote the manuscript.