Multi-Photon Resonances in Josephson Junction-Cavity Circuits

We explore the dissipative dynamics of nonlinearly driven oscillator systems tuned to resonances where multiple excitations are generated. Such systems are readily realised in circuit QED systems combining Josephson junctions with a microwave cavity and a drive achieved either through flux or voltage bias. For resonances involving 3 or more photons the system undergoes a sequence of two closely spaced dynamical transitions (the first one discontinuous and the second continuous) as the driving is increased leading to steady states that form complex periodic structures in phase space. In the vicinity of the transitions the system displays interesting bistable behaviour: we find that coherent effects can lead to surprising oscillations in the weight of the different dynamical states in the steady state of the system with increasing drive. We show that the dynamics is well-described by a simple effective rate model with transitions between states localised at different points in the phase space crystal. The oscillations in the weights of the dynamical states is reflected in corresponding oscillations in a time-scale that describes transitions between the states.


Introduction
Multi-photon resonances at which several photons can be created simultaneously via the down conversion of a high frequency pump are intrinsically highly nonlinear. Such processes are potentially very useful in a variety of contexts such as quantum error correction [1,2]. However, multi-photon processes are also of intrinsic interest. They give rise to quantum states with a variety of novel properties including higher-order squeezing [3,4] and rich periodic structures in phase space, known as phase-space crystals [5,6].
In this paper we explore the dissipative dynamics of an oscillator that is driven nonlinearly so that a handful of photons (up to six) are generated at a time via a resonant process. Our work builds on several recent studies of 3-photon resonances [7,8,9,10,11] as well as investigations of the properties of the eigenstates of the Hamiltonian that display periodic structures in the phase space, known as phase-space crystals, generated at very high-order resonances (typically more than ten photons) [5,12], though the latter focused mainly on the properties of closed systems. We explore how the properties of the oscillator evolve as a function of the applied drive strength, the number of photons that are created together and the strength of the zero-point fluctuations.
For resonances involving three or more photons the system undergoes a pair of closely spaced dissipative phase transitions: a first order one associated with the breaking of rotational symmetry in phase space and then a continuous one linked with a chiral symmetry of the system Hamiltonian. We find that novel features emerge for resonances involving more than three photons, with the system exhibiting a form of bistability in the vicinity of the dynamical transitions. In this regime the steady state consists of a mixture of states in the phase space with different amplitudes and phases. The weights of the different states oscillate as a function of the drive, leading to oscillations in the average occupation number of the oscillator. We examine the eigenspectrum of the Liouvillian of the system and construct a simple effective model for the slow dynamics of the system which reveals corresponding oscillations in the rate at which the system moves between coexisting high and low amplitude states. We also analyse the eigenoperators of the Liouvillian which have a Bloch-like character in phase space, demonstrating how the concept of a phase space crystal can be applied in open systems.
Although multi-photon resonances can be found in a variety of different systems, the exceptionally strong non-linearities achieved in superconducting circuit devices make them especially suited to exploring this kind of physics [1,4]. The specific nonlinearly driven oscillator model we analyse here can be realised by voltage-biasing a Josephson junction (JJ) in series with a microwave cavity. Photons are generated via inelastic tunnelling of Cooper pairs across the JJ and the system can be tuned to resonances where the creation of different numbers of photons is favoured by simply adjusting the bias voltage [13,14,15,16,17,18,19]. Furthermore, the resonances involving between 3 and 6 photons that we investigate are expected to be readily accessible with current devices. However, very similar effects can be achieved using slightly different architectures. For example, closely related Hamiltonians have been engineered by using a time-dependent flux bias instead [20,21].
The rest of this article is organised as follows. Section 2 introduces our nonlinearly driven oscillator model in more detail. Then in section 3 we analyse the classical fixed points of the system which provide a framework for analysing the quantum dynamics. We then examine the properties of the full (quantum) steady-state in section 4, looking in particular at how it depends on the drive strength and the number of photons generated. Next, in section 5 we explore the eigenspectrum and corresponding eigenoperators of the Liouvillian of the system. This allows us to go on to develop a simple description of the dynamics in section 6, as a small set of dynamical states linked by transition rates between them. We summarise and conclude in section 7 and further details for certain aspects of the work are provided in Appendices. consider. An LC oscillator with frequency ω = 1/LC is driven at a frequency ω J = 2eV /h by a voltage-biased Josephson junction, characterised by the Josephson energy E J . The oscillator is assumed to be damped at a rate γ. (b) Sketch of the time-dependent potential in the Hamiltonian. It consists of a harmonic quadratic part and a moving sinusoidal potential [22].

Non-Linearly Driven Oscillator Model
The model circuit we consider is shown in figure 1, it consists of a series combination of a JJ and cavity described by a single LC oscillator with frequency ω 0 = 1/ √ LC (other cavity modes are assumed to be far detuned from resonance), though it could equally be realised with a JJ coupled to a lumped element LC oscillator [19]. The Hamiltonian of the system is given by [14,15] whereâ is the lowering operator of the mode, E J is the Josephson energy of the junction, ω J = 2eV /h is the Josephson frequency set by the bias voltage and ∆ = (2e 2 L/C/h) 1/2 . This Hamiltonian can also be regarded as describing a particle in a quadratic potential perturbed by a travelling wave [22], and as a result the same Hamiltonian emerges in the study of cold atoms in time-varying optical traps [12]. We consider the case where the voltage is tuned close to resonances where the inelastic tunnelling of a Cooper-pair can generate p-photons in the mode, so that ω J pω 0 . Moving to a frame that rotates at a frequency ω J /p, and making a rotating wave approximation (RWA) [14] leads to the Hamiltonian where δ = ω 0 − ω J /p is the detuning from resonance, the colons, ::, indicate normal ordering and J p (x) is a Bessel function: J p (x) = Σ n (−1) n (x/2) p+2n /n!(n+p)!. Although the RWA approximation has been found to be a very good approximation for a wide range of parameters in this system [23,19], we have explicitly checked that it also works well for the parameter regimes we explore here, see Appendix A.
The RWA Hamiltonian has two important symmetries [12]. First, it commutes with the operatorr = exp(iâ †â (2π/p)), which follows from the batch creation of p-photons. In phase space this commutation manifests as discrete rotational symmetry [12,8]. This is the main ingredient of phase space crystals: eigenstates of the Hamiltonian inherit this symmetry, which permits a description in terms of Bloch modes [5,12,7]. The second symmetry is that the operatorĉ =r 1/2 anti-commutes with the Hamiltonain, cĤ RWA = −Ĥ RWAĉ . This is chiral symmetry, half-period rotation leading to a sign change [12].
We include the effects of photon losses from the cavity together with the coherent evolution of the system via a Lindblad master equation, assuming a zero-temperature environment for simplicitẏ with γ the loss rate. In fact, the fixed bias voltage we have assumed is an idealisation: in practice the voltage across the JJ and the cavity is not completely fixed, but subject to fluctuations due, e.g. to the presence of additional impedances in the circuit [15,24].
Here we concentrate on the minimal description of the system dynamics (3), but our approach can be extended to account for the effects of voltage fluctuations and we consider how they affect our results in Appendix B.
The discrete rotational symmetry of the RWA Hamiltonian (though not the chiral symmetry) is preserved in the master equation [25]. In this case, [L, R] = 0 where R · =r ·r −1 . This symmetry does not give rise to a conserved charge as a conserved charge requires the unmet, stronger condition that the collapse operator,â, commutes withr [26].
Within the master-equation description, only two timescales remain when the voltage is tuned to resonance so that δ = 0. Their ratio, E J /hγ, sets the strength of the driving compared to that of the losses. The value of E J can be tuned in-situ by using a slightly more complex SQUID set-up [14] which acts as an effective single junction with a Josephson energy that can be tuned via an applied flux. The quantity ∆ determines the strength of the quantum fluctuations, it measures the magnitude of the zero-point fluctuations in the oscillator flux in units of the flux quantum. Whilst ∆ is fixed within a given experiment, recent developments in device engineering mean that it can be varied over a relatively wide range [19]: from ∆ 1 up to ∆ ∼ 1. Looking at the Hamiltonian (2), we see that ∆ mediates the strength of the Bessel function term. This describes how the presence of photons in the oscillator affect the creation of further photons: for ∆ ∼ 1 this becomes significant even at the level of a single photon. The other parameter that appears is p, the number of photons generated at a time at the resonance with ω J = pω 0 .

Classical Fixed Points
Before looking at the quantum dynamics of the system we examine the corresponding classical dynamics. Identifying the classical fixed points of the system, together with the sequence of bifurcations that arise as the driving strength is increased progressively, provides a very useful framework for analysing the open quantum dynamics.
Equations of motion for a complex classical amplitude α can be obtained from the master equation (3) by simply making the ansatz that the system is in a coherent state [27] (ρ → |α α|). A more sophisticated semi-classical approximation can be developed by transforming the master equation into an equation of motion for the Wigner distribution (WD) of the density operator, W (α, α * ), and dropping derivatives beyond second-order which leads to a Fokker-Planck equation (see Appendix C for details). Although more involved, this route provides a more systematic approach to understanding the classical limit of the quantum dynamics [9].
The ways in which the fixed points evolve with drive strength and the bifurcations that emerge for different values of p are necessarily constrained by the corresponding symmetries. For p = 1 the fixed point which is at zero-amplitude (A = |α| = 0) for zero drive moves away from the origin smoothly with increasing drive (higher E J ). However, for p = 2 the π-rotational symmetry means the fixed point cannot move from the origin. Instead, at a particular value of E J , the origin loses stability [14], with the simultaneous emergence of two stable fixed points in opposite directions-a pitchfork bifurcation [28].
For p > 2 there is always a stable point located at zero amplitude which we call the Dark Point (DP). It cannot be displaced from the origin or undergo bifurcations while maintaining 3-fold or higher rotational symmetry ‡. However the DP is not always the only stable point in these systems.
For p ≥ 3, a set of p new stable points are generated (together with p corresponding saddle points), with a common amplitude A r and phases that differ by 2π/p, via a set of saddle node bifurcations. The bifurcations occur at a threshold drive E (T ) J = ‡ This argument is quite general and applies to other damped systems driven at an overtone above 2 [20,7,11] (hγA 2 r e ∆ 2 /2 )/[pJ p (2∆A r )] with the amplitude given by the smallest nonzero solution to ∆A r J p (2∆A r ) = J p (2∆A r ) (see Appendix D for details). We call these stable points Bright Points (BP), and it is a key characteristic of the multiphoton resonances with p ≥ 3 that these points always coexist with the stable DP. The values of both E (T ) J and A r increase with p, but whilst A r is proportional to 1/∆, E (T ) J instead scales as exp(∆ 2 /2)/∆ 2 .
As E J is increased beyond E (T ) J , the BPs migrate to higher amplitudes (with the corresponding saddle points moving to lower amplitudes) and eventually a second set of bifurcations take place. In this case pitchfork bifurcations occur as each BP splits into three, two stable points that remain locked at a fixed amplitude and thereafter evolve only in phase with increasing drive, and an unstable (saddle) point that continues to increase in amplitude. In contrast to the saddle node bifurcations, these bifurcations occur for all p values [14,29].

Steady State Properties
The strong nonlinearity of the multi-photon resonances mean that numerical methods § must be used to explore the quantum properties of the system. The first question that we address is how the steady-state properties of the system, in particular the occupation number n = â †â , behave and how this compares to the underlying structure of classical fixed points and bifurcations. Then in the next section we investigate the key dynamical time-scales in the problem by calculating the spectral properties of the Liouvillian. Figure 2 shows the photon-number distributions of the steady states as a function of the drive strength, E J , for p = 4, 5, 6 and ∆ = 0.4 and ∆ = 1.0. The plots also show the behaviour of n and the classical fixed points.
The evolution of the classical fixed points with E J follows a very similar pattern for different values of p and ∆. In each case, the first bifurcation occurs at E (T ) J , leading to the new stable BPs appearing well away from zero, followed by the second bifurcation after which the BP amplitude remains fixed. This fixed amplitude is given by [14] A p = z p /(2∆) with z p given by the first zero of J p (z p ), so that z p = 5.318, 6.416, 7.501 for p = 4, 5, 6. Since for all but the largest ∆ and p values, the average occupation numbers match up well with the prediction of the fixed point amplitudes in the limit of large E J , this provides a simple measure of how the occupation numbers in the problem increase with p, but decrease with increasing ∆. Another important effect of increasing p is that it brings the two bifurcations closer together (indeed they merge in the infinite-p limit).
Although the behaviour of the average occupation number comes close to the corresponding classical fixed point for large enough E J in almost all the plots shown in figure 2, the differences are much more marked at lower drive strengths. In particular, the point at which the system displays a threshold, marked by a rapid rise in n differs § We used the QuTiP package [30] for this task. progressively from the classical prediction E (T ) J as p is increased for ∆ = 0.4. For ∆ = 1 the behaviour is qualitatively different with n developing a series of peaks as E J is increased which become more pronounced at higher p . These oscillations are remarkable in that they imply that increasing the drive (or equivalently reducing the damping rate) can lead in places to a reduction of the occupation number.
The photon number distributions shown in figure 2 reveal significant bimodality over a broad range of drive strengths. In each case, a peak is always present for zero We did not find any clear signatures of these oscillations for p < 4.
photons with a second peak emerging around a value associated with the corresponding classical fixed point (i.e. the BP). For ∆ = 0.4 the peaks are well separated, becoming more distinct with increasing p and the weight of the distribution shifts progressively away from the zero-photon peak as the drive is increased. In contrast, for ∆ = 1.0 the peaks are much less well separated and we see that the oscillations in n involve oscillations in the weight of the distribution back and forth between the two peaks.  We trace the evolution of the oscillations in n as a function of E J and ∆ for p = 6 in figure 3. This shows that the oscillations form a series of curving resonances in the E J -∆ plane, developing for ∆ ≥ 0.45, then becoming less distinct for ∆ > 1, suggesting that the peaks in the underlying distribution need to be neither too far apart, nor too close together, for the oscillations to be clearly pronounced. For all values of ∆ the first of the peaks occurs at the same E J /E (T ) J , between the two classical bifurcations, and the peaks eventually disappear for large enough drive strengths.
Further insights can be obtained by looking at the evolution of the system in phase space using the WD of the steady states, as shown in figure 4. The WDs are shown together with the classical fixed points and the flow streamlines (classical trajectories) for a sequence of increasing E J values ¶. In each case the steady state consists of 'blobs' localised about the classical attractors [29,32,33]. The drop in photon population as E J rises is clearly visible in the almost total drop in the magnitude of the WD around the BPs between panels (b) and (c).
The oscillations in n are a feature of the full coherent quantum dynamics of the system. We could find no trace of them using a semi-classical description in ¶ Note that the WDs are almost entirely positive. Panels (a), (b), (c) of figure 4 have negativity volumes [31], of 1 × 10 −2 , 1.2 × 10 −4 and 8.7 × 10 −2 respectively, too faint to display in the figure itself. terms of the Fokker-Planck equation for the Wigner function (obtained by dropping higher order derivative terms in the equation of motion for the WD) which incorporates diffusion on top of the classical dynamics (see Appendix C for details). Furthermore, the coherences in the number state basis play a crucial role. Introducing numberdephasing terms into the master equation to mimic the effects of voltage noise that arises in experiments [34,15,19] with voltage bias JJ-cavity systems, we find that the sharp resonances in n are progressively washed out as the strength of the dephasing is increased (see Appendix B).

Liouvillian Eigenspectrum
The change of the steady state from one localised at the DP to one localised about the BPs is surprisingly complex, associated as it is with an extended region of bimodality in which the weights of the distribution around the different points can oscillate as a function of the drive. To understand more about this behaviour we turn now to the properties of the Liouvillian super-operator, defined in (3).
The Liouvillian (L) transforms one operator into another [25]. It has a set of (notnecessarily Hermitian) eigenoperators, ρ n which satisfy: with λ n the corresponding eigenvalue. As the evolution of the density operator is determined by the Liouvillian, e Lt , these eigenoperators all have a simple time evolution: ρ n (t) = e λnt ρ n (0). The steady state is the eigenoperator with eigenvalue zero. There is only ever one such state for systems like ours with annihilation-operator dissipation [35]. All other eigenvalues have negative real part, so that an arbitrary initial state expressed as a linear combination of eigenoperators converges towards the steady state under time-evolution [36,25]. Time evolution conserves trace and quasiprobability, thus the steady state can be normalised so that its density operator has trace 1 and its WD has integral 1. In contrast the transients all require a trace/integral of 0. Dissipative phase transitions [37,36,25] are associated with the emergence of one or more very slow time-scales which arise when parameters of the system are tuned, leading to a sharp narrowing of the Liouvillian gap, given by |Re[λ 1 ]| where the eigenvalue λ 1 is the (non-zero) eigenvalue whose real part is least negative. Formally, a phase transition occurs in a thermodynamic limit, associated with a divergence in occupation number for nonlinear oscillator systems [25], where the Liouvillian gap vanishes.
The eigenoperators, ρ n , inherit the symmetries of the Liouvillian. In our case, there is a discrete rotational symmetry in phase space, described by the rotation superoperator, R. Since [R, L] = 0, eigenoperators of L will also be eigenoperators of R with Rρ n = e ikn ρ n , with k n an integer multiple of 2π/p. We will see that this symmetry plays a crucial role in applying the concepts of phase space crystals in an open system + . Seen in the lab frame this symmetry corresponds to time-translation symmetry [38].
In the following we start by analysing the behaviour of the eigenvalues of the Liouvillian and the connections to dissipative phase transitions in detail. We also uncover clear connections between the behaviour of the eigenvalues and the oscillations in the occupation number discussed in the previous section. We then go on to explore the properties of the corresponding eigenoperators.

Eigenvalues
The behaviour of the Liouvillian eigenvalues with least negative (but non-zero) real parts is shown in figure 5 as a function of the drive strength for p = 3, 4, 5, 6. These plots illustrate in particular how the emergence and subsequent evolution of slow timescales in the system becomes connected with the oscillations in the average occupation number that develop for increasing p values.
We start by considering the case where p = 3 figure 5(a), where oscillations in n are not seen, and instead the log-scale makes clear that n starts to saturate for J . There is also a clear descent in the real part of several eigenvalues towards zero as E (T ) J is approached, with one eigenvalue (with rotational eigenvalue k n = 0) then rising up again [10] whilst two others (with k n = 0) subsequently stay very small (i.e. with a real part that is very small in magnitude), albeit with a weak but noticeable oscillating component. This behaviour broadly matches what one would expect to see in the vicinity of a first-order dissipative transition in which symmetry breaking also occurs [25], and fits our expectations based on the underlying classical bifurcation. The features get sharper if one decreases ∆, which increases n overall, taking us closer to the expected thermodynamic limit. Furthermore, a second set of eigenvalues start to drop towards zero signalling another dissipative transition around the onset of the second bifurcation. + The rotational invariance of the collapse operators means that the Liouvillian does not inherit the chiral symmetry present in the Hamiltonian. The behaviour becomes more complex as p is increased, see figure 5(b,c,d) and fits less well with the standard behaviour associated with dissipative transition paradigms * . Although we still see a group of eigenvalues dropping towards zero around E J ∼ E (T ) J , as p increases the clear separation of the eigenvalues into those that stay very small and one that rapidly grows again (associated with a first order dissipative phase transition) breaks down. Instead, for p = 6 the set of 6 eigenvalues that drop towards zero stick together and then develop marked oscillations over a broad range of drive strengths. These oscillations match those seen in n , revealing an apparent connection between the steady state behaviour and that of some of the slow dynamical time-scales in the problem. We will explore this connection in more detail in section 6, but for now we * The most natural explanation is that we are in fact moving further away from the thermodynamic limit even though n actually increases with p. Interestingly, this in turn implies that the way in which this limit is approached varies with p.   turn to look in more detail at the properties of the groups of eigenvalues that cluster together in figure 5(d) and properties of the corresponding eigenoperators.

Eigenoperators and Band Structure
The eigenoperators of the Liouvillian reveal structures in phase space which closely resemble the structure of Bloch modes in a crystal. These Bloch-modes are similar to the states identified as phase space crystals [5,12] for nonlinearly driven oscillators that are not subject to dissipation. However there are important differences as well. The modes discussed in the phase space crystal literature are simultaneous eigenstates of the Hamiltonian and state-rotation operator,Ĥ RWA |ψ = E |ψ ,r |ψ = exp(ik n ) |ψ [5]. In contrast, for the dissipative case, we are interested in operators that are simultaneous eigenoperators of the Liouvillian and operator-rotation, Lρ = λρ, Rρ = e ikn ρ [10]. Furthermore, whilst Bloch-like Hamiltonian eigenstates are distributed across the Hamiltonian maxima/minima [5], the Liouvillian eigenoperators we investigate here are distributed across the classical fixed points of the system. In figure 6 we use the Wigner representation to depict the eigenoperators associated with the 7 eigenvalues with smallest (in magnitude) real parts. Note that transforming a Hermitian operator into phase space with the Wigner map produces a real valued field, so that a physical density operator corresponds to a real-valued Wigner function. However non-Hermitian operators give complex WDs in phase space, sometimes called non-diagonal Wigner functions [39,40]. In our case all operators with k n neither 0 nor π have complex WDs, as they must to have the correct (complex) eigenvalue under phase-space rotation. In each case the eigenoperators consist of parts centered at the classical fixed points, but with phases chosen to obey the rotational symmetry.
Taken together, the plots in figure 6 show that the group of points oscillating together in eigenvalue in figure 5 correspond to a set of Bloch modes with differing k n -value, but otherwise much alike. Their close lying eigenvalues suggest that coupling between BPs is very weak so that k n has very little effect on lifetime. Only the k n = 0 mode [mode (1) in figure 6, shown as black crosses in figure 5] deviates significantly from the others: this difference arises because this mode is uniquely able to hybridise with the central DP ( figure 6). This splitting increases progressively with E J (though at a rate that depends strongly on p), resulting eventually in a clear splitting-off of this eigenvalue from the others at higher E J as can be seen in figure 5.
There is an equivalence between the clockwise and anti-clockwise directions in phase space for this Liouvillian, manifest as a degeneracy in eigenoperators of L with equal and opposite rotational eigenvalues, k n . By taking linear combinations inside these degenerate subspaces one could choose an alternative basis, for example one where all the eigenoperators of L were Hermitian giving real WDs. However these eigenoperators would lack the rotational symmetry inherent to the problem, obscuring the important role this symmetry plays (the implications of basis choices are discussed in Appendix E).
In figure 7 we move to a higher value of E J , past the second bifurcation (E J = 2.12E  figure 6 to larger drive values) and a new group which entered view shortly after the second bifurcation and which get progressively smaller (i.e. corresponding to longer lived excitations) as E J is increased further. Plotting the WDs of the corresponding eigenoperators we see that the first set is analogous to the acoustic band of a diatomic crystal, while the other (shorter lived) set parallels the optical band [41]. The key difference being in the former case the two 'atoms' in each unit cell are in phase with one another while in the latter they are in anti-phase. The emergence of the second group of nearly degenerate long-lived transients seen in figure 5 is clearly associated with the formation of these optical modes.

Effective Description
The Liouvillian possesses an (in principle) infinite eigenvalue spectrum. However, as seen in figure 5 a handful of these eigenvalues lead to decay rates of the associated eigenoperators that are much slower than γ/2 and well separated from the rest of the spectrum. This separation of timescales implies that the system relaxes rapidly towards a slow-dynamics subspace of much reduced dimension [36,42]. In this section we show that the eigenoperators within this slow-subspace can be used to construct an effective description of the system's dynamics within phase space which in turn reveals that the oscillations in the occupation number are associated with oscillations in the rate   describing motion of the system between the BPs and the DP. The long lived eigenoperators shown in figures 6 and 7 are concentrated in the vicinity of the stable classical fixed points. This implies the evolution of any initial WD will follow two stages, a rapid decay towards the nearest fixed point(s), and then a slow relaxation in which the weight of the WD is adjusted between the different fixed points to generate the correct balance of quasiprobability between these fixed points corresponding to the steady state.
This suggests that a simple effective model of the system's long term dynamics can be constructed. In this model the dynamics is encapsulated in a set of rate equations describing the rates at which the WD quasiprobability flows between the different states associated with each of the stable fixed points, assuming that each fixed point has rate constants connecting it to its nearest neighbours as illustrated in figure 8(a). This implies four rates: one describing motion from the BPs to the DP, Γ in , one for motion from the DP to the BPs, Γ out , one for motion from one BP pair to a neighbouring pair, Γ l , and finally (above the second bifurcation) a rate for motion from one member of a BP pair to its partner, Γ c . This amounts to a reformulation of the problem from the reduced eigenbasis associated with the slow rates, to a 'physical basis' associated with states localised in phase space.
The rates can be determined by mapping between the reduced eigenbasis and the physical basis. For example, the k n = 0 eigenoperator [see, e.g., figure 6], consists of a central point with a positive amplitude surrounded by satellites of the opposite sign. This mode, like all others except the steady state, has a negative eigenvalue, so that time evolution results in uniform decay. Under time evolution WD quasiprobability is conserved locally [43,44,45], so the decay of the eigenmode must proceed by the positive and negative populations finding one another and annihilating. Using this to write down the rate equation in the physical basis, we can establish that the rate for transitions from a BP to the DP, Γ in and the reverse rate, Γ out , must be related to the corresponding (real) eigenvalue as follows: λ (1) = −Γ in −pΓ out . As shown in Appendix F similar arguments allow all of the decay rates in the simplified model to be determined in terms of the corresponding eigenmode decay rates.
The resulting rates are shown in figure 8 for the case where p = 6. We see that the rate of escape from the DP grows monotonically, and appears to be primarily responsible for the overall upward trend in n with E J . In contrast the rate inwards towards the DP oscillates and is clearly responsible for the peaks/oscillations in photon occupation. We can conclude that the peaks in photon number are driven by the suppression of a process, not the enhancement of one -a reduction in snakes, not an increase in ladders.
The rate describing coupling from one BP-pair to a neighbouring one (Γ l ) is significantly smaller than all the other rates (so much so that we have difficulty resolving it from zero) † †. This is exactly what we anticipated, based on the narrow band-width for the k n = 0 states seen in figure 6. Interestingly, the round-about process where amplitude leaves one BP pair to enter the DP (Γ in ), then leaves the DP to enter another BP pair (Γ out ) greatly dominates over direct movement from one pair to a neighbouring one.
In addition to their impact on the equilibrium photon number these changing rates would have other detectable signatures. In equilibrium the system (for many parameter Except Γ 1 which is extremely small throughout. † † Though, interestingly it can be significantly enhanced when dephasing to model fluctuations in the bias voltage is included, see Appendix B choices) occupies a mixture of the DP and the BPs. If one monitored the photon emission in these regimes one would see periods of high activity, where many photons are detected, corresponding to BP occupation. These bright periods would last an average duration of 1/Γ in before giving way to dark periods of reduced activity which, in turn, would on average last 1/Γ out before the system flipped back to the bright phase. This flickering would imprint itself on, for example, a g (2) measurement [46,19].

Summary and Future Perspectives
We theoretically explored a superconducting system, comprised of a cavity coupled to a Josephson junction. We studied regimes where the dynamics lead to the creation/annihilation of up to 6 photons at a time. The high level of rotational symmetry in phase space had several important effects. In the semi-classics it resulted in an 'invincible' stable point at zero-amplitude, while in the quantum description it resulted in the Liouvillian eigenoperators having a Bloch mode structure in phase space, a structure one could describe as a dissipative phase space crystal.
An unusual quantum effect, where the expected number of photons in the cavity depended nonmonotonically on the drive/loss ratio was discovered and investigated. This behaviour arises from the full quantum coherent dynamics of the system and is not seen in the semi-classical limit described by a Fokker-Planck equation. We found that the oscillations in the occupation number were reflected in complex oscillatory changes in the rate at which the system can return to the vacuum from an excited state.
We hope that our work will stimulate detailed experimental studies of multi-photon resonances in superconducting circuits systems. The main features discussed throughout the paper should be experimentally accessible in a variety of currently available device architectures. For coplanar waveguide cavities like that in [16] the modest ∆ means one needs a high, but attainable, drive-loss ratio (E (T ) J /(hγ) ∼ 3500). High impedance resonators with ∆ ≈ 1 [19], require much lower drive-loss ratios, E (T ) J /(hγ) ≈ 9.7 for p = 6. Resonances below 6 will be more easily accessed. The key idealisations made by our model are assessed in appendices Appendix A (RWA) and Appendix B (constant bias voltage), in both cases the results are not expected to change significantly when relaxing these assumptions. It would be interesting to investigate the extent to which oscillations in the cavity occupation number arise in set-ups that exploit flux rather than voltage bias [20,47,21,48] where clear evidence of several higher-order resonances (and corresponding periodic structures in the phase space) has already been seen [20,47,21].
Finally, we note that from a theoretical perspective, future work is needed to provide a more intuitive understanding of how the oscillations in the photon occupation number arise. The fact that they only occur at relatively large values of the zero-point fluctuations where semi-classical approaches seem to break down makes the problem more challenging, but also more interesting.

Appendix A. Rotating Wave Approximation
Throughout this paper we have made use of the rotating wave approximation (RWA), including only stationary terms in the Hamiltonian within the rotating frame. In this Appendix we check the validity of this approach for multiphoton resonances by comparing with calculations carried out using the full time-dependent Hamiltonian.
For the non-RWA calculations, we assume the initial state is the vacuum, and this is evolved forwards in blocks of 720 drive periods until it shows no significant change over a block. The expected photon number and photon number histograms are averaged over a single period. RWA and non-RWA results for p = 6 are compared in figure A1. It can be seen that even for the relatively high quantum fluctuation strength chosen (∆ = 1) there is quite close agreement. The main difference is that the sharp peak in excitation number, while still clearly present, is somewhat less pronounced in the non-RWA solution.
We note, however, that the RWA eventually fails at very low values of E J (invisible on the scale of figure A1). This is because very far below threshold the resonant p-photon process is suppressed practically to zero, so that the (very far) off-resonant 1-photon process neglected in the RWA is no longer completely negligible in comparison [27].

Appendix B. Voltage Noise
In the main text we assume that the bias voltage, V , is fixed. This is an idealisation of the situation found in experiments where additional impedances in the circuit give rise to fluctuations in the voltage seen by the JJ-cavity system [15,34,24]. In this Appendix we use a simple approximate description [34] which maps the voltage fluctuations onto a fluctuating cavity frequency which in turn leads to photon-number dephasing.
Assuming that the dephasing due to the voltage fluctuations is very weak [34] leads to the modified Liouvillian: withn p =â †â /p. Results obtained using this Liouvillian are shown in figure B1 for the case where γ VN /γ = 0.05, chosen to be similar to that in some experiments [15]. Notice that the sharp features in n have been suppressed by the voltage noise, though the peak structure is still clearly present. Figure B1 also shows the corresponding eigenspectrum of the modified Liouvillian.  Figure B1b shows the impact of the voltage noise on the rates for our simple effective model (see sec. 6). The coupling from one of the BP pair to its neighbours, Γ l is much higher. While it remains the slowest rate it is now large enough to fit a value to it. Recall that, without voltage noise, this rate is so small that we could not reliably resolve it from zero. The increase in this rate is perhaps not surprising as the main effect of photon-number dephasing is essentially angular diffusion.

Appendix C. Wigner transform
This Appendix describes the quantum Fokker-Planck equation used for a semi-classical understanding of the system.
The quantum Fokker-Planck equation describing the evolution in phase space is obtained from the master equation (equation 3) with the following substitutions [39,40]: hence an appropriately symmetrised operator (indicated by the subscript S) becomes the corresponding function of the complex amplitudes. The star product is defined as: with the arrows indicating the direction in which the derivatives act and the exponentiation of the derivatives is understood in terms of the Taylor series. It is worth explaining the nature of these replacements. The density operator (normally represented as a matrix) is replaced with the WD, which contains identical information in an alternative format [49]. The new expression of the Hamiltonian, H, is given by the Wigner transform of the Hamiltonian operator. This is gained by first arranging the Hamiltonian such that its operators are in symmetric order (indicated by the subscript S in equation C.1). Then the appropriate Hamiltonian for the Fokker-Planck model is found by a simple substitution ofâ andâ † with α and α * [49].
For our Hamiltonian moving between normal order (equation 2) and symmetrical order changes the expression very little. The only change is that the factor of exp(−∆ 2 /2) present in the normally ordered version disappears. The easiest way of confirming this is to show that the full, time-dependent, Hamiltonian in equation (1) (which is already expressed in symmetrical order) gains this factor when re-arranged into normal order using the Baker-Hausdorff rule [50]. Although this relationship can also be confirmed in the rotating frame using the expressions in [49].
These substitutions produce the following equation of motion for the WD with H =hδA 2 − E J J p (2∆A) cos(pθ + pπ/2). This Hamiltonian is identical to the classical rotating wave Hamiltonians given in [12,29] to describe equivalent systems.
The terms related to first derivatives in (C.3) (both from the Hamiltonian and loss) are drift terms. The coefficients of these terms give the classical trajectories.
These trajectories are almost identical to the trajectories predicted by the coherent state ansatz, The difference is that the coherent state model posses an extra factor e −∆ 2 /2 appearing next to the terms that derive from the Hamiltonian. The reason for this difference is that in the Fokker-Planck model the drift evaluated at some specific point (α, α * ) reflects the drift motion for the part of the WD at precisely that location. In contrast the coherent state ansatz equation of motion evaluated at (α, α * ) reflects the motion of a coherent state (Gaussian WD) with centre at that location.
The drift terms pull the WD towards the attractors (stable fixed points). However, they compete with the second derivative, diffusion, terms which encourage it to spread out. The third-and-higher derivatives are quantum terms, when they are neglected (an approach known as the truncated Wigner approximation) we obtain a Fokker-Planck equation as discussed in the main text.
In the phase-space description the classical approximation amounts to taking only the first term in the Taylor series of the sin function in (C.3) [39]. In our case the parameter ∆ sets the radial scale of H relative to fundamental scale,h. Thus the parameter ∆ can be thought of as the system's 'quantumness' [12]. Similarly p sets the angular scale, so quantum effects are stronger at higher p (each differentiation of H with respect to angle will bring out another factor of p).
The Fokker-Planck (FP) equation for the WD, including only 1st and 2nd derivatives, was solved numerically using the FIPY package [51]. Expected photon numbers are plotted in figure C1. Lacking higher derivatives these Fokker-Planck solutions are essentially classical, but with added noise. For the lower choices of p these solutions match up almost perfectly with numerical solutions of the full master equation as can be seen in panels (a, b, e, f) of figure C1.
For p = 3 we see some deviation, with the FP solutions remaining near the DP slightly longer. This trend grows more pronounced at p = 4. Importantly the Fokker-Planck solutions lack the peak structures, confirming the quantum origin of these features.

Appendix D. Stability and Bifurcations
In this Appendix we obtain analytic expressions for the threshold E J values at the saddlenode bifurcations which mark the first appearance of fixed points away from the origin (for p > 2). These bifurcations occur when both the drift terms in the Fokker-Planck equation and the determinant of the corresponding Jacobian are zero.
The drift terms of equation C.3 are conveniently expressed in cylindrical coordinates using D is defined so thatẆ = ∇ · ( DW ) + O(∂ 2+ W ), recalling that in cylindrical coordinates ∇ = (1/A)∂ A A (1/A)∂ θ . The first (second) element is the radial (angular) part. The Jacobian related to this drift is: Choosing zero detuning (δ = 0) implies that for low amplitudes the first fixed points to emerge away from the origin (as E J is increased) all lie along specific angles [14]: sin(p(θ + π/2)) = 1. For these angles most of the derivatives in equation (D.2) vanish, only those involving exactly one derivative with respect to θ remain. Setting both the drift and the determinant of the Jacobian to zero we find that the fixed point amplitude at the bifurcation satisfies: and the threshold value of E J at which the bifurcation occurs is This gives the threshold E J at which the FP drift terms show a bifurcation. Notice, however, that the E J at which the corresponding bifurcation occurs using the coherent state anstaz (see Sec. 3) drifts is slightly higher, by the factor e ∆ 2 /2 . It is these values (including the extra factor) that are plotted by the white circles in figure 2. Interestingly, expressions of the form of (D.3) apparently arise quite often [52].

Appendix E. Operator Bases
QuTiP was used to find the eigenoperators and eigenvalues of the Liouvillian [30]. The symmetry under R of these eigenoperators was checked and it was found that, for any eigenoperator with a unique L-eigenvalue λ i , the symmetry was obeyed with rotation eigenvalue either 1 or −1 (k n = 0 or π, the two unpaired modes with no direction of circulation around phase space). The remainder of the numerically discovered eigenoperators occur in sets that are degenerate with respect to the Liouvillian (within numerical precision), which means that any linear combination of eigenoperators from a given set is an equally valid eigenoperator. Thus the operators found numerically are valid, but in the wrong basis to have the rotational symmetry we expect. In order to find a basis of operators that possess this symmetry we first use Gram-Schmidt orthonormalization to find orthonormal operators that span the same space as each of these degenerate sets. Specifically we isolate a set of operators that are L-degenerate to within our numerical precision, {ρ num }, from this we generate the elements of a new basis {ρ GS } by: . (E.1) These basis operators will still not have the required symmetries, but they will be orthogonal allowing us to proceed to the next step, where the rotation super-operator R is determined as a matrix in this GS basis. Writing some linear combination of operators aρ GS 1 + bρ GS 2 + cρ GS 3 .. as a vector B = a b c . . . we see: The eigenvalues of this matrix representation of R are a subset of the eigenvalues of R (specifically the subset that is degenerate with respect to L). The eigenvectors of this matrix represent the symmetry basis-operators in terms of the Gram-Schmidt basis. Thus we have converted the numerical eigenoperators into a new basis where each eigenoperator of L is also an eigenoperator of R. These are the operators whose WDs are depicted in figures 6 and 7.

Appendix F. Physical rates
This Appendix presents the equations linking the physical rates in our effective description of the system dynamics to the Liouvillian eigenvalues. Initially specialising to the regime between the two bifurcations, where there are only p BPs we can express the point-to-point coupling with two differential equations: with A m the amplitude at the m th BP and D that at the DP. We now find the initial choices of D(0), A m (0) that result in all these quantities decaying together as e −λt in response to these differential equations. They can be found by setting λ =Ȧ m /A m =Ḋ/D and re-arranging.
First assume that all A m s are the same and notice that D(0) = 1, A m (0) = −1/p satisfies this relation. This is the k = 0 eigenoperator and has decay rate λ = −pΓ out − Γ in .
Assuming that D(0) = 0 leads to p solutions, one for each way of choosing A m+1 = exp(ik)A m such that m A m = 0. These have decay times λ(k) = Γ l (2 cos(k) − 2) − Γ in . For these modes Γ l produces band dispersion, although for our system we find Γ l to be very small. These solutions correspond to the eigenoperators depicted in figure 6.
Finally we generalise to the case after the second bifurcation, where each BP splits into two. To retain the previous rate definitions we now define A m to be the amplitude across a BP pair. Assuming the pair to be in phase leaves all solutions completely unchanged from the previous (pre-bifurcation) results.
There are now new solutions, corresponding to assuming the two BPs in the pair have opposite phases. Here two things are different. First a new internal decay rate exists, where amplitude from each of the BPs forming the pair annihilates with that from the other. Second, the interactions with the neighbouring pairs acquire a "−" sign, as each BP has an extra −1 factor relative to its neighbours. The new set of solutions have decay rates λ(k) = Γ l (−2 cos(k) − 2) − Γ in − 2Γ c , (these new ones include k = 0).
These relations are exploited in reverse to fit Γ in , Γ out , Γ l and Γ c to the known eigenvalues of the Liouvillian.