Emergent macroscopic bistability induced by a single superconducting qubit

The photon blockade breakdown in a continuously driven cavity QED system has been proposed as a prime example for a first-order driven-dissipative quantum phase transition. But the predicted scaling from a microscopic system - dominated by quantum fluctuations - to a macroscopic one - characterized by stable phases - and the associated exponents and phase diagram have not been observed so far. In this work we couple a single transmon qubit with a fixed coupling strength $g$ to an in-situ bandwidth $\kappa$ tuneable superconducting cavity to controllably approach this thermodynamic limit. Even though the system remains microscopic, we observe its behavior to become more and more macroscopic as a function of $g/\kappa$. For the highest realized $g/\kappa \approx 287$ the system switches with a characteristic dwell time as high as 6 seconds between a bright coherent state with $\approx 8 \times 10^3$ intra-cavity photons and the vacuum state with equal probability. This exceeds the microscopic time scales by six orders of magnitude and approaches the near perfect hysteresis expected between two macroscopic attractors in the thermodynamic limit. These findings and interpretation are qualitatively supported by semi-classical theory and large-scale Quantum-Jump Monte Carlo simulations. Besides shedding more light on driven-dissipative physics in the limit of strong light-matter coupling, this system might also find applications in quantum sensing and metrology.


I. INTRODUCTION
Quantum phase transitions (QPT), both first-and second-order [1] have been at the forefront of physics research for half a century.The original idea of QPTs as abrupt shifts in the (pure) ground state of closed quantum systems as a function of a control parameter applied mostly to condensed matter physics.Dissipative phase transitions (DPT) occurring in the (in general, mixed) steady state of open quantum systems [2][3][4][5][6][7][8][9][10][11][12], however, broadened the scope of phase transitions to encompass mesoscopic and later even microscopic systems, where the interaction with the environment essentially affects the system dynamics.A DPT was first realized experimentally in a Bose-Einstein condensate interacting with a single-mode optical cavity field [13], and DPTs are increasingly relevant to today's quantum science and technology [14][15][16][17].
In view of this success, it is remarkable that in recent years yet another phase-transition paradigm could emerge, namely, first-order dissipative quantum phase transitions.A first-order phase transition means that two phases can coexist in a certain parameter region, like water and ice at 0 °C for a certain range of free energy.Coexistence of phases in the quantum steady state seems paradoxical, since the steady-state plus normalization conditions for the density operator constitute a linear system of equations, that admits only a single solution.That is, given the Liouvillian superoperator L for the Markovian evolution of the system, there exists only a single normalized density operator ρ st that satisfies Lρ st = 0. ( The resolution is that a single density operator can accommodate the mixture of two macroscopically distinct phases expressed as a ratio of the two components.In the water analogy, at 0 °C we could symbolically write with c growing from 0 to 1 as the free energy is increased. Recently, first-order dissipative quantum phase transitions have been found in various systems.One such platform is the clustering of Rydberg atoms described by Ising-type spin models [18][19][20][21][22][23] and realized experimentally [24][25][26].Various other systems of ultracold atoms [27,28] and dissipative Dicke-like models [29,30] also exhibit signatures of a first-order DPT.Other platforms include (arrays of) nonlinear photonic or polaritonic modes [7,[31][32][33][34][35][36][37][38], exciton-polariton condensates [39,40] and circuit QED [16,[41][42][43].In this work we observe and model the scaling and phase diagram of a first-order DPT in zero dimensions, i.e. for a single qubit strongly coupled to a single cavity mode.

II. PHOTON-BLOCKADE BREAKDOWN
The Jaynes-Cummings (JC) model -one of the most important models in quantum science -describes the interaction between atoms and photons trapped in a cavity [44].It is expressed by the Hamiltonian ( = 1) with ω R the angular frequency of the cavity mode with boson operator a, ω A that of the atomic transition with operator σ, g the coupling strength, η the drive strength, and ω the drive frequency.This model yields the prototype of an anharmonic spectrum in the strong-coupling regime, as demonstrated in cavity [45] and circuit QED [46], and with quantum dots in semiconductor microcavities [47].Its strong anharmonicity at single photon levels is the basis of the photon blockade effect [48,49], in analogy with Coulomb blockade in quantum dots or to polariton blockade [50].Photon blockade means that an excitation cannot enter the JC system from a drive tuned in resonance with the bare resonator frequency, or similarly, a second excitation from a drive tuned to resonance with one of the single-excitation levels cannot enter.
This blockade is, however, not absolute, as it can be broken [51][52][53][54] by strong enough driving due to a combination of multi-photon events and photon-number increasing quantum jumps [55].In an intermediary η range, in the time domain the system stochastically alternates between a blockaded, dim state without cavity photons and a bright state in which the blockade is broken and the system resides in the highly excited quasi-harmonic part of the spectrum resulting in a large transmission of drive photons.In phase space, this behavior results in a bimodal steady-state distribution in analogy with Eq. ( 2), with c growing from 0 to 1 with increasing η.This effect has been demonstrated experimentally in a circuit QED system [42].Bistability in the time domain or bimodality in phase space is, however, not sufficient evidence for a first-order phase transition.It is also necessary that the two constituents in the mixture Eq. ( 4) corresponding to the two states in the temporal bistable signal to be macroscopically distinct as is the case in Eq. (2).It has been shown theoretically [51,55], that the photon blockade breakdown (PBB) effect has such a regime, i.e. a thermodynamic limit, where both the timescale and the amplitude of the bistable signal goes to infinity, resulting in longlived and macroscopic distinct dim and bright phases.Remarkably, this thermodynamic limit is a strong-coupling limit, defined as g → ∞, and independent of the physical system size, i.e. the system remains the same JC system composed of two microscopic interacting subsystems.In this limit, the temporal bistability is replaced by hysteresis, where the state of the system is determined by its initial condition, since switching to the other state entails an infinite waiting time.The passage to the thermodynamic limit, i.e. the indefinite increase of g has been termed finite-size scaling [55].
In this work, we demonstrate these additional criteria that clearly signify the observed physical effect as a first-order dissipative quantum phase transition.We demonstrate the finite-size scaling over 7 orders of magnitude towards the thermodynamic limit and back out the phase diagram of a first-oder DPT in zero dimen- sions.We realize this experiment with a superconducting qubit strongly coupled to a bandwidth-tunable microwave cavity mode and find qualitative agreement with large-scale Quantum-Jump Monte Carlo simulations and semi-classical calculations of the phase boundaries.

III. EXPERIMENTAL REALIZATION
Our experimental setting incorporates a transmon qubit [56,57] placed at the anti-node of the standing wave of a 3D copper-cavity, as shown in Fig. 1(a), that can be flux-tuned by applying a B-field via a millimetersized superconducting bias coil mounted at the outside cavity wall.The transmon qubit has a maximum Josephson energy E J,max /h ≈ 48 GHz, charging energy E C /h ≈ 382 MHz and a resulting maximum transition frequency between its ground and first excited states of ω A /2π ≈ 12.166 GHz.When the transmon ground to first excited state transition is tuned in resonance with the cavity mode at ω R /2π ≈ 10.4725 GHz, the directly measured coupling strength between the single photon and the qubit transition is as high as g/2π = 344 MHz, which is only about a factor of 3 below the so-called ultrastrong coupling regime [58].The relatively high absolute anharmonicity between subsequent transmon state transitions is α/h ≈ −418 MHz at this flux bias position.
The cavity has two ports, of which the input pin coupler position is fixed with an external coupling strength of κ fixed /2π ≈ 500 kHz.The output coupler is attached to a cryogenic piezo nano-positioner, which allows for adjusting the pin length extending into the cavity [59].With this tunable coupler the coupling strength can be varied in situ in a wide range κ vary /2π ≈ 20 kHz − 30 MHz.The internal cavity loss at low temperature is κ int /2π ≈ 600 kHz, which is achieved by electro-polishing of the high conductivity copper surface before cooldown to 10 mK in a dilution refrigerator.
All four scattering parameters are measured with a vector network analyzer to calibrate the measurement setup and the cavity properties when the qubit is far detuned from the cavity resonance.
From these fits, we extract all loss rates that add up to the total cavity linewidth κ = κ fixed + κ vary + κ int also indicated in Fig. 1(b).Time-domain characterization measurements confirm that the qubit is Purcell-limited and homogeneously broadened at the flux sweet spot [61], where the measured coherence times are T 1 ≈ 0.5 µs and T 2 ≈ 1 µs.When the qubit frequency is tuned far below the resonator frequency ω A /2π ≈ 6.083 GHz by applying an external magnetic field, the measured coherence times are T 1 ≈ 18.14 µs and T 2 ≈ 0.496 µs, which we attribute to a higher Purcell limit due to the larger detuning as well as a drastically increased flux noise sensitivity.On resonance ω A = ω R , where the following experiments were performed, the energy relaxation is therefore fully dominated by cavity losses.The measured vacuum Rabi peak linewidth changes with and without the qubit in resonance are in agreement with a small amount of flux noise induced dephasing expected at that flux bias position.

IV. PHOTON BLOCKADE BREAKDOWN MEASUREMENT
The photon blockade (and its breakdown) phenomenon most straightforwardly occurs when the two interacting constituents are resonant ω A = ω R .In contrast to the ideal two-level atom limit [51,55], when driven on resonance ω = ω R this does not lead to spontaneous dressedstate polarization [62,63] -a second-order DPT [51], in our experimental situation with three (or more) transmon levels [42] as shown in Fig. 2(a).For low input powers corresponding to less than a single intra-cavity photon on average we observe a vacuum Rabi-split spectrum in transmission, as shown in Fig. 2(a, b) (blue line).No transmission peak is observed at the bare cavity frequency ω R up to intermediate input drive strengths η.This means that a single photon -or even hundreds of photons at the chosen g/κ = 39.1 -are prevented from entering the cavity due to the presence of a single artificial atom.This blockade is observed to be broken abruptly by further increasing the applied drive strength η, which is proportional to square-root of the applied drive power and the corresponding drive photon number.As η is increased by only a finite amount, the transmitted out-put power increases by three orders of magnitude at the bare resonator frequency, as shown in the red spectrum in Fig. 2(b).The central sharp peak in the transmission spectrum corresponds to a time-averaged measurement (determined by the chosen resolution bandwidth) of a cavity that is fully transparent for most of the integration time.This PBB effect can be attributed to the nonlinearity of the lower part of the JC spectrum which is strongly anharmonic [46,64], while the higher-lying part of the spectrum has subsets that are closely harmonic over a certain range of excitation numbers [65] and can hence accommodate a closely coherent state.
In the time domain, with η in the phase coexistence region, the PBB effect results in a bistable telegraph signal, where the system output alternates between a 'dim' state where the qubit-resonator system remains close to the vacuum state unable to absorb an excitation from the externally applied drive, and a 'bright' state where the system resides in an upper-lying, closely harmonic subset of the JC spectrum, cf.Fig. 2(c).The switches between these two classical attractors are necessarily multiphoton events that are triggered by quantum fluctuations.This bistability was shown to be a finite-size precursor of what would be a first-order DPT in the thermodynamic limit (g/κ → ∞) [55], where the bistability develops into perfect hysteresis: the system is stuck in the attractor determined by the initial condition as long as the control parameters are set in the transition domain.
In order to investigate this dynamics qualitatively, we record the real-time single-shot data of both quadratures of the transmitted output field at the bare cavity frequency while applying a continuous-wave (CW) drive tone resonant with the bare cavity, over a range of applied drive strengths.The transmitted radiation is first amplified with a high electron mobility transistor (HEMT) at 4 K followed by another room-temperature low-noise amplifier (LNA), then down-converted with an IQ mixer with appropriate IF frequency and finally digitized with a digitizer.Further this recorded data is digitally lowpass filtered with appropriate resolution bandwidth and down-converted to d.c. to extract the time-dependent quadratures in voltage units.For example, in the case of κ/2π = 8 MHz, the recorded data is 2.88 s long and the final time resolution of the extracted quadratures is 2.5 µs, cf.Fig. 2(c).The selection of an appropriate resolution bandwidth is critical for a number of reasons: (1) to successfully resolve frequent and sudden switching events caused by very short dwell times at high κ values, (2) to maintain a signal to noise ratio that allows to clearly discriminate single shot measurement events without averaging, and (3) to achieve a sufficient total measurement time to resolve long dwell times with the available memory.
From the resulting histograms in phase space, cf.Fig. 2(d-f), which represent the scaled Husimi-Q functions convolved with the added amplification chain noise photon number n amp ≈ 9.2, it can be deduced that for low drive strength the photon blockade is intact (dim phase) with the Q function being centered around the vacuum state.Upon increasing the input drive strength, the Q function becomes bimodal with decreasing weight of the dim state as described in Eq. 4. At high enough η only the bright coherent state is measured.Note that the transformation of the Q-function and hence the steadystate density operator of the system as a function of η is continuous, yet a first order phase transition with a well-defined coexistence region can occur.
A similar conclusion can be drawn from the output power histograms (color map) that trace out a typical bistability curve as shown in Fig. 2(g).The most likely output powers P out and calculated equivalent intra-cavity photon numbers of the empty cavity driven on resonance ncav = Pin ωR as a function of applied input power P in and resulting drive strength η = √ ncav κ/2 are marked with circles.The vacuum, bistability and bright regions are well defined.We find that the derivative of the bright solution obtained at high P in deviates somewhat from the empty cavity response measured when the qubit is far detuned (dashed line).For large g/κ values this is more pronounced and we have observed that this can lead to secondary bistability regions at even higher powers for g/κ 43, which we believe to originate from the multi-level nature of the transmon qubit.In this work we focus only on the bistability occurring at the lowest drive strength.
In Fig. 2(h) we show the measured dwell times of dim and bright states as a function of input power.The average dwell times at each input power are calculated as and the threshold for one of the N switching events during the full measurement duration (with 2.5 µs resolution) is defined at half of the observed full amplitude for the lowest applied η/(2π) ≈ 140 MHz where the bistability is fully developed.Note that in cases of low signal to noise ratio, e.g. for large κ or at high drive detunings shown later, we used a higher threshold that at least exceeds the variance of the output power of the dim state.At low η the statistics is low because the system remains in the dim state for long time scales.As η is increased to η * /(2π) ≈ 167 MHz the measured average dwell times cross at t dim dwell = t bright dwell = t * dwell ≈ 354 µs.We call this the half-filling point where it is equally likely for the system to be found in the dim or bright state, denoted by an asterisk in Fig. 2(h).For even higher drive strength the system prefers to dwell in the bright state, i.e. we observe close to full resonator transmission for most of the time.

V. FINITE SIZE SCALING TOWARDS THE THERMODYNAMIC LIMIT
We set out to experimentally prove the scaling of the measured PBB bistability towards a thermodynamic limit, i.e. to show -despite the underlying microscopic nature of the system -a truly macroscopic behavior as g/κ → ∞ where it has been shown to become a first-order DPT.We also experimentally determine the finite-size scaling exponents of the characteristic time and the corresponding drive strength and intra-cavity photon number for this DPT.
In Fig. 3(a) we show the measured dwell times as a function of drive strength -similar to Fig. 2(h) -for different κ/2π ranging from 18.1 MHz to 1.2 MHz.For each κ value the dwell time in the dim state (blue symbols) decreases with increasing drive strength and that of bright state (red symbols) increases until eventually the system is fully stabilized in the bright state.For each g/κ value, we define a single characteristic time of the process t * dwell at the drive strength η * that leads to half-filling of the telegraph signal, cf. also Fig. 2(h).
Remarkably, for the largest realized value of g/κ = 287, the characteristic dwell time reaches t * dwell ≈ 6 s.This exceeds the characteristic microscopic dissipation times of order κ −1 by a factor ∼ 10 6 and is reminiscent of the emergence of two macroscopically distinct states with strongly suppressed transitions that require a cascade of quantum jumps [55].This strong coupling, high photon number limit where the effect of quantum fluctuations vanishes has been defined as the thermodynamic limit of such a finite-size zero-dimensional system [51].
Finite-size scaling towards the thermodynamic limit means that the characteristic time scales and brightness scale up as a function of g/κ, while the system remains self-similar, which in this case means that it keeps switching stochastically at a fixed filling factor (that we choose 0.5).In Fig. 3(b) we plot the measured increase of t * dwell at filling 0.5 over a range of seven orders of magnitude as a function of g/κ.The behavior follows a strong power law over the full range and the fitted finite-size scaling exponent is t * dwell ∝ (g/κ) 5.4±0.2 .Similarly, the intracavity photon number of the bright state at half-filling increases nearly linearly with n * ∝ (g/κ) 0.96±0.05, as shown in Fig. 3(c), and as a consequence the corresponding drive strength decreases with the square root η * /g ∝ (g/κ) −0.52±0.03 .Here the drive normalization with g is motivated by the 2-level neoclassical theory, where the critical point appears at η/g = 0.5 for ∆ = 0.
The theoretical results (triangles, pentagons and heptagons) shown in Fig. 3(b)-(d) are taken from large scale numerical simulations performed with C++QED: a framework for simulating open quantum dynamics [66].An adaptive version of the Quantum-Jump Monte Carlo (QJMC) method is applied, where a single stochastic quantum trajectory is considered to correspond to a single experimental run [67].The shown data is based on 64 CPU years of simulation time with a Hilbert space dimension of up to ∼ 2 14 − 2 15 (7 transmon levels and 3-5 times n * ).Another computationally demanding aspect is that the required time step is set by the largest characteristic frequency (typically g or η) of the microscopic system (sampled with 1/κ to reduce data volume) whereas the total trajectory needs to cover many times t dwell to obtain sufficient statistics of the macroscopic behaviour.Together this limits the range of numerically accessible g/κ values to the lowest three values investigated.For more details on how we model the system, examples of simulated quantum trajectories and the impact of different transmon dephasing models, cf.Appendix A.
The observed photon number scaling exponent is about half of the analytical prediction of n * ∝ (g/κ) 2 .This is not surprising since the 2-level neoclassical theory does not yield quantitative agreement for the case with a multi-level transmon circuit.Numerical simulations for the lowest three g/κ values taking into account up to 7 transmon levels (heptagons in Fig. 3(c)), agree with the measured linear exponent to within 15%.The absolute value of n * also agrees well (blue) and is further improved when we include qubit dephasing of all transmon levels (red) due to flux noise with the measured γ φ /2π ≈ 50 kHz for the lowest qubit transition.That is, as long as enough transmon levels are taken into account.In the case of just 3 transmon levels the simulated value is about an order of magnitude smaller than the measured one.This highlights the importance and participation of multiple transmon levels in the dynamics of the system.
The dwell time values and scaling shown in Fig. 3(c) g/κ at half filling (asterisks).Fitting the experimental data (dashed lines) yields the exponents t * dwell ∝ (g/κ) 5.4±0.2 for the dwell time, n * ∝ (g/κ) 0.96±0.05for the average resonator photon number and η * /g ∝ (g/κ) −0.52±0.03for the input drive strength.We compare the experimental values (stars) with those obtained from QJMC simulations (polygons with color code) for computationally manageable κ = 11.5, 14.3, and 18.1 MHz.The simulations include 3, 5, or 7 transmon levels both, without transmon dephasing (light blue) and with dephasing increasing linearly with the level number (light red) with γ φ,1 /(2π) = 50 kHz and γ1 = 0. Details of the simulation are in the main text and in Appendix A.
are more robust with regards to the number of transmon levels but we observe a substantial deviation between the measured (5.4) and simulated scaling exponents in the range of 0.9-2.2.In [55] based on a two-level model with finite detuning, the blink-off rate could be calculated from the rate of ladder-switching quantum jumps, and was found to be proportional to κ/ n * , so the waiting time for a blink-off is n * /κ, therefore it scales as g 2 /κ 3 .The numerically determined timescale-exponent in the same work was (κt * dwell ) ∝ (g/κ) 2.2 , which is very close to the analytical value.
The observed scaling exponent is significantly larger which could be due to counter-rotating terms not taken into account in the numerical simulations and due to hybridized decay channels which can invalidate our approach of using separate transmon and resonator decay channels in the master equation.Other potentially participating mechanisms could include transmon ionization [68] or dielectric surface loss saturation [69] that might further stabilize the system -in particular in the bright attractor and thus reduce the blink-off probability.Taken together these effects appear to be driving the system significantly faster to the thermodynamic limit compared to what would be expected from the standard Jaynes-Cummings model.
Importantly and irrespective of the origin of the unexpectedly strong scaling, we observe that the system is always able to relax to the vacuum state eventually -despite the continuous driving.And this vacuum state is then stabilized for seconds by the presence of a single qubit, even for the highest drive strength corresponding to a photon number of n ≈ 10 4 .In contrast to similarly looking fluorescence signals with dwell times on the order of seconds, which have been known in quantum science since the famous first electron-shelving experiments with single-ions in Penning-traps [70,71] due to long lived metastable atomic states, in the present case the measured time scale exceeds all microscopic time scales by up to a factor 10 6 , i.e. the system is very deep in the macroscopic limit which justifies its classification as a finite-size phase transition.

VI. PHASE DIAGRAM
The PBB phase transition has been predicted to give rise to an interesting phase diagram as a function of drive detuning ∆ ≡ ω − ω R = ω − ω A .From the neoclassical equations, valid for pure two level qubit states and γ 1 = γ φ = 0, the bistable region is only expected at finite detunings away from the critical point at ∆ = 0 [51,55].A more realistic prediction is obtained from numerical solutions of the semiclassical Maxwell Bloch equations that include qubit decay (for details, cf.Appendix B), which are shown in Fig. 4(a),(b) for a 2-and a 3-level transmon, respectively.The results are qualitatively different not only in comparison to the neoclassical solution but also between the cases with 2 and 3 level transmons.Most importantly, in the case of 3 levels the bright phase boundary forms a peak rather than a dip at zero detuning and the dim phase is also predicted to exists around zero detuning.
Experimentally we choose a large g/κ ≈ 132.3,where the time scales are long and the phases are very well defined, to back out the phase diagram as a function of ∆.We sweep the drive strength for each chosen detuning and record a trace of single shot time domain transmission data for each parameter combination.The result is shown in Fig. 4(c) and (d), which shows the three regions traced out on the ∆−η parameter plane and the measured dwell time at half filling (dashed lines), respectively.
The phase boundaries (points) are obtained from measured time-domain single shot telegraph transmission data as shown in Fig. 4 (e)-(h) for the range of η at ∆ = 0 indicated in panel (a) with a double arrow.Here we define a threshold (dashed line) as described earlier and count if a single phase switching attempt was successful to cross this threshold within the measurement time.If the answer is yes the corresponding η and δ value pair is assigned to the bistable region boundary.If the answer is no, depending on the measured value (below or above threshold) the parameter combination is assigned to the dim or bright phase.The parameter region where multiple crossings occur is assigned to the bistable region.For each detuning the detection bandwidth and total measurement time has been optimized to be able to determine a sharp phase boundary and to be able to resolve the dwell time over 5 orders of magnitude, as shown in Fig. 4(d).
The raw data in Fig. 4 (e)-(h) reveals an interesting difference between partial phase switching attempts from the dim state, which are quite frequent; and from the bright state, which are rather rare and typically of smaller amplitude (not visible in this data).This asymmetry is not observed in the simulated quantum jump trajectories as shown e.g. in Fig. 5, and its origin is not clear.However, these data point at an additional stabilization mechanism of the bright phase that might also contribute to the stronger than expected scaling towards the thermodynamic limit.
Comparing the theoretical and experimental phase diagrams (a-c), one can see that the semiclassical 2-level model (similar to the neoclassical model) fails to capture the essential features of the experiment.The agreement with the 3-level case is qualitatively much better especially in case of the bright phase boundary.However, the shape of the lower limiting curve is not correctly captured by this theory.In fact, the 3-level semiclassical theory is unable to reproduce the well pronounced dim phase region around resonance that we observe experimentally.In general it reproduces the overall resonance-like dependence on the detuning, and also the asymmetry with respect to the ∆ = 0 line, but no quantitative agreement in terms of shape or absolute values could be obtained.This comparison with the semiclassical theory underlines that the well-resolved spectrum of the strongly coupled transmon-resonator system with more than two transmon levels plays an essential role in our experiment. Δ/κ

VII. DISCUSSION, CONCLUSIONS AND OUTLOOK
It is important to distinguish the presented PBB phase transition and scaling from similar related phenomena.The oldest known such effect is optical bistability, dispersive or absorptive, that is itself a first-order DPT [31,72,73].In the case of the PBB, we are not in the dispersive regime however.The driving is close to or on resonance with the bare resonances of the resonator and the transmon, and the absorption of the latter does not play an essential role either [55].Another related model is the Duffing oscillator that appears in a circuit QED context as a Kerr-nonlinear mode (the transmon) interacting with a linear mode (the resonator) [74,75].Parametric driving can lead to critical behavior [76] and driven nonlinear inductors have exhibited slow classical switching events triggered by low frequency thermal fluctuations on the order of seconds [77].Long bit flip times up to 100 seconds have also been observed in a two-photon dissipative oscillator that is characterized by symmetry breaking of the intra-cavity field phase, but this system does not exhibit a bistability in the photon number [78].In this respect, the transmon, even with many levels considered, is algebraically very different from a nonlinear oscillator when it comes to jump operators because these are not bosonic.This makes an essential difference, as verified by our quantum simulations, where it is possible to try the consequences of different algebras.Our simulations clearly rule out the Duffing oscillator model, which cannot reproduce the phenomenology of the experiment since its bistable behaviour reminiscent of dispersive optical bistability occurs for different parameters and does not exhibit the same scaling towards the thermodynamic limit.Finally, with respect to other recently-discovered QPTs and DPTs in the Jaynes-Cummings or Rabi models [9,[79][80][81], where thermodynamic limits can also be defined in an abstract way, the difference of PBB as firstorder DPT is that the thermodynamic limit is a strongcoupling limit.The well-resolved discrete spectrum of an interacting bipartite quantum system is essential for the effect.
In this paper, we have experimentally followed the finite-size scaling towards the g/κ → ∞ thermodynamic limit with a characteristic time scale ranging over nearly seven orders of magnitude.Just like with a finite-size (non-macroscopic) sample of water, at 0 °C there is a contest of several meta-and even unstable states instead of true phases of liquid and ice [82], in the PBB bistability for any finite value of g, there is a contest of nonmacroscopically distinct dim and bright states.We have experimentally determined the scaling exponent of the bistable switching timescales of ∼ 5.4 ± 0.2 as well as the finite-size scaling exponent of the intra-cavity photon number of the bright state of ∼ 0.95 ± 0.05.We have also experimentally determined the phase diagram of the PBB phase transition and found that the characteristic dwell times drop by orders of magnitude for finite drive detunings.
We have compared these experimental results with large scale quantum simulations based on the QJMC method considering different numbers of transmon levels (3, 5 and 7) and different dephasing models of higherlying levels.This comparison indicates that transmon levels up to at least 7 play an important role in the dynamics, as does the dephasing, since simulations with dephasing, e.g.due to flux noise, have shown better correspondence to the experimental data.Similarly, the comparison with semiclassical phase diagrams have also underlined the important role of higher transmon levels.Even though the full quantum simulations reproduce the observed trends correctly, there are significant differences from the experiment in the measured dwell times, which might indicate the presence of further stabilization mechanisms -in particular in the bright attractor -and calls for a better methods to model such strongly coupled multi-level systems.
Even though the computational resources were substantial, the fully quantum numerical simulations were only suitable to model the three lowest coupling strengths g/κ investigated.This highlights the need for powerful quantum simulators even in the case of comparably simple circuits and in particular to explain how macroscopic phases can be stabilized by individual quantum systems.It is in fact quite surprising that a single transmon qubit can switch back from the bright state -characterized by up to 10 4 intra-cavity photons -all the way to the dim state and stabilize the empty cavity for seconds in the presence of the continuous large amplitude coherent input field -in particular given its limited potential confinement [68].In the future, a fully confined qubit [83] with higher power handling, or larger anharmonicity and superconducting cavities with lower loss could help to explore even more macroscopic phases pushing the characteristic switching timescales from seconds to days.
Besides its fundamental interest as a quantum-classical phase transition, the PBB bistability also promises a few applications.Since single quantum jumps were shown to trigger the switching from the (nonclassical) dim state to the (closely classical) bright state [55], our system may be considered as a quantum-jump amplifier, where ultimately a macroscopic microwave device (out-side the fridge) is getting switched by microscopic quantum events (inside the fridge).An interesting prospect is controlling the switching behavior, that can be envisaged either in a parametric way, but preferably with another strongly coupled quantum system.In the latter case the bistability could act as a quantum readout device with high signal-to-noise ratio.The capability of preparing the system on the verge of a phase switching event could therefore make it applicable in quantum metrology and sensing based on microwave photon counting [84], a new paradigm for the application of first-order DPTs [85][86][87][88][89].
The data and code used to produce the figures in this manuscript will be made available on the Zenodo repository.

H =
) Here, putting h 0 = 0, and assuming the 0−1 transition resonant with the mode (h 1 = ω R ), we obtain a simple form for the bare transmon Hamiltonian, which we list for the first 3 levels: where ∆ an ≡ h 2 − 2h 1 is the anharmonicity of the third level, which is related to the charging energy.
Let us turn to dissipation, which we describe with the Liouvillian with the following three dissipative channels: Resonator decay, L mode This is described by the jump operators L − = 2 (n th + 1) κ a and L + = √ 2 n th κ a † .Here n th is the number of thermal photons, which can be neglected in our system, so the second kind of quantum jumps (absorption of thermal photons) does not exist.
Energy relaxation of the transmon, L relax In analogy with the coupling to the resonator mode, we assume that this occurs only as transitions between adjacent levels.It is described by the jump operators L u+1→u = √ γ u+1→u |u u + 1|.In the simulation, we take γ u+1→u equal for all levels, and we identify it with γ 1 in cQED.
Dephasing of the transmon, L dephase This is also defined separately for all transmon levels, and its jump operator for level v is so it simply flips the phase of level v by π.Modeling the behavior of the dephasing for different transmon levels is nontrivial.We consider three possibilities: 1. γ φ,v = 0 for all v.This is only to get a theoretical baseline of dephasing-free behavior.
2. Linear growth as γ φ,v = v γ φ /8 following the above convention, as expected for flux noise due to the higher flux gradient of higher levels.Here the dephasing of the first qubit transition γ φ /2π = 50 kHz is taken from the measured vacuum Rabi linewidths and the independently measured cavity linewidth κ.
3. Dephasing proportional to the charge dispersion of the transmon levels [64].
Example trajectories for the three possibilities are displayed in Fig. 5.It is apparent that model 3 leads to very noisy trajectories that do not reproduce qualitatively the experimentally observed behavior of stabilized attractors.Therefore, we omitted this possibility from the quantitative comparison presented in the main text.
In the simulation, for each physical parameter set, several trajectories are run with different random number generator seeds.Relying on the assumption of ergodicity, these trajectories are concatenated for a single long trajectory for each parameter set, which is then used for dwell-time statistics.Since each trajectory is started from the ground state, this method has a bias toward the dim state (breaching of ergodicity), which is the stronger, the larger the dwell time with respect to the simulation time.
The full quantum simulations were implemented within the C++QED simulation framework (http:// github.com/vukics/cppqed),and took about a year on a 64-core virtual cluster defined within an OpenStack Cloud environment (http://science-cloud.hu/).From the master equation ρ = [H, ρ]/(i ) + Lρ we can derive equations for the expectation values of the operators a and σ uv = |u v|.In the case of a two-level system, this simply reproduces the Maxwell-Bloch equations, with the added complication of the qubit dephasing.Here, we list the equations for a three-level transmon with states |g , |e and |f , which still leads to an algebraically tractable scheme.In this case, 6 equations are needed for a complete system: α =(i∆ − κ)α + η − g 1 s ge − g 2 s ef (B1a) ṡge =(i∆ − γ 1 ) s ge − g 1 (s ee − s gg ) α − g 2 s gf α * (B1b) ṡgg =γ 1 s ee − 2g 1 {α * s ge } (B1c) ṡef =(i[∆ − ∆ f ] − [γ 1 + 4(γ φ,1 + γ φ,2 )]) s ef + g 2 (s ee − s f f )α + g 1 α * s gf (B1d) ṡee =2 g 1 {α * s ge } − 2 g 2 {α * s ef } − γ 1 s ee + γ 1 s f f (B1e) ṡgf =(i[2∆ − ∆ f ] − [γ 1 + 4γ φ,2 ]) s gf − g 1 α s ef + g 2 α s ge (B1f) Here α = a , and s uv = σ uv .The system is completed with the completeness relation s gg + s ee + s f f = 1.We are interested in the steady state, which can be obtained by zeroing the left hand side of the equations, that leads to an inhomogeneous nonlinear set of equations.We do not need to solve the full set of equations.Instead, we can obtain a single implicit equation for only the intensity |α| 2 as follows.First, we define the complex dispersive shift then, from Eq. (B1a) in steady-state we express α explicitly.As we will show below, Σ depends only on powers of |α| 2 , and not on other combinations of α and α * .Therefore, the equation for the intensity can be written as What we have to show for the validity of Eq. (B3) is that the solutions of s ge and s ef have the form of an |α| 2 -dependent expression multiplied by α.The polarizations can be expressed as functions of the populations multiplied with α from Eqs. (B1b), (B1d) and (B1f) in steady state.When these solutions are substituted into the steady-state population equations (B1c) and (B1e), the factor α in the solutions together with α * in those equations give an |α| 2 .Hence, the populations can be expressed from these equations as functions only of the intensity, and when these are substituted back into the solutions of the polarizations, we obtain the necessary form for these latter.
A typical solution of Eq. (B3) exhibiting bistability is displayed in Fig. 6.The semiclassical theory is inferior to the full quantum-trajectory solution described in the main text and Appendix A in at least two respects: 1. Dealing with (possibly multi-valued) steady state solutions, it does not provide information on timescales.
2. Whereas the set of three complex polarizations and two populations in Eq. (B1) give a complete picture of the transmon in itself, the mode is represented only by a single amplitude.This means that the theory cannot account for nonclassical states of the mode and transmon-mode entanglement.

FIG. 1 .
FIG. 1. Experimental realization.(a) Schematics of the experimental device consisting of a superconducting transmon qubit fabricated on a silicon substrate that is placed at the antinode of the fundamental mode of a 3D copper cavity.The cavity has a fixed length port (red) and an in-situ variable length pin coupler port (blue).(b) Measured cavity transmission spectra with the qubit far detuned for different coupler positions (color coded) together with a fit to Eq. (5) (dashed) and the extracted κ/2π.

Figure 1 (
b) shows transmission measurements fitted with the scattering parameter S 21 derived from the Input-Output theory of an open quantum system[60]

1 FIG. 2 .
FIG. 2. Observation of photon blockade breakdown at g/κ ≈ 39.(a) The Jaynes-Cummings ladder for a three-level atom illustrating the PBB effect in the frequency domain: single-photon (blue) and multi-photon transitions (red) are indicated according to the measured spectrum at the Rabi-split frequencies and near resonance, respectively.(b) Measured cavity transmission spectra for ωA = ωR for two applied external drive strengths 1.05 MHz revealing a typical vacuum Rabi spectrum and 167.01 MHz where a sharp peak at ωR is observed.(c) Measured cavity output bistability at ωR in the time domain indicating the dwell times of the bright (ton) and dim states (t off ). (d-f ) Measured quadrature histograms (proportional to cavity Q-functions convolved with amplifier noise) for the dim phase at η/2π = 105 MHz, for the bistable region with equal probability at 167 MHz, and for the bright phase at 210 MHz, respectively.(g) Measured histograms of the output power (arranged vertically and probability is color coded) as a function of input power.The maxima indicated by circles trace out a typical bistability curve, cf.Fig. 6 in Appendix B. (h) Extracted average dwell times in the dim and bright state (Eq.(6)) as a function of η.The error bars represent the standard error that is extracted from 5 sections of the full data set.The dwell-time and drive strength corresponding to half-filling where t dim dwell = t bright dwell are indicated with an asterisk.

FIG. 3 .
FIG. 3. Finite size scaling towards the thermodynamic limit.(a) Measured dwell times in the dim (blue) and bright states (red) are shown as a function of η/g for different values of κ.Black stars denote half-filling values where the probabilities to be in the dim or bright state are equal.The dashed line is a power law fit.The error bars represent the standard error that is extracted from at least 3 sections of the full data set.(b)-(d) Measured scaling laws as a function of control parameter g/κ at half filling (asterisks).Fitting the experimental data (dashed lines) yields the exponents t *dwell ∝ (g/κ) 5.4±0.2 for the dwell time, n * ∝ (g/κ) 0.96±0.05for the average resonator photon number and η * /g ∝ (g/κ) −0.52±0.03for the input drive strength.We compare the experimental values (stars) with those obtained from QJMC simulations (polygons with color code) for computationally manageable κ = 11.5, 14.3, and 18.1 MHz.The simulations include 3, 5, or 7 transmon levels both, without transmon dephasing (light blue) and with dephasing increasing linearly with the level number (light red) with γ φ,1 /(2π) = 50 kHz and γ1 = 0. Details of the simulation are in the main text and in Appendix A.

FIG. 4 .
FIG. 4. PBB phase diagram at g/κ ≈ 132.Phase diagram on the ∆ − η plane obtained from semiclassical Maxwell-Bloch equations for (a) two and (b) three transmon levels with values for the parameters γ1 and γ φ according to the experiment.In panel (b), the noisy lower boundary of the bistable region is a result of numerical errors.(c) The PBB phase diagram with boundaries obtained from the experimental data (points) and the half-filling drive strength (dashed line).(d) Measured dwell times at half filling as a function of drive detuning.Error bars are extracted from the mean and standard error of the measured dim and bright dwell times.(e)-(h) show the experimentally observed telegraph signals that define the phase boundaries at ∆ = 0 in the range of η/2π from 85.6 to 115.5 MHz as indicated by a double arrow in panel (c).
Appendix B: The Maxwell-Bloch equations and their solution for the intensity

FIG. 6 .
FIG.6.Typical bistability curves obtained from Maxwell-Bloch equations for the intra-cavity intensity as a function of the drive strength for g/κ ≈ 132 and detuning ∆/(2π) = −10 MHz for in case of 2 and 3 transmon levels in panels (a) and (b), respectively (γ1/(2π) = 1 kHz, γ φ,1 /(2π) = 50 kHz).The bistable region is characterized by a threefold (for certain parameters even fivefold) solution but one (two) that show a decreasing intensity with increasing drive strength are unphysical.The empty-cavity solution (blue line) is shown for zero detuning.