Finite-size scaling of the photon-blockade breakdown dissipative quantum phase transition

The time evolution of a quantum system undergoing a dissipative, first-order quantum phase transition is studied. We consider the photon-blockade-breakdown phase transition, which takes place in the driven Jaynes-Cummings model with strong coupling between the two-level system and the harmonic oscillator. For a certain range of drive strength, the stationary solution corresponds to a bistability of classically distinguishable states. By unraveling the stationary solution into quantum trajectories, we resolve the nature of coexistence of phases. We show that the bistability solution develops into a first-order phase transition in the thermodynamic limit. We calculate the finite-size scaling exponent numerically. Even in the thermodynamic limit, the stability of phases originates from the discrete spectrum of a small quantum system. The nonclassicality of the photon statistics leads to the possibility of switching between the phases.


Introduction
Quantum phase transitions (QPTs) are abrupt, symmetry-breaking changes in the quantum state of a large quantum system at a critical value of an external control parameter. The QPT takes place at zero temperature and thus it generally refers to a change of the ground state [1]. The ground states on the two sides of the critical point have different symmetries and correspond to different macroscopic observables when a quantum measurement is performed. The concept of QPT can be straightforwardly extended to dissipative (open) quantum systems. Assume that the system is excited by an external coherent field, e.g., by laser irradiation. The coherent driving competing with the various damping processes sets the system in a dynamical equilibrium which is not the ground state, moreover, it may be far from a thermal state due to energy flowing through the system. Clearly, the steady state will depend on external parameters. Just like the ground state of closed systems, the stationary state of driven-dissipative quantum systems can also exhibit non-analytic behaviour as a function of some control parameter [2,3,4]. All this can happen at zero temperature. The critical point is marked by diverging quantum fluctuations and corresponds to a symmetry breaking change of the steady-state [5,6]. Note that the presence of dissipation in the open system amounts to information leakage to the environment. Therefore, if such a transition occurs in the quantum system, the different steady states under continuous measurement are associated with different measured, 'classical' outputs. In this sense the quantum states are 'amplified' to phases with different macroscopic observables even in the case of a small system with only few degrees of freedom.
Besides continuous phase transitions, theoretical and experimental works pointed out recently that first-order phase transitions can also have their counterpart in open quantum systems [7,8]. First-order phase transitions are benchmarked by the coexistence of phases. Instead of a critical point separating two phases, here there is a critical region in which bistability and hysteresis occur. How can the coexistence of 'macroscopically distinct' states be accommodated into quantum theory? The state of a driven-dissipative system is expressed by a density operator. The steady-state density operator, which obeys a linear algebraic equation, must be a continuous function of the external control parameter of the system. At the same time, the density operator can represent a bimodal distribution, a mixture of two macroscopically distinct components. There can be then a 'rift' separating the two 'phases', meanwhile the transition of the density matrix can take place smoothly from one phase to the other when tuning the control parameter. It is the weights in the mixture which vary continuously as a function of the control parameter: one of the components vanishes gradually with the simultaneous growing of the weight of the other. The complete transition can happen in a finite range of the control parameter, which is the bistability range.
Such a transition can be observed in connection with the photon-blockadebreakdown phenomenon. Consider a finite discrete quantum system, such as a two-level atom, very strongly coupled to a harmonic oscillator. Because of the strong coupling, Figure 1. (a) In a quantum system consisting of a finite-level system strongly coupled to a harmonic oscillator, the lower part of the spectrum is anharmonic, but there are always higher-lying harmonic subsets of the spectrum. (b) The Jaynes-Cummings spectrum is a prototype of this behaviour (this part of the figure is to scale). Panel (i): the anharmonic part of the spectrum. An excitation tuned to a 1-photon transition misses the second rung of the ladder by a significant amount. That tuned to a 2-and 3-photon transition, misses the first rung and the first and second rung, respectively. Panel (ii): closely harmonic part of the spectrum. An excitation tuned to the transition between |6, − and |7, − is close to resonance also with the transition between |7, − and |8, − (the detuning is invisible on the figure's scale).
the low-lying levels form an anharmonic spectrum [9] such that an external coherent driving nearly resonant with the bare oscillator's frequency is out of resonance with respect to all the possible transitions from the ground state. The transitions are thus blocked by off-resonance and the system is confined into the ground state. This effect has become known as 'photon blockade'. Above a certain intensity of the driving, however, the system can be excited via higher-order (multi-photon) transition processes into the higher-lying part of the spectrum. Since one of the constituents is a harmonic oscillator, in the high-lying part there are harmonic subsets of the spectrum (cf. fig. 1), ladders of equidistant levels which can host a coherent state with large amplitude and welldefined phase. Following a low-probability multi-photon transition into this range of the spectrum, the system gets stabilized into such a semiclassical trapping state by the competition of coherent driving and decay. Such an excitation, i.e., the breakdown of the photon blockade takes place in the form of a bistability. In a finite interval of the drive strength, the steady-state density operator of the system is the mixture of the 'dim' and 'bright' phases, i.e., the ground state and a highly excited coherent state of the oscillator.
There is a number of quantum systems where the blockade of excitations due to off-resonance is broken down for strong driving, quite analogously, in the form of a phase transition. This has been studied with regards to the creation of Rydberg-state atom clusters described by Ising-type spin models [10,11,12,13,14,15] and in arrays of nonlinear cavity systems [16] or in a single semiconductor cavity polariton mode [17,18,19,20]. Bistability has also been observed in the vacuum-superfluid transition in a dissipative Hubbard model [21]. The effect is thus quite generic. The photon-blockade breakdown is special, on one hand, because the degrees of freedom involved are kept at a minimum, i.e., only a single atom is coupled to a single bosonic mode. On the other hand, because here the thermodynamic limit is achieved in an infinitely strong coupling limit. The anharmonicity of the spectrum plays an essential part in our thermodynamic limit, in contrast to previously-reported thermodynamic (classical) limiting cases of quantum phase transitions in the Jaynes-Cummings and Rabi models [22,23,24,25]. The bistability associated with photon-blockade breakdown has already been pointed out [7,26,27] and has recently been experimentally demonstrated in the case of an artificial three-level atom [8].
This paper is devoted to complete and justify the analogy between the photonblockade-breakdown effect and a classical first-order phase transition by means of a finite-size scaling analysis. This is needed because the bimodal steady-state density operator, demonstrated in previous works, does not tell us about the stability of the phases. It allows only for identifying the two phases and their relative weights. Think of, for example, the fluorescence measurement of a driven V-type three-level atom, the system that was used to demonstrate quantum jumps for the first time [28]. The atom can be excited into a long-lived metastable state suppressing the fluorescence. On the other hand, when the atom spontaneously decays back to the ground state, an electricdipole-allowed transition gets driven by a laser, resulting in bright fluorescence. Dark and bright periods alternate according to the internal state of the atom. In the measurement process, the atomic quantum state is translated to a strong classical signal, however, these are not robust in the sense that single spontaneous quantum jumps lead to the macroscopic change of fluorescence intensity. The characteristic dwell times in the dark and bright periods of the blinking signal are determined by a microscopic timescale, namely, the lifetime of the metastable atomic state. In contrast, for a phase transition, there must be a 'thermodynamic' limit in which these components become stable -in the sense that no single quantum jump or other microscopic event can tip the system from one phase to the other -and the initial condition determines the final state.
The thermodynamic limit, following the suggestion of Carmichael [7], is taken as infinite strong coupling between the two-level atom and the harmonic mode, and use the coupling strength as the scaling parameter. Then we construct an appropriate finite-size scaling: we determine the scaling of the other parameters along which any given system with finite parameters can be 'scaled' up to the thermodynamic limit while keeping the system self-similar in a well-defined sense. The finite-size scaling reveals whether the stability of the phases increases, that is, whether their characteristic lifetime tends to infinity with increasing scaling parameter. If so, which we will show to be the case here, the bistability in the given realization of the system can be considered the finite-size approximation to what is a genuine phase transition in the thermodynamic limit.
The rest of the paper is structured as follows. In section 2 we recall the model and the basic ingredients of the photon-blockade-breakdown effect, together with the equations governing the underlying classical phase transition. In section 3 we present the time evolution of the system in steady state, which we show to be a telegraph signal. We undertake a detailed numerical analysis of timescales and the statistics of dwell times in the attractors. In section 4, the finite-size scaling of the parameters is determined, and section 5 is devoted to the analysis of the stability of phases. We shortly comment on the role of atomic decay in section 6 and then conclude in section 7.

The photon-blockade-breakdown phase transition in a nutshell
The simplest possible system exhibiting the photon-blockade-breakdown effect is composed of a two-level system coupled to a harmonic oscillator. The two-level system can be an atom or artificial atom, whereas the oscillator can represent a single lossy mode of the radiation field or a longitudinal mode of a stripline resonator. We describe this interaction within the driven Jaynes-Cummings model, i.e., using the electric-dipole coupling and the rotating-wave approximation: where ω M is the angular frequency of the mode with boson operator a, ω A is that of the atomic transition with step-down operator σ, and g is the single-photon Rabi frequency expressing the coupling strength. The external coherent drive is characterized by the amplitude η, here taken to be real for simplicity, and the frequency ω. Assuming resonance between the mode and the atom, i.e. ω M = ω A , and going into a frame rotating at the drive frequency, one gets a virtually time-independent Hamiltonian, where the detuning ∆ = ω − ω M is a tunable parameter of the drive. The mode is that of a high-finesse resonator and is subject to loss. Similarly, the two-level system can have decay through spontaneous emission. These incoherent processes can be modelled by Liouvillian terms in the master equatioṅ In the rest of the paper we will consider the case κ γ, most importantly, γ = 0. The mode relaxation parameter κ defines the microscopic timescale of the problem.
For weak drive strengths η g, the excited eigenstates of the Hamiltonian are close to the Jaynes-Cummings dressed states with n = 1, 2, . . ., and the energy levels are In the strong coupling regime g κ, the level shifts ± √ n g with respect to the bare frequencies exceed significantly the linewidth ∼ κ, so the system cannot be excited out of the ground state |g, 0 by a near-resonant driving ∆ ≈ 0. This is the photon blockade effect. ‡ Obviously, there is a strong-enough drive strength η such that the power broadening exceeds the detuning, whereupon the system can be excited into a state with large quantum number n. Recently, it has been argued that the transition between the two possibilities takes place in the form of a phase transition [7,26].
According to Carmichael, the underlying classical theory of this phase transition is provided by the so-called neoclassical equations, that are the mean-field approximation to the Heisenberg equations of motion in the special case of γ = 0. The mean-field equations are obtained by replacing the operators with c-numbers, e.g. a → α. Hence the means of operator products in the nonlinear equations are reduced to the products of means. The theory leads to the self-consistent equation (valid for ∆ < 0) where we introduced the parameter N scale = g 2 /4κ 2 . There can be different solution sets. The various domains are depicted in the phase diagram in fig. 2. Below the lower boundary, the phase is the photon blockade regime, whereas above the higher boundary, the system is highly excited. In between, eq. (6) has multivalued solution indicating bistability. The coordinates were chosen the tunable variables, i.e., the frequency and amplitude of the driving, which can serve as control parameters of the phase transition. Only the ∆ < 0 half-plane is shown, the positive-detuning part being the same mirrored to the ∆ = 0 axis. The neoclassical result suggests an appropriate thermodynamic limit [7] and a corresponding scaling of system parameters. On fixing the timescale to the microscopic one, κ = 1, we see that a characteristic photon number is expected to scale as ∼ g 2 . Hence, the thermodynamic (infinite-system) limit is the strong coupling limit g → ∞.
Simultaneously, the drive amplitude must also be scaled up. The first guess, cf. the prefactor on the right hand side of eq. (6), would be η → ∞ such that η/g is kept invariant. This is why the quantity η/g is used for the drive amplitude on the vertical axis of the phase diagram. With this scaling, the lower boundary of the bistability phase becomes indeed independent of g: the three curves coincide perfectly. On the other hand, the upper boundaries reveal a further dependence on g. That is, the solution sets of the neoclassical equations are not invariant under the transformation of g, η → ∞ with η/g = const. Later, in sec. 4 we will identify the non-trivial scaling of η which leads to a properly defined thermodynamic limit of the system.
Nevertheless, since the lower boundary of the neoclassical bistability domain is invariant, and also the upper boundary does not vary strongly at ∆ = −5κ for the g ‡ In the literature, the 'photon blockade' sometimes denotes the effect when the first excited state with a single photon can be excited resonantly, but further excitations are suppressed due to off-resonance. This is analogous to the effect of Coulomb blockade to some extent. Here we use the term in a more general sense, where the system is blocked in the ground state, and no photons at all can be transferred to the system. As the numerics becomes very unstable when approaching the far-off-resonance closing point of the bistability region, the closing of the red curves on the left is inferred by extrapolation (dashed segments). The common closing point on the right for all g values is the spontaneous dressed-state polarization critical point. The cyan star denotes the workpoint chosen in this paper: it is at this detuning and around this drive strength that we are going to study the bistable solution for several g values. Magenta dashed line: the inference (8) for the lower limit of the bistability region.
values shown in fig. 2 and used in this paper, in the plots the drive strength η is given in units of g.
As a side remark, let us note that there is a critical point at η/g = 1/2 specific to the resonant driving case ∆ = 0, where the lower and upper limits of the bistability region converge. It separates the solution α = 0 with population inversion increasing from −1/2 to 0 from the one with σ z = 0 and increasing α as the drive strength η is increased further. This result is in accordance with that of the full quantum treatment which can be pursued to an analytical solution in the resonant case [29]. It shows that the quasienergies coalesce in this critical point. This critical behaviour was identified as the spontaneous dressed-state polarization by Alsing and Carmichael [30].
Numerical simulations of the full quantum problem defined by eq. (3) confirm qualitatively the phase diagram based on the neoclassical equations. The existence of a bistability regime close to resonance |∆| g has been confirmed [7,26]. The stationary density-operator solution of the master equation is then a statistical mixture of two states: the 'dim' state where the field is close to the vacuum and the atom is in the ground state, and a 'bright' state which consists of a highly excited coherent state of the field. On increasing the drive strength η at a fixed detuning, the relative weight of the two components is continuously varied such that the probability of the bright component goes from zero to 1 in a finite range of η [8]. The steady state is hence a continuous function of the parameters, however, there is a 'rift' between the two components of the mixture: these are classically discernible states.
The stability of the bright component for high-enough drive strength can be understood intuitively as a balance of coherent drive and spontaneous decay in a harmonic oscillator. The frequency separating adjacent dressed states |n + 1, − and In the limit of large photon numbers n → N , this tends to a harmonic ladder in which the coherent drive η and the decay κ compete to create a coherent state with amplitude where the self-consistent equation for the photon number N was obtained by taking the absolute value squared of the amplitude α. The smallest drive strength for which this equation can be satisfied is in the case of 'resonance', i.e. when the expression in the parentheses in the denominator vanish: ∆ = −g/(2 √ N ). The self-consistent solution is then N = (η/κ) 2 , from which the minimum drive strength follows as η min g κ 2|∆| .
This law is drawn in the thick magenta line in the phase diagram in fig. 2, and coincides quite accurately with the lower border of the classical bistability domain. It is remarkable that the solution of a classical equation obeys a law extracted from intuitive consideration of quantized energy levels. On further increasing the drive strength, the photon number increases, thereby leading to some detuning ∆+g/(2 √ N ) = 0 in eq. (7). However, there may still be a self-consistent solution. Because the lower part of the harmonic ladder is missing (that part of the spectrum being the anharmonic photonblockading part), there is no deterministic evolution into such a self-consistent solution. Nevertheless, the displacement operator corresponding to coherent driving contains multi-photon excitation processes, so the excitation into the near-resonant regime with a self-consistent solution can take place in a probabilistic manner.
The upper limit of the bistability domain cannot be determined from an argument as simple as the above for the lower limit. The reason is that the more we increase the drive strength η, the more the quasi-energy levels, i.e. the true eigenvalues of the Hamiltonian (1), differ from the dressed levels of the η = 0 Jaynes-Cummings model, since they are getting dressed also by the coherent drive [29]. However, the analytical form of the quasi-energy levels in the finite-drive strength case is not known for ∆ = 0, only for the case of ∆ = 0. Nevertheless, it is clear that the appearance of an η-dependence of the energy levels makes that the η/g scaling suggested by both eqs. (6) and (7) is disrupted for large η values, which puts a finger on the g-dependence of the upper curves that we observe in fig. 2.

The telegraph signal
Not only the classical phase diagram, but even the steady-state density operator solution of the full quantum problem defined by eq. (3) does not describe all the relevant aspects of our phase transition. In the following we prove by numerical simulation that the components of the mixture become robust classical attractors in the thermodynamic limit. To this end, we need to unravel the density operator into the time domain, so that we can extract the dwell timescales, we can show their divergence, and we can determine the relevant exponents. To this end, we use the quantum trajectories generated by the Monte-Carlo wavefunction method. In principle, the ensemble average of many trajectories yields a density operator evolving in time towards the steady-state one. However, in the ergodic case (which will be our assumption here), the temporal averaging of the stochastic state vectors along a single long quantum trajectory yields the same steady-state density operator. Therefore it makes sense to consider a trajectory as an actual evolution under continuous measurement with an ideal photodetector.
Obviously, due to the large difference in the photon numbers, the components of the mixture correspond to very distinct output signals. The photon number is continuously monitored via the photons outcoupled into the κ loss channel. The classical distinguishability of large photocurrent versus dark counts amounts to a projection of the quantum state into only one of the components at a time. This is shown in fig. 3, where the instantaneous photon number along quantum trajectories is plotted for various coupling strengths g. The bright and dark periods alternate sharply in the form of a telegraph signal. Detailed analysis of the statistical data shows that the presented signals are indeed very accurately (even to the limit of numerical accuracy) described by a telegraph process. § This means that such trajectories have essentially three parameters: the amplitude of the bright period and the rates of blink on and off, µ and λ, respectively, since a telegraph process is nothing else than the composition of two temporal Poisson processes with exponential waiting-time distribution. Hence, the waiting time for a blinkon (the inverse of the blink-on rate µ) equals the dwell time in the dim period, the same being true for blink-off and the bright period. The characteristics of the telegraph process are summarized in table 1.
The trajectories for different coupling constants g in fig. 3 are generated using the drive strength η/g = 1/4 kept invariant, while fixing the detuning ∆ = −5κ. Hence all these curves correspond to the single point denoted with the cyan star in the phase diagram fig. 2 within the bistability domain. The photon number increases with increasing g. The numerics shows that the bright state has the photon number N ≈ g 2 /(2∆ 2 ), the corresponding dotted straight lines fitting nicely on the noisy § In particular, what we do is to define a binary signal from the somewhat noisy trajectories whose model is depicted in fig. 3, simply by assigning the value 1 to the time instants where the photon number is higher than half of the temporal average and 0 to the others. On this binary signal X(t) we verify the fulfillment of the relation var{X} = X (1 − X ), that is characteristic of the telegraph process (cf. table 1). We find agreement up to 10 −14 precision.   fig. 2). The dotted lines in panels a) and d) represent the estimate g 2 /(2∆ 2 ) for the photon number, following from eq. (7). Panel b) shows the phase of the field along the green trajectory of panel a) (g = 50κ), assuming a coherent state. The dotted green line is the field phase expected from eq. (7), that is, the phase of , substituting the bright-state photon number N = 50 that can be read from panel a). We see that the coincidence with the simulated field phase in the ON periods is remarkable. Panel c) displays the Mandel-Q parameter along the same trajectory, exhibiting a larger nonclassicality for the OFF periods than the ON periods, which we will discuss in section 5. Throughout this work, the numerics was performed with C++QED: a C++/Python framework for simulating open quantum dynamics [31,32,33].  numerical measurement record. This N is twice the photon number at the lower boundary of the bistability range of the neoclassical theory. Interestingly, this analytical estimate satisfies the simple classical relation (7) with high accuracy, which suggests that the intuitive picture of near-resonantly driving an equidistant ladder holds even at such drive strength far above the lower boundary of the bistability domain. The same is suggested by fig. 3(b), where the field phase is plotted together with an estimate from the same eq. (7), once more to yield excellent agreement (cf. the figure caption for details).
Since N scale ∝ g 2 , the photon number when measured in units of N scale proves to be invariant for the different curves. On the other hand, the dwell times in the attractor states increase significantly with increasing g. This reveals that there is a thermodynamic limit in which the phases become robust, the telegraph signal would disappear and it would be replaced by a hysteresis-like behaviour in a genuine first-order phase transition. Figure 3(d) shows that the addition of a small amount of atomic decay (γ = 0.01κ) does not lead to a qualitative change of the conclusions above. Note that the atomic decay noticeably decreases the dwell time in the ON period but not the photon number.

Filling factor and the scaling of the drive strength in the thermodynamic limit
In the bimodal steady-state density operator, the weights of the components vary through the bistability domain, in particular, on increasing the drive strength the 'dim' state component vanishes gradually in favour of the 'bright' state. In the time domain, the telegraph signal manifests this form of transition through the variation of the filling factor which is theoretically F = µ µ+λ . The trajectories in fig. 3 show that the filling factor varies for the telegraph signals with different g. This is confirmed in Fig. 4. Panel (a) shows the monotonous increase of F as a function of η/g for a set of g values. The curves are shifted with respect to each other. This dependence is made explicit in panel (b), where g is varied while keeping η/g at various fixed values. The filling factor is constant in the range of smaller g values, whereas there is a decrease of F in the range of larger coupling strengths. For example, the green line represents a closely constant filling factor around 2/3 up to g ≈ 50κ, but if g is increased further, the filling factor drops.
It follows then that the scaling of the drive strength η such that η/g is kept invariant does not preserve the self-similarity of the telegraph signal. Scaling to the thermodynamic limit means that the system keeps a self-similar behaviour, only scaled up in time and brightness of the ON state. Since the latter two scale with the single parameter g, there is the parameter η left at our disposal to ensure self-similarity during the scale-up, which, in the case of such a simple process as the telegraph one, cannot mean anything else than keeping the filling factor constant. As the most obvious choice, we pick the case when the filling factor is 0.5.
Therefore, in order to determine the correct scaling of η in the thermodynamic limit, we have to find the drive strength for each g value (denoted η * (g) in the following) that leads to half filling. Looking at fig. 4(a), we observe that an approximate linear interpolation (linear fit) captures appropriately the behaviour of the F(η) curves for the different g values in the range of interest. Hence, we use such an interpolation to Figure 5. The drive strength η * leading to half filling of the telegraph signal as a function of g, i.e., the equal mixture of the two phases in the steady-state density matrix. Solid red curves represent the phase boundaries of the bistability domain from the classical theory.
determine η * (g). As an example, in the plot we show for the red curve (g = 70κ) the corresponding η * . In fig. 5, we plot the function η * (g), embedded in the bistability domain that we have calculated similarly to fig. 2. As it turns out, a linear fit reproduces quite exactly the behaviour of η * for g 50κ, which indicates that the correct scaling of the drive strength in the thermodynamic limit is η ∼ g 2 . Figure 6 presents the characteristic timescale of the bistable blinking process. This is defined as the inverse of the sum of blinking rates µ + λ, which is extracted from an exponential fit on the numerically calculated temporal self-correlation of the signal (cf. table 1). We use the units of κ so that the large difference is manifested in the figure: the characteristic times are orders of magnitude above the microscopic timescale κ −1 . On increasing the drive strength η in the presented range, the average dwell time decreases because of the increase of the rate µ of jumping upward, in conjunction with the increase Figure 6. The characteristic timescale of the bistable blinking process given by the inverse of µ + λ. The points in panel a) correspond to the values of η * for the different g values where half-filling of the telegraph signal is achieved. In panel b), the solid red curve designates the timescale at η * as a function of g, that is, the timescale in the correct finite-size scaling when the telegraph process is kept at half filling during the passage to the thermodynamic limit. The dotted line shows a log-log linear fit on the values leaving out g = 100. The exponent resulting from the fit is roughly 2.2. Panels c) and d): timescale with finite atomic decay γ = 0.01κ. of the filling factor. Points represent the values of η * where the telegraph signal has halffilling.

Dwell times
From the point of view of the thermodynamic limit, the dependence of the characteristic timescale on g is the most relevant ( fig. 6(b)). We show the increasing timescale over two orders of magnitude of the coupling constant g for various values of η/g. This g range is given by computational limitations, nevertheless, it is enough to demonstrate the power-law scaling of the increase of the timescale and to determine the exponent. To obey the correct self-similar scaling detailed in section 4, we have to find the timescale for the drive strength η * (g) for the different g values. To this end, we again use a linear fit on the curves of η-dependence, which captures the behavior quite correctly (cf. dotted lines in the panel (a), with the η * values depicted with big dots of the corresponding colour). The timescale change under finite-size scaling τ (g, η * (g)) is shown in the thick, solid, red line in fig. 6(b). A linear fit in the log-log scale leads to the numerical estimate of 2.2 for the finite-size scaling exponent of the characteristic time. The point g = 100κ was omitted from the fitting because the dwell times were systematically underestimated due to truncation of the trajectories for the long but finite simulation time. Figure 7 is devoted to the numerical analysis of the rates of jump upward (blinkon) and downward (blink-off), µ and λ, respectively, which sheds light on the physical processes leading to switching between the robust classical attractor states in the finitesize system. These are calculated from combining the above-discussed characteristic timescale with the filling factor. For both rates, the increase of g implies a reduction, in agreement with the expectedly growing stability of phases on approaching the thermodynamic limit.
The downward process is initiated by a single photon loss (detection) event. This is because there is a chance that the state residing in the ladder |n, − gets projected into the ladder |n, + under such an event [7], since Once such a jump occurs, there is a downward cascade of photon escapes, because while our (red-detuned) drive was closely resonant with the high-lying part of the |n, − ladder, resulting in an approximately coherent state in this part, it is off-resonant with the |n, + ladder (which would be resonant with the drive blue-detuned with the same amount), so on this ladder there is no drive to compensate the photon loss. The passage If we kept η constant when performing the limit g → ∞ (cf. the bleached lines in fig. 6(b)), the exponent would be about 2.5, with very little dependence on the value of η. downward consists of a quick cascade of many jumps amounting to an exponential decay.
Although a single ladder-switching quantum jump initiates the blink-off, however, the likelihood of such a rare jump vanishes completely in the g, N → ∞ thermodynamic limit. The rate of ladder switching scales as κ/N , ¶ and as we saw above, the brightstate photon number scales as g 2 in the thermodynamic limit. This would suggest the exponent 2 for the finite-size scaling of the blink-off timescale, which is verified in the inset of the right panel of fig. 7 to very good accuracy. This downward jump rate λ is largely independent of η, that only slightly influences the photon number in the bright phase.
On the other hand, the switching from the dim to the bright phase is induced by the external driving and thus is sensitive to η, as can be seen in the left panel of fig. 7. The upward process is suppressed by the off-resonance of low-lying quasi-energy levels (anharmonic part of the spectrum). Therefore, it is easy to understand that the larger the coupling g, the larger the shift of the levels from resonance and the smaller the jump rate. In the dim phase, the state of the bosonic mode is close to the vacuum, however, there must be a small deviation from that due to the driving. The state is a non-classical superposition with positive Mandel-Q parameter. The Mandel-Q parameter measuring the nonclassicality of the field state is defined as where the averaging can be performed either as a quantum average on the actual state vector of the field to reflect the nonclassicality at a given time instant, or also over time. The Mandel-Q parameter is zero for a classical field state (coherent state). In fig. 3(c), Q is calculated as a time-dependent quantum average, and we observe that the nonclassicality is stronger in the dim phase than in the bright one. This is consistent with our picture that the bright phase consists of an approximately coherent state on the manifold |n, − . The state in the dim period has the property that the projection of the wavefunction after a photon detection increases the weight of a high-excitation component. While in the dim state there is a negligible amount of excited photon component generated by the η driving, triggered by a single quantum jump (which is very rare on account of the very low dim-state photon number), it can grow in an exponential runaway process for subsequent photon detections. So the blink-on also takes place in the form of a cascade of quantum jumps.
In fig. 8, we plot the field Mandel-Q parameter averaged over time, but only in the dim periods. The overall trend is that the nonclassicality of the field in the dim phase increases both with increasing η and g, and the dependence flattens out for large g at a value close to 1. Hence, the dim phase remains nonclassical also in the thermodynamic limit. Let us try to model the dim state of the field in the following form: where ε is a small number and |ϕ is a state orthogonal to the vacuum state. The photon count rate κ a † a = κ ε 2 N |ϕ can be made very small with ε → 0, where N ϕ = ϕ|a † a|ϕ is the mean photon number associated with the component |ϕ superposed on the vacuum. At the same time, the Mandel-Q parameter of the state (11) is found to be independent of ε in the limit ε 1: one finds Q |Ψ = Q |ϕ +N |ϕ . That is, it depends only on the properties of the state |ϕ . Note that (i) N |ϕ > 1 because the vacuum component is missing from photon number expansion of |ϕ , and (ii) the Mandel-Q parameter is always limited below 1 in fig. 8. These two observations imply that the state |ϕ is a nonclassical state with negative Mandel-Q factor.

The role of atomic decay
The condition γ = 0 is essential in the neoclassical theory, since the derivation of the transcendental eq. (6) relies heavily on the fact that the length of atomic pseudo-spin is conserved, σ x 2 + σ y 2 + σ z 2 = 3/4. Allowing γ = 0 hence leads to qualitatively different behaviour since this conservation law is broken.
Yet, for very small γ values, the behavior of the full quantum model does not seem to be qualitatively affected: as we have shown in the course of this study in passing, the photon-blockade-breakdown effect can be observed in the case of a finite but small γ.
The effect of an atomic decay in the bright state, i.e., on a high-lying coherent state can be assessed similarly to the above: This means that in the event of an atomic decay, there is a 1/2 probability of a ladder switch. Hence, a γ on the order of κ would wipe out the blinking effect that is the central theme of this paper, since a blink-on would immediately be followed by a collapse of the bright state. The blinking effect manifests itself most clearly when γ = 0, the case that we considered mostly here. If γ κ, which case we exposed in passing in figs. 3, 4 and 6, the atomic decay emerges as a competing timescale for high-enough bright-state photon numbers, resulting overall in smaller characteristic times and filling factors. With respect to the phase transition, the system is interesting as long as the microscopic timescale of spontaneous emission is longer and thus is dominated by the shorter macroscopic timescale τ of the phase stability.

Conclusions
We showed that the bistability observed numerically in the driven Jaynes-Cummings model with very large coupling g in a certain domain of the drive parameters is the finite-size precursor of a first-order phase transition in the thermodynamic limit. Let us emphasize that the thermodynamic limit does not lead to a classical system with many degrees of freedom. In contrast, it can be defined by properly redefining the values of the parameters in the same model of a coupled two-level system and a boson mode. The thermodynamic limit of large coupling constant g amounts to enhancing the anharmonicity of the discrete spectrum of the quantum system. This is the mechanism needed for the stabilization of the 'dim' phase and is thus responsible for the bistability.