Josephson Photonics with Simultaneous Resonances

Inelastic Cooper pair tunneling across a voltage-biased Josephson junction in series with one or more microwave cavities can generate photons via resonant processes in which the energy lost by the Cooper pair matches that of the photon(s) produced. We generalise previous theoretical treatments of such systems to analyse cases where two or more different photon generation processes are resonant simultaneously. We also explore in detail a specific case where generation of a single photon in one cavity mode is simultaneously resonant with the generation of two photons in a second mode. We find that the coexistence of the two resonances leads to effective couplings between the modes which in turn generate entanglement.


I. INTRODUCTION
Circuits in which voltage biased Josephson junctions (JJ) are combined with microwave cavities provide an ideal platform for exploring a wide range of microwave photonics. All of the voltage energy associated with tunneling Cooper pairs must be transferred into photons and the properties of JJ-cavity systems can be tuned over a wide range either in-situ or by design, [1][2][3][4][5][6] . Furthermore, the energy transferred by tunneling Cooper pairs into microwave modes can be tracked by monitoring either the resulting dc current or the microwaves leaking out of the circuit 1 . Recent experimental, 1-6 and theoretical work 7-18 has explored a wide range of ways in which JJ-cavity systems can be used to generate non-classical microwave states.
Energy exchange between charge carriers and microwaves in JJ-cavity systems is concentrated at resonances where the energy lost by a given Cooper-pair is commensurate with that of the photons in one or more microwave mode(s) [1][2][3] . Such resonances can be selected by simply tuning the voltage and are modelled theoretically using a rotating wave approximation (RWA) which leads to a convenient time-independent Hamiltonian for the system 9, 10,12,19 . The simplest resonances involve a single mode and can be exploited to provide a single photon source 4,5 , although higher order resonances in which two or more photons are generated within a particular mode also occur 9,18,20 .
Resonances involving two modes (realised, e.g., within the same cavity or in two separate cavities in series with the JJ) can be used to produce entangled photons, via processes in which photons in both are generated simultaneously via a single tunneling process 4,6 . The effective coupling between modes generated by the JJ also supports resonances where Cooper pair tunneling is accompanied by an exchange of photons between modes, processes which could be exploited to engineer efficient heat engines 21 .
Despite the very wide range of possibilities offered by JJ-cavity systems, so far attention has generally focused only on cases where a single photon generation/exchange process is resonant. In this paper we instead consider situations where two or more distinct resonant processes can occur at the same time, leading naturally to competition between them. Here we show how the theoretical formalism used to obtain time-independent Hamiltonians for single-resonance problems can be generalised to address cases with multiple co-existing resonances. We introduce a compact analytic description of the resulting RWA Hamiltonians and show that it leads naturally to an efficient description of the system's classical dynamics.
We illustrate our analysis by investigating in detail a specific example of competing resonances: a two mode system where a single tunneling Cooper pair can generate either two photons in the first mode or one photon in the second mode (see Fig. 1). We find that the quantum dynamics doesn't produce a clear 'winner' in the competition between resonant processes, instead they can coexist with similar strengths. Furthermore, the coexistence of the resonances generates effective couplings between the modes which can lead to significant entanglement.
The rest of this article is organised as follows. We start by introducing the theoretical model for the JJ-cavity system in Sec. II. In Sec. III we show how special functions can be used to obtain compact expressions for RWA Hamiltonians describing competing resonances and the corresponding classical description is derived in Sec. IV. Then in Sec. V we explore the quantum dynamics that arises for the example with two co-existing resonances. Finally, we conclude in Sec. VI.

II. MODEL SYSTEM
We consider a system of N harmonic modes, with individual frequencies ω 1 , . . . , ω N , in series with a JJ and with a voltage bias V applied, as sketched in Fig.1(a). The modes could be different harmonics within one or more microwave cavities 1-3,6,22 , or they could be realised as lumped element LC-oscillators 4 . The circuit can be described by the following time-dependent Hamiltonian 9 where ω J = 2eV / ,â n are the annihilation operators of the modes, ∆ n the zero-point displacement (determined by the corresponding mode capacitance, C n , and inductance, L n ) ∆ n = (2e 2 L l /C l / ) 1/2 and E J is the Josephson energy of the junction. Almost all of the parameters in this circuit can be varied, either through circuit design 4 (ω n , ∆ n ), or in-situ within a given device e.g. via a change of voltage 1 (ω J ). The value of E J can be tuned in-situ by using a parallel combination of two JJs (SQUID) and applying a flux bias.
The time-dependence makes Eq. (1) a difficult Hamiltonian to work with. In cases where only a single mode is included, resonances where ω J is an integer multiple of the mode frequency can be described by an approximate time-independent Hamiltonian obtained via a rotating wave approximation 4,9,10 . A similar method was applied to study two-mode systems with ω J chosen to match the sum of the mode frequencies, defining a single resonance 3,12,19 . We will now consider how this approach can be generalised to problems involving a wider set of modes and allowing for cases where more than one process can be resonant.
Multiple resonances involving a set of N modes arise naturally when their frequencies, and that of the drive frequency ω J , are all commensurate. For convenience, we shall start by assuming that all of the frequencies can be expressed as integer multiples of the fundamental (lowest) mode frequency ω 1 : i.e. the values of q l = ω l /ω 1 with l = 1, . . . , N and p = ω J /ω 1 are all positive integers. Resonances in the system associated with the inelastic tunneling of a Cooper-pair across the junction 23 are then described by vectors m, with N integer components that satisfy N l=1 q l m l = p, with positive (negative) components m l describing the gain (loss) of |m l | photons in the l-th mode. In cases where more than one such vector can be found the system has competing resonances.
For the simple competing resonance illustrated in Fig.  1b we have N = 2 and ω J = ω 2 = 2ω 1 , hence p = 2 and the set {q} = (q 1 , q 2 ) = (1, 2). We can think of this as a competition between two resonances, as to lowest order in the number of photons created/destroyed, creation of either one photon in mode 2 or two photons in mode 1 are both resonant. However, the behavior described by Eq. (1) is rather more complex, and higher order processes involving an exchange between the modes must also be accounted for. In fact, all vectors of the form m (k) = (2k, 1 − k) satisfy the resonance condition with k = 0, ±1, ±2, . . . . This illustrates the basic problem in dealing with competing resonances: as soon as there are two modes with frequencies that are both commensurate with ω J , direct processes in which just one mode, or the other, is excited by inelastic tunneling are accompanied by a whole host of others in which photons are exchanged between the modes. This is a manifestation of the complex mode-mode coupling that the Hamiltonian (1) gives rise to.
In the following we will consider systems where the resonance condition(s) are met up to some small detunings, δ l , such that ω l = (q l /p)ω J + δ l (with q l , p positive integers and q 1 = 1 as before). We proceed by transforming into a rotating frame via the unitary transform: The RWA is then made, assuming that terms that retain a time-dependence in the rotating frame can be neglected. This is equivalent to assuming that only the terms describing (close to) resonant processes need to be retained. The simplest way of expressing the resulting Hamiltonian is to simply pick out the matrix elements in the number state basis that have no time dependence in the rotating frame 11 . For the multi-mode case we can do this formally via a filter which selects only the relevant timeindependent terms. This results in the following recipe for the RWA Hamiltonian with the filter, E, defined by the relation where |n = |n 1 , n 2 , . . . , n N is an N -mode Fock state. The sum over n runs over all states whilst the other sum is over the vectors m belonging to the set S that satisfy the resonance constraint, l q l m l = p, whilst also having n l + m l ≥ 0 for all l. Hence for the 2-mode competing resonance where ω J = ω 2 = 2ω 1 , the set S is over the vectors m (k) = (2k, 1 − k), leading to the states |n 1 + 2k, n 2 +1−k , with k an integer within the range −n 1 /2 ≤ k ≤ n 2 + 1.
In addition to the coherent drive represented by Eq. (3) a model of the system dynamics must also include the inevitable photon leakage from the modes. This could represent unwanted losses, coupling to collection lines or a mixture of the two. For simplicity we assume a standard zero-temperature Lindblad master equation 24 with γ l the loss rate for mode l.

III. SPECIAL FUNCTION FORM OF HAMILTONIAN
The filtering out of the resonant terms to produce a power series embodied by Eq. (3) is a convenient route for numerical calculations, but it is difficult to connect with simpler approximate descriptions based, e.g. on a coherent state ansatz (see Sec. IV below) in particular. Instead it is convenient to derive compact functional forms for the power series of operator terms left after the RWA has been implemented, an approach which is facilitated by the use of normal ordering.
In the simple case of a single mode system where ω 1 = (ω J /p) + δ, the Taylor series of a Bessel function can be identified in the normally-ordered expansion that follows after the RWA is made. This leads to the compact expression 9,10 where : . : indicates normal order, J p (x) is a Bessel function of the first kind of order p andẼ J = E J e −∆ 2 /2 is the renormalised value of the Josephson energy 4,25 .
For single-resonance circuits containing multiple cavities, the normally ordered operator power series in the RWA Hamiltonians can be written as products of Bessel functions, one for each cavity involved in the process 12,19 . We now generalise this approach to find a compact representation for the RWA Hamiltonian for situations where two or more resonances compete. To do so we go back to consider the full Hamiltonian [Eq. (1)] transformed into the rotating frame using Eq. (2). Making the RWA by discarding the time dependent terms leads to complicated algebra which significantly complicates an analytic derivation, but this difficulty is sidestepped by instead introducing the RWA with an integral over the fundamental period, T = 2pπ/ω J , The time dependent terms average to zero over a period, so this integral form is fully equivalent to directly discarding those terms. Formulating the RWA Hamiltonian in this form provides a straightforward way to express it in terms of special functions which can be defined via integrals as we now show.
To simplify Eq. (7) for the N mode system, we introduce special functions denoted Z, defined via the generating function with the colons indicating normal ordering, as usual. The ..x N ) withx l in our case a mode raising or lowering operator, up to a constant factor. The N superscript indices, {q}, together with the subscript index, p, together fully encode the resonance conditions that will need to be incorporated in the reformulation of Eq. (7). These functions are essentially multi-dimensional Bessel functions 26,27 , but with minor modifications to incorporate complex and operator arguments more readily. As with the single mode case, normal ordering removes all ambiguity from the corresponding power series involving operator arguments.
For a single mode (N = 1) case the Z function for the c-number argument Ae iθ is just an ordinary Bessel function multiplied by a phase factor, Z The two mode (N = 2) version of the Z function is closely related to the 2D generalisation of the Bessel function 27 . Many properties of these functions, such as Taylor series representations, derivatives and relational properties are derived in Appendix A. These relations prove to be surprisingly simple, allowing expressions involving the Z functions to be manipulated quite straightforwardly.
A useful integral representation of the Z functions is obtained by setting y = exp(it) in Eq. (8), then inserting a factor of (1/2π) π −π dt exp(−imt) on both sides of the equality. Noticing that on the left the integral reduces to a Kronecker δ-function 27 , one finds Returning to the RWA Hamiltonian, the expression is simplified by splitting the cosine in Eq. (7) into a sum of exponentials, each of which is rearranged to achieve normal order, and then identifying the integral representations of the Z functions, Eq. (9). The Hamiltonian can therefore be expressed aŝ wherex l = 2i∆ lâl and we have redefinedẼ J = E J exp − N l=1 ∆ 2 l /2 . Although apparently rather abstract, Eq. (10) facilitates analytic manipulations, as we demonstrate in the next section.

IV. COHERENT STATE ANSATZ
A coherent state ansatz can be used to obtain a simpler approximate description of the system's dynamics 9,30,31 . The idea is to assume that each mode is in a coherent state, ρ α = N l=1 |α l α l |, described by a complex amplitude α l . Substituting this into the master equation (5) leads to a set of equations of motion for the amplitudes, the fixed points of which provide a valuable framework for understanding the dynamics of the system 9,18,32,33 . This approach can be thought of as providing an essentially classical description of the dynamics as (quantum) fluctuations in the amplitudes are neglected 34 . When applied to systems with a unique resonance, the resulting fixed point amplitudes have been shown to provide an increasingly accurate way of predicting properties like the average occupation numbers of the modes as the strength of the quantum fluctuations (measured by ∆ l ) are reduced 14 . However, one key limitation is that no information is provided about the way in which the density operator spreads between two or more coexisting stable fixed points 18 .
To apply the coherent state ansatz to the competingresonance case we evaluateα l = Tr(â lρ ) using Eq. (5), exploiting the relation [â, f (â,â † )] = ∂f (â,â † )/∂â † , and then substituting in ρ α . The analytic properties of the Z functions (discussed in Appendix A), and in particular their derivatives [see Eq. (A5)], make this a straightforward calculation. The amplitudes are thus found to obey the following coupled set of equations: The classical fixed point(s) are obtained by solving the algebraic equations for α l obtained by settingα l = 0 for all l. The stability of these points (and hence their role in the system's long time dynamics) is determined by the eigenvalues of the Jacobian matrix, If any of the eigenvalues are positive the classical fixed point is unstable. The fact that the Z-functions can be differentiated and evaluated fairly easily is invaluable in locating the stable fixed points of the classical system. In particular, Eq. (9) provides a convenient way of carrying out the numerical evaluations of Z functions with complex arguments that arise when calculating the fixed points and the Jacobians needed to determine their stability.

V. EXAMPLE: TWO MODE COMPETITION
Having obtained formal expressions for the RWA Hamiltonian in cases where competing resonances exist, we now look in detail at the specific example of the twomode problem with p = 2 and {q} = (1, 2) (see Fig. 1b). Our main aim is to gain insight into how competing resonances can affect the quantum dynamics of the system, but the analysis also serves to illuminate the very general formulations presented in the preceding sections.
On-resonance, the RWA Hamiltonian (10) for our two mode system with competing resonances takes the form We note that one can use the properties of the Zfunctions detailed in Appendix A to re-express this as an infinite sum over products of Bessel functions of different orders, or equivalently as an operator power series with three nested summations, but the resulting expressions are unwieldy. However, some insight into the interaction between the modes can be gained by analysing the parts of the Hamiltonian related to the lowest order processes. Including just the terms up to 4 th order in the creation/annihilation operators: withn l =â † lâ l . We can see that two qualitatively rather different effects are present. The first two lines of (14), together with their corresponding Hermitian conjugates, describe processes where photons are added/removed to just one of the modes, but in both cases the effective rates are modified by the photon populations of both modes. This provides a form of nonlinear coupling similar to that arising, for example, in cavity optomechanics 35 . In contrast, the last line of (14) consists of a more direct form of interaction involving the conversion of quanta between the modes, though it occurs at 4 th order in the operators.
It is interesting to note that this effective interaction is very different to that which arises for a single resonance where the energy from a Cooper pair generates photons in two modes simultaneously 6,12,19,29 . In the latter case, the interaction contains terms which are bilinear, reducing to a degenerate parametric amplifier to lowest order in the operators.

A. Fixed Point Analysis
For our two mode case, the equations of motion for the mode amplitudes that follow from (11) arė , with x j = 2i∆ j α j . To obtain the corresponding fixed points we use standard optimisation methods, evaluating the Z functions through numerical integration (9). It is possible to instead proceed by splitting the Z functions into sums over products of Bessel functions (see Appendix A). However, direct use of Z functions, evaluated by integration, has a number of advantages. Firstly, it readily scales to higher dimensions (more modes) 36 . Secondly, it avoids the subtleties of working out where to truncate the (in principle infinite) summations that arise. Indeed, we found the integration method to be much faster in our calculations. The fixed points are given by pairs of values α 1 , α 2 , the amplitudes of which are shown in Fig. 2 as a function of the drive strength, E J , normalised by a threshold value E (T ) J (defined below). Initially there is only one stable solution with the amplitude of mode 2 (which is resonantly driven ω 2 = ω J ) growing linearly at first whilst the amplitude of mode 1 (ω 1 = ω J /2) remains zero throughout. This fixed point represents the case where mode 2 wins completely in the competition between resonances. Indeed, the behavior of α 2 for this fixed point matches exactly what one gets with a single resonantly excited mode 9 : it grows more slowly with increasing E J and its amplitude eventually becomes locked to a constant value (at E J /E (T ) J 1.69). At larger drive strengths the picture changes significantly with a second stable fixed point emerging. A saddle-node bifurcation occurs atẼ J / γ ≈ 6.87, where we have assumed γ = γ 1 = γ 2 . Since the first mode can now become excited, we use this bifurcation point to define the threshold value for the drive strength, E (T ) J . The bifurcation is collective: the amplitudes of both modes change abruptly. The new stable solution has nonzero amplitudes in both modes, though with that of mode 1 significantly larger than that of mode 2. The threshold occurs at a higher drive strength than that which is required to excite a single mode at the two-photon resonance 9 , and hence one can think of the presence of the resonantly driven mode 2 as tending to suppress the excitation of mode 1.
There are in fact two bifurcations that occur simultaneously at the threshold, although the two are identical up to phases leading to pairs of fixed points with matching amplitudes, leading to only one set of curves in Fig. 2 37 . Interestingly, the amplitude in mode 2 of the new stable points initially drops with increasing drive, until it touches zero for E J /E (T ) J ∼ 1.055, after which it grows again. Seen in the full phase space the complex amplitude of the fixed point moves continuously through the origin. We can think of this second stable fixed point as representing a case where mode 1 wins the competition between resonances, winning completely for

B. Quantum Steady State
We now move on to examine the full quantum dynamics of the mode competition, using numerical solutions of the master equation (5) obtained using the QuTiP package 38 . Figure 2 compares the steady state expectation values n 1 and n 2 with the stable fixed point amplitudes. Although the connection between these quantities is apparent at low E J (for mode 2 in particular), it is no longer clear after the bifurcation which leads to bistability with the emergence of the second stable fixed point.
A much clearer understanding of the quantum behavior can be obtained by looking instead at the joint number state probability distribution of the two mode system, given by P (n 1 , n 2 ) = n 1 , n 2 |ρ ss |n 1 , n 2 with ρ ss the steady-state density operator. Several examples of P (n 1 , n 2 ) for different choices of E J are shown in Fig.  3, overlaid with the locations of the corresponding classical fixed points. Whilst there can be more than one stable classical fixed point for a given parameter set, the quantum dynamics always have a unique steady state solution. We see that at low E J the probability distribution is peaked around the location of the only classical stable fixed point, albeit with a significant spread due to quantum fluctuations. For E J > E (T ) J , the probability distribution becomes bimodal with peaks roughly concentrated around the locations of the two co-existing classical stable fixed points. Interestingly, these two peaks have a rather different character: the one corresponding to high occupation of mode 1 (and low occupation of mode 2) is much more diffuse than the one corresponding to high occupation of mode 2 (and low occupation of mode 1). Nevertheless, the overall message is clear: the mode competition has no overall winner in the quantum regime. Instead, both of the classical solutions are represented within the quantum steady state.

C. Mode Correlations
Finally, we examine the correlations that develop between the two modes that ensue as the quantum system combines the two very different outcomes apparent in the bistability of the fixed points. We will look at amplitude correlations within and between the modes and then quantify the entanglement that is generated.
The bimodal number-state distributions that emerge at larger E J values indicate that the photon populations have become anti-correlated 39 . The detection of a photon from one mode means that it is less likely that one will be found in the other. Such effects can be quantified using second order correlation function 3,16,20 The auto-correlations for each mode (i = j = 1, 2) and the cross-correlations (i, j = 1, 2) are shown in Fig. 4. The auto-correlations are what one would expect for uncoupled modes at low E J / γ, with more complex behavior emerging at larger drive strengths. For mode 1, photons are always produced in pairs (ω J = 2ω 1 ) and assuming rare (uncorrelated) pair creation events implies 20,40 g 11 (0) ∼ 1/(2 n ), which matches the low E J behavior very well. For mode 2, photons are produced one-ata-time (ω J = ω 2 ) and a modest anti-bunching of the photons is expected at low E J / γ, taking into account the non-linearity 20 g 22 (0) remains slightly higher than this estimate (at low E J ) ij (0). Full lines are from numerical calculations, dashed lines are low-EJ/ γ estimates for g (2) 11 (22) (0) discussed in the text. and drifts higher still with increasing E J . This is a result of coupling to the other mode which opens up the possibility of a range of higher order processes that tend to promote bunching, e.g., one in which inelastic Cooper pair tunneling generates two photons in mode 2 whilst simultaneously annihilating two photons from mode 1.
The cross-correlation, g 12 (0), remains less than unity throughout indicating the expected anti-correlations. No clear connection to the behavior of the classical fixed points is apparent here, though there is a minimum in g (2) 12 (0) within the bistable region. Furthermore, the anticorrelation means that the Cauchy-Schwartz inequality g (2) 11 (0)g (2) 22 (0) ≥ g (2) 12 (0) is never violated here. This is in contrast to a single resonance where photons are created in pairs with one in each of two modes 12,16,19 , thereby generating positive correlations.
The effective interactions generated by the competing resonances do not just generate anti-correlations in the two modes, they are also able to entangle them even though they are purely nonlinear, involving only terms that are 3 rd order or higher in the creation/annihilation operators [see Eq. (14)]. To demonstrate this we use the log-negativity as a convenient measure of entanglement 6,29 , defined as where the negativity, N , is the absolute value of the sum of the negative eigenvalues of the partial transpose of the density operator 41,42 . A logarithmic negativity exceeding zero is sufficient (though not a necessary condition) to identify a state as entangled. The behavior of E N (ρ) as a function of the drive is shown in Fig. 5. We find that the logarithmic negativity initially grows smoothly with the drive strength, later going through a maximum (before the threshold is reached) and then a minimum, but remaining non-zero throughout. The values of the logarithmic negativity achieved are not especially small given the higher order nature of the processes that give rise to the correlations [see Eq. 14]. The peak in Fig. 4b is only about a factor of two less than ln 2 which is the upper bound achievable in the two-mode squeezed state produced by a coherent parametric amplifier interaction 43 , which is bilinear in the operators.

VI. CONCLUSIONS
We have explored the quantum dynamics of systems in which inelastic tunneling of Cooper pairs across a voltage biased JJ excites a series of microwave oscillators via two or more competing resonant processes. The competing resonances arise when the mode frequencies and the Josephson frequency (set by the bias voltage) are commensurate. The competition between the resonances can be described by a simplified time-independent Hamiltonian using a RWA, following the approach used for cases with a single resonance. However, the resulting Hamiltonians are rather complicated and unwieldy, even for systems with just two modes. The very strong nonlinearity of the system, together with commensurable mode frequencies, mean that a large number of processes that couple the modes together need to be accounted for. We introduce a compact and efficient technique for analysing such RWA Hamiltonians using normal ordering and a generalised special function. We illustrate the utility of this approach by showing how it can readily be applied to obtain simplified (classical) equations of motion for the amplitudes of the modes.
We also explored in detail a simple example in which two resonances compete in a two-mode system. Two stable classical fixed points of the system emerge, each one associated with a different one of the two competing resonances clearly 'winning'. The quantum dynamics reveals a more complex situation in which bistability emerges naturally with contributions from both resonances evident in the steady-state density matrix. Furthermore, although the effective interactions between the modes in the presence of competing resonances are purely nonlinear, they are sufficient to generate a significant amount of entanglement.
It would be interesting to investigate how competing resonances evolve in cases involving more than two modes in the future. Unfortunately, straightforward numerical solutions of the quantum dynamics become less and less tractable as the state space grows with the number of modes. However, the compact formulations of the multimode RWA Hamiltonians developed here should prove a useful starting point for developing analytic approximations.