Rare quantum metastable states in the strongly dispersive Jaynes-Cummings oscillator

We present evidence of metastable rare quantum-fluctuation switching for the driven dissipative Jaynes-Cummings oscillator coupled to a zero-temperature bath in the strongly dispersive regime. We show that single-atom complex amplitude bistability is accompanied by the appearance of a low-amplitude long-lived transient state, hereinafter called `dark state', having a distribution with quasi-Poissonian statistics both for the coupled qubit and cavity mode. We find that the dark state is linked to a spontaneous flipping of the qubit state, detuning the cavity to a low-photon response. The appearance of the dark state is correlated with the participation of the two metastable states in the dispersive bistability, as evidenced by the solution of the Master Equation and single quantum trajectories.


I. INTRODUCTION
Fluctuation-induced switching between metastable states of driven quantum nonlinear oscillators interacting with their environment, in principle lacking detailed balance, constitutes one of the most general and intricate physics problems, intimately linked with the problem of quantum activation [1][2][3][4]. In this framework, the Duffing model provides the simplest description of a self-interacting nonlinear oscillator involving one quantum degree of freedom, where the Fokker-Planck equation (FPE) is sufficient to provide an exact treatment of quantum fluctuations. Such a formulation ultimately yields a steady state which presents notable differences from a Gaussian distribution, as one expects in linear FPEs. The FPE can be solved exactly in the steady state, since the conditions for detailed balance are satisfied at zero temperature [5]. It has recently been shown that a coherently driven system with two quantum degrees of freedom, i.e., a transmon qubit coupled to a resonant cavity mode, both connected to a dissipative environment, may still be amenable to a FPE description subject to an adiabatic elimination of the fast decaying cavity field amplitude [8].
Historically, the problem of defining switching rates between states of classical nonlinear dissipative systems is long standing. The driven Van der Pol oscillator is a very characteristic case subject to a description where an effective potential V (x) can be devised as a function of the driving parameters for the nonlinear drift term * t.mavrogordatos@ucl.ac.uk of the FPE with constant diffusion. We can then define the forward (associated with an energy gap Q + ) and backward (associated with the gap Q − ) jump rates for the phase fluctuations, with a Kramers-type dependence: When we analyze the single-atom Jaynes-Cummings (JC) model, we face a complexity which transcends the difficulty of solving a nonlinear FPE, as in the aforementioned oscillators. In particular, it is impossible to define an FPE to study the switching dynamics, which is a consequence of the non-perturbative nature of light-matter interaction in the strong-coupling regime [12,13]. In contrast to the dispersive optical bistability, discussed in [14] as a special case of a dissipative system with a potential [15], for the case of single-atom JC bistability we cannot formulate a suitable potential function yielding the various attractors in the phase portrait. The inability to obtain a potential force may lead to large deviations from the optimal path minimizing the action for two degrees of freedom [1,16].
In this paper, we report on a metastable state, called the dark state, with very low intracavity amplitude and intense qubit fluctuations which is not predicted by the mean-field equations. We identify this state as resulting from quantum bistability involving two degrees of freedom [18], and we depict the state on a plot of the associated quasi-distribution in the coherent phase space for the cavity field, as well as the associated distribution in the Bloch sphere (see Figs. 1 and 12). We find that the dark state is (1) rare, appearing only on the longer time scales after the transient period, (2) very noticeable and strongly fluctuating in the Bloch sphere with regards to the qubit observables, (3) long lived, compared to the typical time scales of cavity and qubit dissipation, and (4) fragile to fluctuations yet more 'resilient' than the unstable state of mean-field dispersive bistability.
In Sec. II we define the system Hamiltonian and main methods used for analyzing quantum bistability. In Sec. III we first approach with an approximate mapping to the Duffing oscillator for small drive strengths, by examining a perturbative expansion of the driven dispersive JC Hamiltonian, with renormalized parameters to account for the cavity-atom coupling. As the driving power is further increased, the Duffing approximation breaks down because both the qubit and cavity participate in the bistable switching. Subsequently, in Sec. IV we study the switching dynamics in single quantum trajectories, linking our discussion to the mean-field and neoclassical predictions. While the quantum-fluctuation switching takes place in the steady state, the qubit flips and brings the cavity mode out of resonance, with a very low photon excitation.

II. MODEL AND METHODS
We will first provide a brief account of the properties of dispersive complex amplitude bistability for varying drive strength and frequency. In a frame rotating with the drive frequency ω d , the Hamiltonian describing the interaction of a damped two-level atom (qubit), with inversion operator σ z and raising (lowering) operators σ + (σ − ) with bare resonant frequency ω q , and a driven damped cavity mode (with photon annihilation and creation operators a and a † respectively) with bare frequency ω c , reads [19]: where ∆ω c,q = ω d − ω c,q , g is the atom-cavity coupling strength and ε d is the drive amplitude (or strength). The cavity is coupled to a thermal bath at zero temperature inducing a photon loss rate of 2κ, while the qubit relaxation rate is denoted by γ (due to both radiative and nonradiative processes, such as quasi-particle formation).
The system density matrix obeys the Lindblad master equation (ME) [19]: which is solved numerically via exact diagonalization, as well as unraveled into single quantum trajectories. The steady-state solution of (2) ρ ss is used to calculate the average value of the observables considered, as O ss = tr(ρ ss O), where O is a system operator. A normalized conditional state unraveling the ME is subject to an evolution obeying the Stochastic Schrödinger Equation (SSE): where D 1 is the drift term, D 2 is the diffusion term (both functions of Lindblad operators), and dW is a real increment (for more details, see [20,21]). The density matrix ρ k (t) = |ψ k (t) ψ k (t)| is used to calculate the average value of system observables as O(t) = tr(ρ k (t)O) [13,20]. In all of the cases considered here, the initial state ψ k (t = 0) is a pure state with ψ k (0)|a † a|ψ k (0) = 0 and ψ k (0)|σ z |ψ k (0) = −1, unless explicitly stated otherwise. The properties of steady-state bistability were not affected by a change in the initial conditions for single quantum trajectories. Convergence with respect to the time step in the evolution as well as in the number of states comprising the truncated Hilbert space for the cavity has been ensured. The strongly dispersive regime is defined by an atomcavity detuning δ ≡ |ω c − ω q | of the order of (and usually greater than) the coupling strength g 2κ, γ, alongside its relation to the coherent drive strength: max(2κ, γ) < ε d g 2 /δ, which takes us beyond the linear dispersive regime. For the detuning that we are considering here, ∆ω q > g and 2κ < ∆ω c < g 2 /δ. For all cases discussed in this work, δ ω c + ω q while ω c , ω q g, so that the rotating wave approximation (RWA) can be safely performed (see also Fig. 1 of [22]) for the number of photons involved in the steady-state response. For the average photon number n ss in the steady state, we would typically have g n ss + 1 ≤ 0.1 ω c . The standard ME can then adequately describe many cavity QED and circuit QED experiments [23], both at resonance and in the dispersive regime (for a direct comparison between theory and experiment see e.g. [24] and [18]).
In our treatment we have not included the phase destroying term (γ φ /2)(σ z ρσ z − ρ) in the ME [12]. Such a term produces the decay coefficient γ + 2γ φ for the qubit coherence, erasing even more rapidly the memory of the initial state in the averaged system response. Qubit dephasing would affect the lifetime and fluctuations of the states of quantum bistability together with the scaling constants of the mean-field equations, both present when the (energy) decay coefficient γ is already taken into account. For the limiting cases considered later on (see Fig.  11), we would take γ → 0 together with γ φ → 0.

III. FROM THE EFFECTIVE DUFFING OSCILLATOR TO THE FULL JC NONLINEARITY
A. The effective Duffing oscillator The quantum Duffing oscillator is a precursor of the JC nonlinearity. After applying the dispersive transformation [25,26], the Hamiltonian of Eq. (1) in the dressedcavity Duffing approximation for δ g reads (up to quartic order in the small parameter g/δ): In the above expression we have kept only linear terms with respect to g/δ in the transformed drive term, while we have set σ ± = σ ± = 0, taking σ z = σ z ≈ −1 (see Eq. 3.15 of [25]). Based on Eq. (4), we can extract the Wigner function for the effective Duffing oscillator [27,28], calculated via the generalized P -representation (see the Appendix for a full derivation): where 0 F 1 (a; x) and 0 F 2 (a, b; x) are generalized hypergeometric functions of the variable x with parameters a, b.
Here, c = (κ − i∆ω c )/(iχ) with the renormalized detuning ∆ω c = ∆ω c + (g 2 /δ)σ z − (g 4 /δ 3 )(2σ z + 1) and The Wigner distribution function of Eq. (5) is valid for 4 N ss g 2 /δ 2 1 (here the subscript ss denotes the steady state, and N = a † a+σ + σ − is the number operator of system excitations; see [25] and [29] for more details) and can be used for the calculation of the intracavity field moments in the complex plane as opposed to the four-dimensional space of [6]. Valuable information can be extracted from the effective Duffing oscillator model for low drive strengths, where σ z ≈ σ z ≈ −1. The expression of Eq. (5) for the Wigner function of the renormalized Duffing oscillator predicts a variety of critical points surrounding the dim state (in agreement with the low-amplitude bistability plots presented in [28]). The steady-solution of the ME yields an almost identical distribution, capturing the same amount of nodes in very similar positions. In the regime of low intracavity amplitude, quantum fluctuations are essential for the onset of complex amplitude bistability, determined by the scale parameter δ 2 /(4g 2 ). With increasing drive strength, where the number of system excitations approaches the scale parameter, the Duffing approximation becomes inapplicable, and the full JC dynamics with two quantum degrees of freedom must be taken into account. For the driven JC oscillator, the semiclassical bistability region, characterized by one unstable and two metastable states (a dim state with lower photon occupation and a bright state with a higher photon number), is constructed in the drive parameter space (∆ω c /κ, ε d /κ) from the Maxwell-Bloch equations [19]. The latter are known to yield solutions that exhibit overlap between different domains of attraction and chaotic behavior [30,31]. An alternative construction can be carried out from Hamilton's equations of motion for time scales during which the qubit degrees of freedom can be considered as constants of motion for γ/(2κ) → 0 [29]. Taking now into account the quantum fluctuations, Fig.  1 shows the onset of complex amplitude bistability extracted from the ME solution for the steady-state intracavity field for a constant drive detuning and varying strength. The two metastable states, the dim (D) and the bright (B), exchange probability as the drive strength increases, being connected via the dispersive excitation spiral along which we can also discern the unstable state (U). Together with the increase of the intracavity photons, we can also observe the progressive separation between the dim and the dark states, which will also be demonstrated when the ME is unraveled into single trajectories. A good separation between the dim and the dark states is shown in Fig. 1(d). In the following, we will also discuss the behavior of the qubit observables when the perturbative approach leading to the Hamiltonian of Eq. (4) cannot be applied.
B. Mean-field and quantum trajectories away from the critical point C1 In order to understand the origin of the dark state we have at first ignored the cavity-qubit quantum correlations. While assessing single-atom dispersive bistability, the authors in [29] present a construction in which they depict the semiclassical bistability region as a 'leaf' in the phase space when γ → 0, opening at one critical point (C 1 ) in agreement with the effective Duffing oscillator, and closing at another (C 2 ) when the drive is at resonance with the significantly excited cavity. Quantum fluctuations of both the qubit and cavity field alter the overall shape of the leaf, bringing about non-equilibrium dynamics where quantum noise cannot be treated perturbatively. As the drive strength and the intracavity photon number are further increased, the quasi-distribution function of Eq. (5) fails to adequately describe the quantum dynamics, which now involve the qubit more actively. The region of coexisting states with probabilities of the same order of magnitude marks the boundary of the region where quantum fluctuations are important. Solving the ME for the system density matrix in the steady state, and tracing out the qubit degrees of freedom, we can investigate how the cavity bistability builds up with varying drive frequency and power. The meanfield predictions and the full quantum treatment are in closer agreement outside the bistability region. Inside this region, conversely, we expect a first-order quantum phase transition boundary marked by coexistent semicoherent states, clearly indicated by the Q function plots we present in Fig. 2. One metastable state has a low photon mean a † a , called 'dim' while the other one has a higher cavity excitation, called 'bright'. In Panels 1 and 2 of Fig. 2 we are traversing the steady-state quantum bistability region by varying the drive frequency at constant drive strength. The dim coherent state gives its place to the bright one while crossing the first-order transition line. In Panel 3 we present evidence of the dark state in the averaged response, where we are able to discern a center and a saddle point. The state appears to be adjacent and linked to the dim metastable state. In Panel 4 we attempt a pictorial analogy to the case where the photon loss rate 2κ is of the same order of magnitude as the spontaneous emission rate γ [going further away from the zero system size, i.e. the limit γ 2 /(8g 2 ) = 0], revealing that the bright and dim states are joined in probability transfer as a consequence of increased spontaneous emission. Furthermore, the variation of γ has an important effect on the steady-state distribution, resulting in the persistence of the dim state for the same drive strength (compare frames 1 and 4), as opposed to the low-power bistability. In that respect, a high photon number limiting behavior for this system far from equilibrium can be defined through the intracavity amplitude n scale = [δ/(2g)] 2 for which the nonlinearity in the response can no longer be treated perturbatively [32]. In the strongly dispersive regime, the presence of the small term g/δ precludes the divergence of nonlinearity at low intracavity amplitudes in the a priori absence of spontaneous emission and dephasing (γ, γ φ =0), as deduced from [32]. We will now focus on the switching behavior of a qubit coupled to a driven resonant cavity mode in the regime of dispersive bistability far away from the critical point C 1 in the semiclassical bistability leaf. Qubit switching is revealed by single quantum trajectories after calculating the reduced qubit density matrix, ρ k Q (t) = tr c {|ψ k (t) ψ k (t)|}. The dynamical organization of the quasi-coherent distributions is depicted in the Bloch sphere, which serves as an equivalent representation of the phase space for the complex cavity amplitude. Figure 3 shows the build-up of bistability in the qubit amplitude projected on the equatorial plane for constant driving power and varying driving frequency. Each point in the scatter plots corresponds to one time instant within a single quantum trajectory, generated by solving numerically an SSE under the diffusive approximation using an explicit weak scheme [20,21].
This representation is analogous to the development of cavity bimodality depicted in the Wigner function plots presented in Fig. 1. The dim state exhibits a concentration around the south pole of the Bloch sphere while the bright state approaches the equatorial plane with decreasing drive-cavity detuning ∆ω c . As ∆ω c increases, the bright state distribution moves towards the south pole, consistent with the approach of the Lorentzian lineshape outside the bistability leaf, where only one state (and consequently one distribution) is expected. An example of the mean-field distorted Lorentzian profiles within the bistability region are given in Fig. 4.
In Fig. 4 we show the mean-field dynamics for γ = 0 alongside the nonlinear cavity response in the steady state (see Ch. 11 of [19] for the relevant equations). Here we can observe the decreasing nonlinearity for decreasing coupling strength, marking the transition from a skewed Lorentzian curve to a linear response function. A solid line with an arrowhead intersects the blue curve at three points: the dim, the unstable and the bright states from bottom to top respectively. Mean-field (Maxwell-Bloch) bistability is present in the transient evolution as well, evidenced by the two distinguished modes in the Inset (c), exhibiting anew the limit cycle approach we encounter for the qubit in one single quantum trajectory. The Insets (a) and (c) focus on the approach of the dim state (marked with a point), with varying coupling strength.
Focusing now on a quantum trajectory, in Fig. 5 we present a field quasi -distribution function for the reduced cavity density matrix ρ k C (t) = tr Q {|ψ k (t) ψ k (t)|} at two time instants t 1 , t 2 during a period of macroscopic switching of the coupled cavity-qubit system to the bright state. We find that the three semiclassical states coexist along a spiral during a switch 'up' to the bright state. This figure, moreover, shows that the dim and the unstable (semiclassical) state are connected by two probabilityflow paths, similar to the development of bistability we have seen in Fig. A.1, and the dark-dim state connection we have depicted in Fig. 2 (c). A faint peak is just about visible in the third quarter of the phase space, establishing the radius of the spiral span due to the external drive. It is remarkable that all this information can be extracted from the instantaneous quasi -distribution function for the intracavity field alone, after the qubit degrees of freedom have been traced out.
The cavity nonlinearity manifests in a nonperturbative fashion only in the region of large photon numbers in comparison to n scale , far within the semiclassical bistability region. In Fig. 6 we show the variation of the moduli of the complex cavity amplitude and the qubit projection | σ − | as a function of the normalized The solid curves depict Lorentzian and skewed-Lorentzian profiles for the same drive amplitude, ε d /γ = 100: blue, skewed Lorentzian, with peak at level (2) further from resonance (∆ωc = 0) for g/δ = 0.14 and 2κ/γ = 12; red, skewed Lorentzian, with peak at level (2) closer to resonance for g/δ = 0.87 and 2κ/γ = 12; green, skewed Lorentzian, with peak at level (1) for g/δ = 0.14 and 2κ/γ = 22; black, Lorentzian with peak at level (1) for g/δ = 0.042 and 2κ/γ = 22; orange, Lorentzian with peak at level (2) for g/γ = 3347, δ/g = 0 and 2κ/γ = 12. The ratio ε d /κ, giving the empty cavity amplitude, determines the plateau for the Lorentzian peaks marked by (1) and (2) respectively. The peak of the orange curve indicates the bare cavity frequency. The three insets (a, b, c) represent solutions of the Maxwell-Bloch equations of the phase space x − y of the intracavity field for varying time, corresponding to the drive frequency ∆ωc/κ = 56.83. In (a) g/γ = 3347 and in (b) g/γ = 1000. In (c) we are plotting the intracavity field for g/γ = 3347 but for a longer time, showing the approach of the second semiclassical state (bright) in a limit cycle fashion. The dashed lines indicate the unstable branches.
drive phase space. The deformation of the Lorentzian shape is accompanied by a line in the phase space where the complex amplitudes of the two metastable states cancel coherently, as we can see in Fig. 6. In that region of pronounced quantum fluctuations, the correlation function g (2) ss (τ = 0) attains its maximum, a behavior which is also a discerning feature of the Duffing oscillator [6].
A limiting behavior, in the sense discussed in [32], is achieved for g → 0 (sending n scale to infinity). This is a weak-coupling limit [with the co-operativity parameter C = g 2 /(κγ) remaining constant] for which the fluctua-tions vanish and the ME results are in close agreement with the mean-field predictions. Nonlinearity manifests itself markedly differently, however, in the case of drivecavity resonance, where ω c = ω d : along the line ∆ω c = 0 the semiclassical amplitude bistability region closes up in contrast to the Duffing oscillator phase diagram [29]. The emerging phase bistability for δ = 0 is associated with spontaneous symmetry breaking alongside a second order phase transition, more evidence of the JC nonlinearity [13,32]. Amplitude bistability is again recovered in the presence of spontaneous emission for γ 2κ, with a sufficiently large saturation parameter γ 2 /(8g 2 ), The dispersive cavity excitation spiral. Quasidistribution function Q(x + iy) for two time instants t1 (a) and t2 (b) with t2 > t1 during a switch to the bright metastable state. Parameters: ε d /κ = 16.67, g/δ = 0.14, γ/(2κ) = 1/12, g/γ = 3347. as shown in [33]. In our case, however, this number is vanishingly small, and n scale sets the dominant scale for the manifestation of dispersive nonlinearity. The meanfield states of dispersive bistability are sensitive to both the initial conditions and the value of γ. Different phasespace portraits for the intracavity field are depicted in Fig. 7 for varying γ and initial conditions. In Fig. 7(a), the initial conditions have been selected to be very close to the dim steady-state amplitudes. For γ = 0 and γ = 0 we see an approach to the dim-state fixed point, which is however lost if the initial conditions change. A new state is approached in Fig. 7 (c) for γ = 0, different from the metastable states of dispersive bistability, which are recovered for γ = 0. We are guided by the mean-field results to explore the attributes of the dark state for three different drive strengths using single quantum trajectories and the (averaged) exact ME results, depicted in Fig. 8 (Panels I-III). The frames (c) and (d) attest that the dark state is characterized by intense qubit fluctuations which follow quasi-Poissonian statistics, to which the frame (d) of Panel II testifies. At the same time, there is clear evidence of the dark state from the ME steadystate cavity distribution in frame (b) of Panel II showing a particular excitation path linking the dark to the bright state, which is the only one anticipated by the Maxwell-Bloch equations. The dark state coexists with the dim state in frames A and B of Panel III in Fig. 8 until the dim state vanishes completely into the excitation probability path for increased drive strength. A quick look at the Bloch sphere of Fig. 9 suffices to convince us of the departure from the mean-field predictions. According to the Maxwell-Bloch equations, we would expect to find the qubit vector lying solely on the southern hemisphere FIG. 8. Drive parameter phase diagram and the dispersive JC bistability. Bistability leaf produced from Hamilton's equations, assuming σz is a constant of motion (see Fig. 2(a) of [29]). The two critical points C1,2 are marked together with three vertical cuts (I-III) following the transition from C1 to C2. Each of the three cuts corresponds to a panel given underneath for a different value of the cavity detuning. Panel I: of the Bloch sphere. Interestingly, these fluctuations are described by quasi-Poissonian statistics as well, with a mean inversion in the northern hemisphere.
To conclude this section, we will present the dim state and the two nodes identified as the dark state away from C 1 , in connection to Fig. 1. Figure 10 shows the Wigner function in a drive region where the qubit participates significantly in the dynamics, exhibiting large fluctuations in the Bloch sphere when transitioning between the dim and the dark states.
In the following section we will investigate the meanfield dispersive bistability in the absence of spontaneous emission, prompted by the behavior we have encountered in Figs. A.1 and 10. A change is heralded by a new scaling parameter relevant for the development of nonlinearity, namely δ 2 /(4g 2 ), as the equation for dispersive bistability shows:

IV. SWITCHING DYNAMICS IN SINGLE QUANTUM TRAJECTORIES
Let us now seek some evidence of the dark state within the bistable switching itself, as a result of the quantum fluctuations. We have performed the ME unrav-  Fig. 1(d). The dim state is denoted by D, the unstable node by d1 and the center by d2. Parameters: ∆ωc/κ = 9.167, g/γ = 600, 2κ/γ = 12.
eling through numerically solving Stochastic Schrödinger Equations (SSEs) using the second-order weak scheme in the diffusive approximation, as devised by Platen (see Ch. 15 of [20]). The presence of the dark state is associated with an intense fluctuation having a spectral content on the left of the drive frequency, far beyond the spectral peaks of the bright and the dim states. The spectrum of the coherent fields a(t) and σ − (t) is expected to be asymmetric with respect to the drive, because of the presence of dissipation. The dark state reinforces this asymmetry. When γ = 0 the lifetime of the dark state is significantly prolonged and comparable to that of the metastable states [ Fig. 11 (a)]. As the spectra of frames (b,c) in Fig. 11 evidence, when the quantum fluctuation switching involves the dark state, there is a peak in the spectrum located at −g 2 /δ in the rotating frame (excluding the time evolution independent of cavity-qubit coupling), corresponding to the qubit flipping from a state with σ z = −1 to a state with σ z = +1. At the same time, the appearance of the bright state is a key element to the switching. In Fig. 11(d) we encounter a situation where only the dark and the dim states are present, in which case there are no peaks on the left of the drive tone. The qubit flipping brings the cavity mode out of resonance, far detuned from the drive frequency. In the low-amplitude dispersive regime, the cavity response is a Lorentzian centered at ∆ω c = g 2 /δ. Hence, for ∆ω c > 0 and higher drive strengths, the appearance of the dark state can be construed as a spontaneous projection to a state with σ z ≈ +1 and very low intracavity excitation, which is not an expected mean-field solution for γ = 0. As the dark state is visited (with one center and one unstable node, as shown in Figs. 10 and 16), the qubit inversion transitions from σ z ≈ +1 to σ z ≈ −1, which explains the observed intense fluctuations seen in the quantum trajectories. The significance of the low-power regime is further ascertained by the analytical expression for the Wigner function we encountered for the Duffing oscillator with one 'active' quantum degree of freedom, which captures a variety of nodes apart from the Maxwell-Bloch states (see Panel I], whereas with increasing γ the dim state dominates and the dark state appears short-lived after the bright metastable state, as depicted in Panel II of Fig.  11. The histograms of that panel provide information on the lifetime of the dark state, which changes from about 20 to about 150 cavity lifetimes, on average, with diminishing spontaneous emission rate. We will now focus on the phase-space representation of the dark state during a transition involving a metastable state of the Maxwell-Bloch bistability. Regarding the salient features of the cavity amplitude quasi-distribution, we are already familiar with the spiral rotation in the phase space following the de-excitation path in the JC ladder at resonance in the presence of dissipation [see Fig. 3(b) of [32] for resonance]. We are also acquainted with squeezing in the quadrature along the mean-field direction, from resonance fluorescence [12,34] as well as from the Duffing oscillator [13]. Switching among metastable states means another swirl in the spiral established by intracavity (g) and intercavity (ε d , κ, γ) coupling, combining features of resonance fluorescence and decaying optical oscillations. On the one hand, such a representation reveals the statistical mixture of semi-coherent states involved in the switching itself, and on the other hand it provides details on the excitation path followed in the JC ladder. At a particular time instance during the decay of the unstable state to the bright state, probability accrues at the highly-excited cavity state while the bottom part of the spiral becomes more pronounced (in accordance with the fully-averaged results of Fig. 1). At that time we expect to find the qubit amplitude following a trajectory that en- circles the bright state (with σ z closer to zero than in the region of the low-excitation critical point) in a limit cycle fashion, as shown in Figs. 12 and 9. While this spiral is described, the dark state is occasionally visited, when the qubit vector is found on the north pole of the Bloch sphere. In Figure 12 we show joint bistability for the cavity and qubit, when the bright and dim state distributions are significantly separated. The entanglement entropy attains its highest values during the occupation of the dark state [see Fig. 12(d)], which is consistent with the breakdown of the Duffing approximation and the description provided by the Maxwell-Bloch equations, due to the active participation of both quantum degrees of freedom. In most cases, the dark state follows the transition from the bright to the dim metastable state, with a lifetime which is much shorter than the duration of the two metastable states in the presence of spontaneous emission, yet significant in comparison to the cavity and qubit lifetimes, as Fig. 12 evidences (for a more detailed account on the switching, see Secs. IV A and IV B). In Fig. 13, we are following the thread from the appearance of the dark state up to the establishment of the bright metastable state, in a regime where the Maxwell-Bloch equations predict only the occurrence of the latter (see Fig. 8). However, as the ME results suggest, the dark state complex comes about together with a center and an unstable point in the dynamical evolution with changing relative position between them during the quantumactivated switching (see Fig. 16). The dim state reappears in Fig. 13(d) as the bright state is about to be reached, presenting itself as a relic of Maxwell-Bloch dispersive bistability. This finding supports the argument that the dark state is fragile and subject to rare quantum fluctuations, coexisting with the mean-field metastable states during the switching, and vanishing over the much longer lifetimes of the latter. The coexistence of the dark with the dim state is also verified by the ME results.
In light of the lengthening of the dark state lifetime in the absence of spontaneous emission, as shown in Fig.  11, the limit γ/(2κ) → 0 deserves a special consideration. For low drive strengths, the quasi-probability distribution for the cavity field remains essentially unchanged as γ → 0, pointing to the fact that the nascence of complex amplitude bistability is not related to the scale parameter γ 2 /(8g 2 ) → 0, as is the case for the absorptive bistability at resonance [33], but rather to δ 2 /(4g 2 ).
A deviation from the equilibrium configuration has already been pointed out for a multiple-atom saturable absorber on resonance, following an adiabatic elimination of the atomic variables [35], where a recourse to the Gaussian probability distribution is sought for the calculation of moments, apart from the transition region where the zero-delay second-order correlation function g (2) ss (τ = 0) diverges. In the dispersive regime, as the number of system excitations increases, not only are we unable to assume σ z = σ z ≈ −1, but the coupling to the environment and drive field (and consequently the entire ME), are rescaled to account for the actively participating system degrees of freedom (see Sec. III of [25]). The dark state appears as a result of joint quantum bistability, having roots in the region of the critical point C 1 where the qubit dresses the cavity with a weak nonlinearity (see the Appendix). The lifetime of this quasimetastable state is heavily dependent on the spontaneous emission rate, responsible for significant mixing between the various states participating in the switching. When γ = 0, transitions between the qubit states occur via the cavity through the Purcell decay, with a weaker mixing, resulting in a close to a hundredfold increase in the participation of the dark state in the dynamics (for a discussion on ladder switching at resonance, see Sec. V of [32]). Both for zero and non-zero spontaneous emission rates, the dark state persists past the bifurcation point of semiclassical bistability (the characteristic S-shaped curve), where the Maxwell-Bloch equations predict the sole presence of the bright state. With increasing drive strength the dim and the dark states exchange probability and the former eventually dissolves into quantum fluctuations (see Figs. 1 and 8).

A. Dark state for increasing drive strength
The dark state is characterized by very low photon occupation, and appears to follow closely the switching from the dim to the bright metastable state and vice versa (see for example Panel I of Fig. 8, as well as the Bloch-sphere phase portrait of Fig. 14). The associated histogram reveals a quasi-coherent state rather than a thermal state, as we can observe in Fig. 14. In sharp contrast with the predictions of the Maxwell-Bloch bistability, the dark state populates the upper half of the Bloch sphere. Figure 15 depicts a similar phenomenon to the one observed in Panel III of Fig. 8, where the dim state dissolves into the quantum fluctuations. In a detailed focus, the saddle point and the center of the dark state are shown in Fig. 16.
In our treatment so far we have linked the ME results to single quantum trajectories for the cavity and qubit, showing explicitly the appearance of the dark state. This quasi-metastable state is strongly related to bistable switching between the states of mean-field dispersive bistability. This is not a necessary condition though. In many instances (see Figs. 8 and 1), we have shown bistable switching for increased drive power, where eventually the dim state disappears and only the dark state remains. In that sense, the dark state can be considered on an equal footing as the metastable states of dispersive single-atom bistability, a pure result of quan-tum fluctuations (see Fig. 11 for γ = 0) and remaining last in the excitation spiral after the disappearance of the dim state with increasing drive strength (see frame (d) of Panel II in Fig. 8). With increasing spontaneous emission rate, it comes about as a rarer fluctuation state with a shortened lifetime testifying to its fragility with respect to decoherence. At resonance, the states of phase bimodality associated with the limit γ/(2κ) → 0 persist even in the presence of spontaneous emission (see Figs. 3 and 4 of [36]).

B. The neoclassical equations
The Maxwell-Bloch equations with γ = 0, also called the neoclassical equations [32], predict two states lying close to the two poles of the Bloch sphere, one stable and one unstable, both having a very low photon occupation. At resonance, near the limit of zero-system size, γ 2 /(8g 2 ) = 0, the neoclassical states and the states of absorptive bistability become structurally unstable [13]. A similar conclusion can be drawn for our case in the dispersive regime in the absence of spontaneous emission. We should point out here that in the dispersive regime, the above limit refers to the lower bound of the Purcell contribution (a second-order effect, see Section IV(B) of [37] for further discussion) γ P = κ(g 2 /δ 2 ) which is typically one to two orders of magnitude smaller than the linear cavity decay rate in the dispersive regime. Even in the case of resonance, the limit γ → 0 has only a formal meaning since in the absence of spontaneous emission no switching can occur between the JC excitation ladders (see Sec. 5 of [36]).
In the steady state, and since ∆ω q = 0, we obtain with |µ| = g|α| In the drive parameter regime under consideration, we observe that (the subscript d denotes the dark state) when |µ d | 2 < (g|α d |/|∆ω q |) 2 1, the following inequality follows: 7a], which suggests that the neoclassical states have an occupation below the level of one photon and also implies that one of the two neoclassical states is very close to the north pole of the Bloch sphere. As γ 2κ everywhere, the neoclassical state appears in the switching and coexists with the mean-field states of dispersive amplitude bistability. The Maxwell-Bloch states as well as the pair of dark states then become fragile to fluctuations as γ → 0, in a fashion similar to the resonant case (see the relevant discussion in Sec. 16.3 of [13], and [36]). In the former reference we read, "When γ is close to zero, the relaxation time to these states [the steady-states of absorptive optical bistability] becomes extremely long. The limit of zero system size is in this sense structurally unstable. It follows that near to this limit all of the mentioned states [i.e. steady-states of absorptive optical bistability and the neoclassical states] are quasi-stationary and fragile to fluctuations." Note also that the neoclassical solution with ζ > 0 is unstable with respect to fluctuations, which is yet another indication of the transient character of the dark state in our quantum simulations. The contour plots of the quasi -distribution functions evidence the presence of two distinct states in the phase space, with very low |α|, as predicted from Eqs. (7). The dark state can be verified experimentally through direct Wigner tomography [38] or via observing the qubit vector close to the north pole in the Bloch sphere for a time greater than max{1/(2κ), 1/γ}, noting at the same time the strong entanglement between the cavity and qubit.

V. CONCLUDING DISCUSSION
In this work we report on the appearance of a metastable state in the strongly dispersive regime, which is not predicted by the Maxwell-Bloch equations. We have investigated the rôle of quantum fluctuations, which induce bistable switching in the driven dispersive Jaynes-Cummings model with weak spontaneous emission, in the appropriate region of the drive strength and frequency where the Maxwell-Bloch equations predict steady-state bistability. The breakdown of the Duffing approximation and the appearance of terms that are higher order than quartic in the field operators multiplied by qubit operators certainly suggest that an FPE cannot be formulated, with the qubit playing a very active rôle in the cavity nonlinearity. While some typical instances of quantumfluctuation switching, such as the decay of the mean-field unstable state, manifest as statistical mixtures of semicoherent cavity photon states with varying weights, can we expect that the two-level atom will manifestly break the classical picture for increasing drive powers within one quantum trajectory? The appearance of the dark state seems to yield a preliminary "yes", highlighting the importance of the neoclassical equations combined with quantum fluctuations that are responsible for organizing the asymptotic dynamics when γ/(2κ) 1. The origin of this state brings us closer to the low-excitation dispersive regime (below the critical point C 1 ), where the qubit can be considered a 'spectator' and is not actively involved in the switching dynamics. Fragile to fluctuations, the dark state is intimately linked to the qubit-cavity interaction, since the entanglement entropy of the two oscillators increases drastically while the state lasts in the trajectory. Inasmuch as its lifetime is concerned, it may be deemed a quasi-metastable state coexisting with the states of dispersive bistability, which also reveals itself after the various quantum trajectories have been averaged, in contrast to the unstable mean-field state.
The data underlying this work is available without restriction [39].  In the Appendix we derive the basic results applying to the Duffing approximation and the associated Wigner distribution we present in Section III. In terms of the generalized P-representation we can write (A1) Substituting the steady state solution for the Duffing oscillator [6] Recognizing the following integral representation of the Bessel function: with C being the Hankel path starting at −∞, encircling the origin in an anticlockwise fashion, and returning back to −∞, the final result reads [27,28] W (α) = N e −2|α| 2 J c−1 ( (A5) (N is the normalization constant). We note that the function is everywhere positive. We normalize it through the condition which leads to (with c = d * ) (A7) The integral over φ evaluates to 2πδ kl and the integral over ρ yields Taking into account these two results we finally arrive at .
(A9) We will now prove that the probability density function p(n) is normalized. From the given steady-state photon number probability density function, we have ∞ n=0 p(n) = S 0 with: But the sum inside the brackets equals unity, by virtue of the normalization of the binomial distribution with p = 1/2. Hence, the above sum reads with so that finally We will now use the Wigner quasi-distribution to calculate the symmetrically ordered operator moments (a † ) n a m S as follows: (a † ) n a m S = 2 π Γ(c)Γ(d) Despite the fact that the Wigner function for every steady state of the driven dissipative Duffing oscillator is positive, the departure from a Gaussian distribution is obvious in the region of bistability, which is also reflected in the expressions for the various moments of the intracavity field. The latter are highly nonlinear functions of the drive strength, as is also verified in Fig. 1 of [6]. The corresponding steady-state photon number probability distribution function can be written as .
(A18) The above expression shows the deviation from the Poissonian distribution of a coherent state, with increasing drive power.
In Fig. A.1, we display the Wigner quasi-distribution functions for the intracavity photon field (with α = x + iy) in a driving region where the qubit is not significantly excited and we can set σ z = σ z = −1. The exact master equation predictions and the Duffing approximation of Eq. (4) are in good agreement, showing the development of low amplitude bistability alongside the departure from the Gaussian shape of a coherent state.

(c) and
A.1(d), the mean-field analysis with 2κ/γ = 12 predicts only one state with n ss ≈ 1.56, captured in the Wigner function plots by the peak of the squeezed state centered at α ss ≈ 1.025 − 0.92i, exhibiting negligible variation with changing γ. The last frame of Fig. 1 shows that the dim state is clearly distinguished from the dark-state pair. This separation occurs at a drive strength for which the Maxwell-Bloch equations predict the existence of the bright state only, similarly to the situation we have encountered for the low-amplitude bistability approximated by the effective Duffing nonlinearity, as we have seen in Fig. A.2.