Josephson photonics with a two-mode superconducting circuit

We analyze the quantum dynamics of two electromagnetic oscillators coupled in series to a voltage biased Josephson junction. When the applied voltage leads to a Josephson frequency across the junction which matches the sum of the two mode frequencies, tunneling Cooper pairs excite photons in both modes simultaneously leading to far-from-equilibrium states. These states display highly non-classical features including strong anti-bunching, violation of Cauchy-Schwartz inequalities, and number squeezing. The regimes of low and high photon occupancies allow for analytical results which are supported by a full numerical treatment. The impact of asymmetries between the two modes is explored, revealing a pronounced enhancement of number squeezing when the modes are damped at different rates.


I. INTRODUCTION
It has long been known that the current flowing through a voltage-biased mesoscopic conductor can provide an extremely sensitive probe of its electromagnetic environment [1][2][3][4] . The current-voltage characteristics of a tunnel junction placed in series with a transmission line resonator is a particularly well-studied case 1,2,5 . The transmission line resonator contains a series of welldefined harmonic modes whose presence opens up inelastic current channels leading to characteristic features in the dc current flowing through the junction 5 . The advent of high-Q superconducting resonators whose quantum state can be measured with great precision 6 together with the development of hybrid devices which couple nonmetallic conductors to resonators 7,8 , has led to a renewed interest in the interaction between tunneling electrons or Cooper pairs and harmonic modes. Whilst earlier experiments 5,9 on mesoscopic conductors coupled to electromagnetic resonators focussed on how the harmonic modes affect the current in a regime where the modes themselves are close to thermal equilibrium, more recent experimental [10][11][12][13] and theoretical work [14][15][16][17][18][19][20][21][22][23][24] has begun to investigate how the current influences the resonator state and to explore the dynamics of systems where the resonator is far from thermal equilibrium.
For a Josephson junction which is biased with a subgap voltage, V , the relationship between the dc current and the energy pumped into the electromagnetic environment is particularly simple as all of the energy associated with a tunneling Cooper pair must be absorbed by the environment, 11 . When the Josephson junction is placed in series with a transmission line resonator a dc current is expected when the ac Josephson frequency ω J = 2eV / matches one or more of the mode frequencies in the transmission line. Experiments using low-Q resonators 5,11 have demonstrated that when the individual harmonic modes remain close to thermal equilibrium, they lead to well-defined peaks in the dc current whose heights and widths can be calculated using perturbation theory. In contrast, a high-Q resonator can be excited to far-from- equilibrium states containing many photons 13 which are predicted to display intriguing non-classical features such as number squeezing 20,21 . This new field of Josephson photonics combines typical processes known from quantum optical set-ups with those known from charge transfer physics in highly versatile devices.
In this article we consider a voltage-biased superconducting junction whose ac Josephson frequency is tuned to excite two electromagnetic modes simultaneously (see Fig. 1). Signatures of such processes have been observed in the dc current flowing through Josephson junctions coupled to low-Q resonators and can also be understood within a perturbative approach as the modes remain close to thermal equilibrium. While we address this domain as well, our main focus here lies in the regime where the power transferred to the resonator modes is sufficient to drive them into far-from-equilibrium states while still displaying strong quantum properties. Note that the system we consider here differs from those used in recent experiments to produce photon pairs 25,26 in that the energy comes from a dc voltage.
Starting from a simple model Hamiltonian which describes the effect of the Cooper pairs on the oscillators through a highly non-linear ac drive at the Josephson frequency, we use a rotating wave approximation to derive an effective time-independent Hamiltonian which we use to analyse the quantum dynamics of the oscillators. Although the full behavior of the system can only be uncov-ered by numerical solutions of the quantum master equation, we find that approximate analytical descriptions are available in the two regimes of low and high photon occupancy. In the former one a perturbative treatment in the Josephson energy applies while in the latter explicit results are obtained by linearizing about the classical fixed points which provide a faithful description of the quantum dynamics when the zero-point fluctuations of the oscillators are small.
The excitation of the two oscillators shows a clear threshold as a function of the Cooper pair pumping rate. Earlier work, which investigated the quantum dynamics of a single mode 20-22 driven by a voltage-biased Josephson junction, showed that non-classical features in the state of the oscillator such as number squeezing (sub-Poissonian photon statistics) occur very generally. For the two-mode system, we also find significant numbersqueezing occurs in the states of the individual oscillators, especially in the above-threshold regime where the oscillators are strongly excited. Interestingly, when the damping rates of the oscillators are very unequal, the lessdamped oscillator displays much stronger strong number squeezing than is ever found for a single-oscillator system. Provided that the quantum zero-point fluctuations are not too small, the number squeezing is strong enough to lead to negative regions in the Wigner function.
This work is organised as follows. We introduce our theoretical model in Sec. II and analyse its low photon limit in Sec. III and its semi-classical dynamics in Sec. IV. Sections V and VI explore the quantum dynamics of the system in the below and above threshold regimes, respectively. Finally, Sec. VII contains a discussion and the conclusions. The Appendix contains further details on some of the calculations described in the main text.

II. MODEL SYSTEM
We consider a system consisting of a Josephson junction in series with two LC oscillators, A and B with angular frequencies ω a and ω b across which a voltage V is applied (see Fig. 1). The two oscillators could both be modes of a single superconducting resonator in which a Josephson junction is embedded between the ground plane and center conductor 13,20,27,28 (See Ref. 20 for a detailed derivation of the Hamiltonian for this case), but the system could also be realized using modes of two different electrical resonators 5 . The effective Hamiltonian of the system takes the form where E J is the Josephson energy of the junction, a and b are the lowering operators of the oscillators with frequencies ω a and ω b respectively, and ω J = 2eV / . The parameters ∆ a(b) quantify the strength of the zero-point fluctuations of the oscillators, Here we analyze the case where the system is operated close to the resonance that occurs when the voltage energy lost by a single Cooper pair traversing the circuit matches the energy required to simultaneously create one photon in each of the LC oscillators, ω J = 2eV / = ω a + ω b . We assume that the modes are not degenerate so that ω a = ω b . This means that the resonance at ω J = ω a + ω b does not compete with processes in which two photons are absorbed by just one of the modes.
We examine the behavior of the system as a function of the Josephson energy which describes the strength of the Cooper pair tunneling. The value of E J can be thought of like a pumping rate for the oscillators: as it is increased the oscillators will be more strongly driven, become more strongly excited and behave more non-linearly. In practice E J can be varied in an effective single-junction by forming two junctions in parallel and applying a tunable flux in the SQUID loop that they form 20,29 .
The strengths of the quantum fluctuations parameterised by ∆ a , ∆ b , also play a very interesting role in determining the dynamics of the system and we will examine how the behavior is modified when they are varied. For systems where a Josephson junction is embedded in a superconducting resonator designed to have a very high-Q the quantum fluctuations will typically be very small ∆ a(b) ≪ 1 . However, significantly stronger quantum fluctuations have very recently been engineered in low-Q resonators coupled to tunnel junctions 30 and it may be possible to combine stronger quantum fluctuations with higher Q values in the future.

A. Rotating wave approximation
The explicit time-dependence in the Hamiltonian complicates the analysis of the corresponding dynamics significantly. However, close to the resonance we are interested in, ω J ≃ ω a + ω b , only some of the terms will play an important role and these can be picked out by a rotating wave approximation (RWA).
We proceed following the approach in Refs. 20-22. We move to a rotating frame, applying a unitary transformation of the form U (t) = e iωaa † at e iω b b † bt where we definẽ ω a +ω b = ω J , and make a RWA in which we neglect all of the rapidly oscillating terms in the rotating frame. The resulting effective Hamiltonian takes the form, where the colons imply normal ordering of the operators, For sufficiently low photon numbers (such that 2∆ a a † a , 2∆ b b † b ≪ 1) we can expand the Bessel functions in Eq. (2) to lowest order. In this limit the system reduces to a non-degenerate parametric amplifier 31

B. Quantum master equation
The two oscillators are assumed to be weakly damped at rates γ a and γ b which in general will not be the same. We therefore assume that the quantum master equation of the system takes the standard quantum optical form in the T = 0 limit 31 dρ dτ where we adopt dimensionless units of time . In an actual experimental realization of the JJoscillators system in Fig. 1, the damping of the oscillators (due to photon decay from the resonators) is not the only source of dissipation. Indeed, the existence and impact of local voltage fluctuations at the JJ can be seen in the broadening of the spectrum of emitted microwave radiation 11,21 . The existence of such fluctuations necessitates including explicitly an extra degree of freedom for the number of Cooper pairs N transported across the junction in the model. In the effective Hamiltonian, Eq. (2), the a † b † + ab term then gets replaced by e iη a † b † + e −iη ab , where e ±iη = N |N N ± 1|. Local voltage fluctuations are included by an additional dissipator in (4) which in the simplest version takes the form 21 describes, how to treat the corresponding quantum master equation in the extended JJ-resonator space.
However, it turns out that only certain observables sensitively depend on the strength of these fluctuations, characterized by γ J , for example the spectral broadening. For other observables, such as the photon occupation and photonic correlation functions that are of relevance for this work, the impact of local voltage fluctuations is basically negligible since experimentally one typically has γ J ≪ γ a,b (see for example Ref. 11). Then, formally, the Hamiltonian (2) is regained by putting γ J = 0 so that the phase operators e ±iη simply appear as phase factors which can be removed via the gauge transformation e iη/2 a † , e iη/2 b † → a † , b † . Note that this reflects a phase invariance of the RWA Hamiltonian (2).

C. Relevant observables
The basic structure of the RWA Hamiltonian [Eq. (2)] in which photons are always created (or destroyed) jointly in the two oscillators and the linear damping that we assumed in formulating the master equation lead to a simple connection between the occupation numbers of the two modes n a(b) = a † a(b † b) and the average dc current, I dc , flowing through the junction that can be obtained from an energy balance argument without the need to work with a current operator. Since each Cooper-pair that contributes to the dc current must create exactly one additional photon in each of the oscillators, the requirement that the energy gain and loss rates balance tells us that where in this case we have returned to dimension-full units.
The quantum nature of the photonic states in the oscillators is captured by photon correlation functions such as , g ab (0) = a † ab † b n a n b (6) and the Fano factors Whilst these two types of correlation functions are closely related to each other, they are nevertheless useful to characterize the photonic states in opposite regimes of parameter space. In the regime of weak driving and low photon occupation deviations from the case of a driven harmonic oscillator are best seen in the g (2) functions. Namely, with increasing driving amplitudeẼ J , the photon distributions for the number states in the cavities evolve from Poissonian distributions with almost empty cavities towards distributions peaked around finite mean occupations n a , n b . In this case the g (2) (0) functions (6) sensitively indicate deviations from the linear regime g (2) aa(bb) (0) ≡ 1 with g (2) ab (0) = 0 capturing growing cavitycavity correlations. In the opposite regime of strong driving, nonlinearities may substantially influence the widths of the peaks for photon occupations (energy fluctuations) as properly measured in the Fano-factors (7).
In the following, we will first focus on the regime of low photon occupancy, where charge transfer through the Josephson junction occurs sequentially (Coulomb blockade regime) and analytical results can be obtained via a perturbative treatment in the drive amplitude E J . In the opposite domain of large photons numbers in the cavities, a semi-classical type of approximation applies, though, with substantial quantum fluctuations. In both domains, the magnitude of the parameters ∆ a , ∆ b , i.e. the strength of the ground state fluctuations of the respective cavities, plays a decisive role: they rule the impact of nonlinearities in the former regime and control the impact of quantum effects in the latter one. Corresponding findings will be supported by numerical solutions of the stationary states according to (2) and (4).

III. FEW-PHOTON LIMIT
The physics of the system described by the Hamiltonian (2) and the master equation (4) is at its simplest when it is driven so weakly that excitations in the resonators will relax to equilibrium well before a new excitation occurs. In that regime, very few photons, n a/b ≪ 1, reside in the resonators on average. Transport across the junction in turn is in the (dynamical) Coulomb-blockade regime, where subsequent Cooper-pair tunneling events occur uncorrelated with some tunneling rate. While the charge flows uncorrelated, the photons exhibit correlations already at the weakest driving. Now, for the present set-up one derives from the full RWA Hamiltonian )/2 and where n b follows by replacing r → 1/r. In the lowest order in the driving strength this reduces to with the superscript indicating the leading order in E 2 J and with n (0) b again following from r → 1/r.
For the correlations we focus on the symmetric case γ a = γ b at resonance so that n a = n b = n. Then, based on the full Hamiltonian (2) one can show the general relation bb (0) which implies g (2) ab (0) = 1 2n bb (0) (11) with n as given in (8). Now, working to order E 4 J , one finds 8 .
(12) Two types of correlations are encoded in the above g (2) (0) functions. The most obvious ones stem from the common excitation process of photons in the two resonators. They are therefore already present in the parametric amplifier limit of the Hamiltonian (3) and well understood for that case, see e.g. Ref. 32. A convenient tool to characterise them is the noise reduction factor 26 However, the perfect correlation of the excitation process leads to perfectly correlated occupations in the oscillators with a noise reduction factor NRF = 0 only for the undamped case γ a = γ b = 0. For any finite photon lifetimes in the cavities, the decay out of the two cavities occurs uncorrelated which according to (11) always implies in the stationary state and for the symmetric situation NRF = 1/2. Further correlations in the light field are caused by the back-action of the resonator occupations on the photon creation processes. Generally speaking, the existence of photonic excitations in the resonators can either increase the probability of further excitations, similar to a stimulated emission effect, or it can hinder further excitations. Formally, these effects are encoded in the transition matrix elements of the RWA-Hamiltonian (2) between neighboring oscillator states, where the nonlinearities of the Bessel functions enter. If charge quantization of the Cooper-pair current is significant, the parameters ∆ a/b become large, so that the nonlinearities already appear at the few photon level. For the case of a single resonator, it was shown in Ref. 21 that ∆ 2 = 2 can completely suppress transitions to higher occupations and reduces the resonator effectively to a two-level system, thus operating as a perfect single photon source. The behavior of the correlation functions in the two-mode case is shown in Fig. 2. While a non-zero g (2) aa (0) requires oscillator a to be populated up to the second excited state by two successive photons, this need not be the case for oscillator b as it can relax before the second photon arrives. Consequently, as seen in (12), g (2) aa (0) = 0 at ∆ 2 a = 2, but not at ∆ 2 b = 2. The general result (11) also reveals that the classical Cauchy-Schwartz inequality for photon intensities is always violated in the quantum case, i.e., g (2) aa (0) g (2) bb (0) ≤ g (2) ab (0) .
Namely, introducing the parameter ǫ = g (2) bb (0)/g (2) aa (0) the violation of the inequality requires [−g (2) aa (0)](1 − √ ǫ) 2 ≤ 1/n which always applies since g (2) aa (0), n ≥ 0. Accordingly, emission of photons from the cavities occurs in a correlated way for all driving strengths and photon occupations. In the next Section we ascribe to the individual photon states in the cavities respective amplitudes (energies) and phases. One then sees that these states are correlated through their phases due to the simultaneous creation process in the transfer of a single Cooper pair.

IV. SEMI-CLASSICAL DYNAMICS
For large photon occupancy in the cavities a semiclassical type of approximation applies. The simplest semiclassical description of the dynamics of the system is obtained from the equations of motion for a and b which follow from Eq. (4), making the replacements a = α, b = β and treating expectation values of products of operators as products of expectation values. Hence we findα whereδ (a,b) = δ (a,b) / √ γ a γ b . Obtained in this way, the factors of e (∆ 2 a +∆ 2 b )/2 embodied in E c J that appear in these equations are accidental: they would not be present if we had instead chosen to use a symmetric ordering for the operators when deriving the Hamiltonian. However, Eqs. (15) and (16) would also arise from a simple-minded ansatz in which we assumed that the density operator of the system is just a product of the coherent states ρ(t) = |α(t) α(t)| ⊗ |β(t) β(t)|, in this approximation the factors of e (∆ 2 a +∆ 2 b )/2 would arise naturally. Using amplitude-phase coordinates for the two oscillators, α = Ae −iφa and β = Be −iφ b , and introducing the total and relative phase variables ξ ± = φ a ± φ b , Eqs. (15) and (16) take the forṁ where we used the Bessel function identity, J 2 (z) + J 0 (z) = 2J 1 (z)/z, and have defined δ (±) =δ (a) ±δ (b) .
Further,  A, B). The behavior of the system is determined by the fixed points of the amplitudes A 0 , B 0 and the total phase ξ + 0 . Since the relative phase does not appear on the righthand side of any of these equations its fixed point value is arbitrary. For simplicity, we concentrate on the on-resonance case δ (a) = δ (b) = 0 in our analysis.
The amplitude equations lead to the fixed point con- The second equality in Eq. 23 leads to the energy balance condition B 0 = rA 0 . From the equation for ξ + , we see that fixed points arise when either cos ξ + 0 = 0 or This latter condition is independent of E J and hence leads to a locking of the amplitudes at particular values as a function of E J , something which is an important characteristic of the dynamics in the single-oscillator system 20 . For symmetric oscillators (r = 1 and ∆ a = ∆ b ) F + = 0 implies J ′ 1 (z) = 0 with z = 2∆ a A 0 = 2∆ b B 0 which has a first solution at z = 1.841 20 .
Thus we identify three possible fixed points for the system: a zero-amplitude one, one given by the conditions cos ξ + = 0 and (from Eq. (23)) and a third solution for which the amplitudes lock to values where Eq. (24) is satisfied (together with the condition B = rA) and the total phase is be given by Eq.
. We can look for small amplitude solutions to Eq. 25 (∆ b rA 0 , ∆ a A 0 ≪ 1) by expanding the Bessel functions and retaining the lowest order terms in A 0 , Thus we see that a non-zero amplitude solution only exists for E J > E c J . Thus E c J has a simple physical interpretation: it is the value of E J at which the oscillators reach the threshold for non-zero amplitude oscillations.
Taking into account the stability of the fixed points, we find that as E J is increased from zero the amplitudes remain zero until the system reaches threshold at E J = E c J , after which the amplitudes grow smoothly according to Eq. (25) with the global phase locked to ξ + 0 = π/2. For a sufficiently large E J , which we define as E c2 J , a bifurcation occurs as the amplitudes become large enough to satisfy Eq. (24) and the amplitudes then lock, becoming independent of E J .
In the next two sections we will examine the quantum dynamics of the system in the below and above threshold regimes.

V. SUB-THRESHOLD DYNAMICS
In the sub-threshold regime (E J < E c J ) the semiclassical fixed points have zero amplitude (A = B = 0). In this case we can gain some insight into the behavior of the system by approximating the Hamiltonian of the system by its lowest order terms, i.e. setting H RWA = H When this approximation is made the Hamiltonian is quadratic and the equations of motion for the moments take a rather simple form. Solving these equations, we find in the steady-state We note in passing that the result for n a reduces to the one derived in (9) in leading order in E J /E c J . Simplified in this way, the linearized description leads to a Gaussian steady-state Wigner function which takes the form 36,37 30) where C = (n a + 1/2)(n b + 1/2) − |µ| 2 and µ * = − ab . This is a mixed state which combines two-mode squeezing and thermal-like fluctuations 37 . The Wigner function of the individual oscillators is obtained by integrating over the phase space of the other one leading in either case to a thermal distribution. Thus for oscillator a, for example, we have W a (α) = 1 π(n a + 1/2) exp − |α| 2 (n a + 1/2) . The full behavior of the average energy of oscillator a, n a , obtained by solving the master equation numerically 33 , is shown in Fig. 3 for symmetric oscillators (r = 1, ∆ = ∆ a = ∆ b ). The divergence in n a which the linearized analysis predicts for E J → E c J [Eq. (27)] never occurs in the full quantum problem as higher order terms in the RWA Hamiltonian always saturate the energy gain. As ∆ is increased the saturation occurs at progressively lower values of the photon number whilst the range of E J /E c J values for which the linearized calculation is accurate becomes smaller and smaller.
The fluctuations in the energy of the oscillators, described by the Fano factors F a(b) (7) change rather more dramatically with ∆ a . The thermal Wigner function obtained from the linearized calculation [Eq. (31)] predicts the simple relationship between Fano factor and photon number associated with thermal states, F a(b) = n a(b) + 1, leading to growth in F a(b) as E J /E c J increases and again there is a divergence at threshold. For small values of ∆, the full quantum dynamics follows a similar pattern though with saturation in F a(b) at the threshold leading to a peak rather than a divergence. In contrast, for larger ∆ values the behavior is completely different: the value of F a(b) drops monotonically as E J /E c J is increased and its behavior contains no signature of the threshold at E c J . The change in the behavior of F a as ∆ a is increased is reminiscent of quantum optical systems like the laser 38 , which in the 'thermodynamic' limit of weak atom-photon couplings display clear thresholds (accompanied by a signature peak in the Fano factor) whose properties can be understood in terms of an analogy with classical phase transitions, but which for sufficiently strong couplings behave quite differently without clear signatures of a threshold 38,39 .

VI. DYNAMICS ABOVE THRESHOLD
Above threshold the oscillators become strongly excited though this does not mean that their states become classical. As in the case of the single-oscillator system 20 , strong number squeezing (marked by a Fano factor below unity) occurs even at large average occupation numbers. As in the sub-threshold regime, the behavior of the system in the limit of very small zero-point fluctuations, ∆ a , ∆ b ≪ 1, can be captured within an approximate description which linearizes about the semi-classical fixed points of the system, but for larger zero-point fluctuations numerical solution of the quantum master equation becomes essential. We start by exploring the general properties of the steady-states of the individual oscillators in the above-threshold regime for symmetric oscillators and the role played by the size of the zero-point fluctuations before going on to examine how asymmetry alters the behavior.

A. Symmetric Oscillators
For symmetric oscillators (r = 1, ∆ a = ∆ b = ∆) the steady-state properties of the two oscillators must be the same and there is a very simple scaling to the semiclassical fixed point amplitudes obtained in Sec. IV: the value of 2∆ a A 0 is a function of just E J /E c J , see (23). This scaling provides a convenient way of comparing the average oscillator occupation n = n a = n b (obtained by solving the master equation numerically) for different values of ∆ with the semi-classical prediction, as shown in Fig. 4a. We solved the master equation using standard numerical methods 33 ; for smaller values of ∆ we carried out quantum trajectory simulations, whilst for larger ∆ we were able to solve for the steady-state of the master equation directly because the state-space required was rather smaller. Indeed, the strong suppression in the magnitude of the oscillator occupation number as ∆ is increased (there is a reduction by a factor ∼ 100 in going from ∆ = 0.1 to ∆ = 0.6) is the most significant feature in Fig. 4a, which is captured by the 4∆ 2 scaling. Figure 4a also shows that the semi-classical amplitudes provides a very good description of the oscillator occupations for ∆ ≪ 1. For ∆ = 0.1 we see that there are small deviations from the semi-classical predictions which become apparent just above threshold and near the bifurcation that occurs at E c2 J = 2.5E c J . As the size of the zero-point fluctuations is increased, these small devia- tions grow much larger and spread out over a much wider range of E J /E c J values. Nevertheless, the semi-classical amplitude continues to provide a useful estimate of the full quantum results even for ∆ = 0.6.
We now turn to the fluctuations in the occupation numbers of the oscillators, described by the single mode Fano factors, F = F a = F b . The value of F decreases progressively the further above threshold we go as shown in Fig. 4b. For very small ∆, F is strongly elevated close to threshold (the other side of the peak in F seen below threshold), but decreases rapidly with increasing E J /E c J leading to substantial number state squeezing with F ∼ 0.5 before the bifurcation at E c2 J . For larger ∆ values there is no peak around threshold and F < 1 throughout though the lowest values are slightly larger than those obtained for very small ∆.
The simple semi-classical analysis in Sec. IV can be extended to describe fluctuations (at least to lowest order) by essentially adding a noise term to the equations of motion for the amplitudes, Eqs. (15) and (16), so that they become Langevin equations. Formally, such Langevin equations can be derived within the framework of an approximate semi-classical approach known as the truncated Wigner approximation, as we show in Appendix A. We again make the change to amplitudephase variables and then linearize about the fixed point values to obtain expressions for the amplitude fluctuations δA 2 = (A − A 0 ) 2 which can be related to the Fano factor in a simple way F a ≃ 4 δA 2 (details of the calculation are provided in Appendix A).
The comparison of the semi-classical and quantum calculations of the Fano factor shown in Fig. 4b shows that the semi-classical Fano factor, which is a function of E J /E c J alone in the symmetric case, can be thought of as giving the low-∆ limit. As ∆ is increased the deviations from the semi-classical value get stronger around threshold and the bifurcation at E c2 J = 2.5E c J as well as spreading over a wider range of E J /E c J in much the same way as for the oscillator occupation. Note that the semiclassical calculation predicts a Fano factor which tends to 0.5 as the system tends to the bifurcation, E J → E c2 J . This matches the lowest Fano factors found for the oneoscillator system which occurs as the system tends towards an above-threshold bifurcation at the 2-photon resonance 20 .

B. Asymmetric oscillators
We now consider what happens when the oscillators are no longer entirely symmetric. We start by considering the case where the zero-point fluctuations of the modes remain the same (∆ = ∆ a = ∆ b ), but the damping rates are different r = 1 and then go on to consider the general case where ∆ a = ∆ b and r = 1.
The effect of asymmetric damping on the average occupation numbers of the oscillator (shown in Fig. 5), is twofold with both effects following from the underlying semi-classical dynamics discussed in Sec. IV. Firstly, the bifurcation which occurs at E c2 J is pushed to larger values of E c J . Secondly, the average energies of the modes become unequal in proportion to the underlying asymmetry in the damping, n b = r 2 n a . Figure 6 shows the effect of asymmetric damping on the occupation number fluctuations for different values of ∆. What is striking here is that the fluctuations become asymmetric and the Fano factor becomes significantly lower than 0.5 in the less damped oscillator. The lowest values of F are achieved well-above threshold, close to the bifurcation at E c2 J for small-∆, though for larger ∆ values the minimum F is at a lower value of E J as the increase in F associated with the bifurcation starts to occur at progressively smaller values of E J /E c J as ∆ is increased. Above the bifurcation the value of F settles down to a steady, but rather higher value.
The semi-classical calculation predicts a minimum value of F ≃ 0.1 for the small-∆ limit when r = 1/3, substantially lower than any of the Fano factors predicted for the single-oscillator system 20 , and this value continues to decrease for smaller r. This suggests that the asymmetric two-oscillator system may provide a very effective route to preparing a particular mode in a strongly non-classical state at large photon numbers. As F → 0 the state of the oscillator must eventually become a pure Fock state and so one naturally expects to find negative features in the Wigner function for very small values of F . However, the presence of negative regions in a Wigner function is not simply a function of F , but also the average occupation number n : as one goes to larger average oscillator occupation numbers, smaller and smaller values of F are required to form negative regions.  Fig. 6). For ∆ = 0.6 there is strong evidence of negativity in the Wigner function whilst it is almost washed out for ∆ = 0.3 since although the Fano factors are very similar, the latter has a much higher average occupation number.
Finally, we examine the behavior in the regime where ∆ a = ∆ b . Figure 8 shows examples of the behavior of the occupation numbers and Fano factors of the two oscillators in this case. Interestingly for r = 1 whilst energy balance means that n a = n b , the fluctuations in the two modes are no longer the same. When r = 1 the occupation numbers of the two oscillators spilt according to the usual relation, n b = r 2 n a and the fluctuations become even more asymmetric. Indeed, the minimum values of the Fano factors, are lower than those in the cor- responding cases where ∆ a + ∆ b takes the same value, but ∆ a = ∆ b .

VII. DISCUSSION AND CONCLUSIONS
We have analyzed the quantum dynamics of two electromagnetic oscillators coupled to a voltage biased Josephson junction. We considered the case where the voltage across the junction was tuned so that the energy lost by a Cooper pair crossing the circuit matches the sum of the photon energies of the two oscillators. In this regime the oscillators are pumped by the flow of Cooper pairs and can become strongly excited. Using a rotating wave approximation, we derived an effective timeindependent Hamiltonian for the system and explored the behavior it gives rise to under a wide range of conditions using a mixture of numerical and analytic approaches to solve the master equation. We use a perturbative approach to obtain analytic results for the regime where the occupation of the oscillators is low while in the opposite regime of large occupation numbers a semi-classical approach provides an effective description.
The steady states of the oscillators display signatures of non-classical behavior over a very wide range of conditions with sub-Poissonian photon statistics found in both the low and high occupancy regimes. The strength of the zero-point fluctuations in the oscillators, ∆ a(b) , plays an important role: as these are increased the overall excitation level of the oscillators tends to move towards lower photon numbers whilst the signatures of non-classicality are enhanced. The ratio of the damping rates of the two cavities, described by r = r a /r b , also has an interesting effect on the behavior of the system. The photon numbers in the two oscillators are related in a simple way, n b = r 2 n a , as one would expect. However, the quantum fluctuations (e.g. measured by the Fano factors F a(b) ) also become unequal in the asymmetric case, r = 1. Indeed we find that the Fano factor in the less-damped oscillator can become low enough to lead to significant negative regions in the corresponding Wigner function.
Strong correlations between the two oscillators are to be expected in the regime we consider given the fact that the tunnelling Cooper-pairs excite photons in each of the two oscillators simultaneously. The violation of the classical Cauchy-Schwarz inequality for the photons in the two oscillators, g (2) ab , indicates that the corresponding two-mode states are non-classical. It would be natural to also investigate the entanglement between the two oscillators. However, this is complicated by the fact that in practice local voltage fluctuations, even when weak, would be expected to have a very strong influence on phase dependent correlation functions such as ab which can be important in determining the level of entanglement. This is in contrast to the observables such as photon occupation numbers and correlation functions which we have focussed on here which, as remarked in Sec. II B, are expected to be only very weakly affected. We plan to address the issue of inter-oscillator entanglement in a future work using a form of the master equation where the effects of voltage fluctuations are explicitly included 21 . where and a corresponding expression for h (b,a) . The noise terms obey the correlation functions Using Eq. (A5) we obtain the steady-state variances . (A13) Recalling that α and β are phase space variables of a Wigner function, we can connect these variances to quantum averages: A 2 = a † a + 1/2 and A 4 = (a † a) 2 + a † a + 1/2. For fixed points where A 0 ≫ 1, corrections of order A −2 0 can be neglected, leading to the simple result, and there is of course a corresponding relation for F b . The Langevin equation for δξ + takes the forṁ where η ξ + (t)η ξ + (t ′ ) = 2Dδ(t− t ′ ) with 2D = r/(4A 2 0 )+ 1/(4rB 2 0 ). Hence we find (δξ + ) 2 = D/F + (A 0 , B 0 ).
Note that as the system approaches the bifurcation at E J = E c2 J , F + (A 0 , B 0 ) → 0 implying that the total phase fluctuations within this linearized approach diverge.