Cavity-field distribution in multiphoton Jaynes-Cummings resonances

We calculate the cavity-field distribution in the Wigner representation for the two-photon resonance of the weakly driven Jaynes-Cummings (JC) oscillator in its strong-coupling limit. Using an effective four-level system, we analytically demonstrate the presence of steady-state and transient bimodality which breaks azimuthal symmetry in phase space. The two steady-state peaks are located at opposite positions and do not correspond to the two-photon amplitude of the driven transition. The developing bimodality is portrayed in parallel with the evolution of the intensity correlation function for the forwards-scattered photons, before being finally contrasted to the few-photon steady-state and transient phase-space profiles for the cavity mode in the JC model driven on resonance.


I. INTRODUCTION
The two-level behavior of an optical cavity mode strongly coupled to a single atom and excited near a vacuum Rabi resonance has been demonstrated about three decades ago as a paradigmatic system relying on the √ n nonlinearity of the Jaynes-Cummings (JC) model and exhibiting photon blockade [1]. This effect was reported a few years after single-atom absorptive optical bistability had been ascertained in a regime lying at the "interface between the quantum limit, in which quantummechanical noise invalidates the semiclassical prediction of bistability, and the classical limit, in which quantum noise is a negligible perturbation on semiclassical results" [2]. Only a handful of years until nowadays have lapsed, however, since the experimental observation of two-photon blockade in the strong-coupling limit of cavity QED [3] and the demonstration of the photon blockade breakdown in circuit QED by means of a first-order dissipative quantum phase transition [4]. This experimental verification of criticality in zero dimensions has been very recently followed by a theory for the critical scaling of steady-state observables as the critical point of the second-order phase transition associated with photon blockade [5] is approached [6]. Apart from its rather obvious link to photon antibunching and sub-Poissonian statistics [7][8][9][10][11][12], which arise as distinct quantum features, photon blockade has also been associated with a more intricate property, namely the existence of a strong-coupling "thermodynamic limit" in single-atom QED, where the relevant system size parameter -revealed for a vanishing spontaneous emission rate -is proportional to the square of the light-matter coupling strength. Upon approaching this limit, quantum fluctuations are not a"negligible perturbation on semiclassical results". Instead, multipho- * Email address: th.mavrogordatos@gmail.com ton resonances of higher order appear, turning increasingly sharper before they saturate for a growing drive amplitude, while the discrepancy from the semiclassical nonlinear response persists (see also Sec. IV-C of [5]).
In their experimental report on the properties of the two-photon resonance for one trapped atom in a highfinesse optical cavity, Schuster and collaborators mention that they "also analyzed the nonlinear theory of optical bistability and found that it is inconsistent with all of the measurements presented", identifying the operation on the lower branch of bistability where the nonlinear response is weak [13]. Shortly after, in one of the exemplary experiments of circuit QED, spectroscopic evidence for multiphoton dressed states along the JC ladder was provided by heterodyne transmission measurements in a setup using a superconducting qubit coupled to a microwave field [14]. Therein, Bishop and collaborators appeal to an effective two-level system to model the supersplitting of the vacuum Rabi resonance, bringing us back to the early ideas on photon blockade developed in [1] and their two-mode extension implemented in the experiment of [11].
Bringing these pieces together, in our work we aim to bridge the gap between the aforementioned quantum and classical limits focusing on a multiphoton resonance, a configuration where the classical nonlinearity is intertwined with quantum fluctuations in a fundamental way that is pictured in phase space. For the special case of the two-photon transition the link is provided by lowamplitude quantum bistability, or rather bimodality, evidenced in the quasiprobability distribution of the intracavity field in the strong-coupling regime. After reducing our description to a minimal four-level model in Sec. II, we first calculate the steady-state Wigner function at quasiresonant excitation of the transition, and we then relate the transient bimodality to the evolution of the intensity correlation function for the forwards scattered photons in Sec. III. We then find in Sec. IV that the multi-modality accompanying a multiphoton resonance disappears when the JC oscillator is resonantly excited with a weak drive amplitude. Instead, in the nonperturbative treatment one speaks of two counter-propagating wavepackets in close connection with the initial Gaussian shape of the state naturally produced by the coherent driving. In the last part of our discussion, reserved for the Conclusions, we point to a possible experimental observation of the change occurring in the multiphoton quantum nonlinearity when we transition from the photon-blockade regime to resonant excitation.

II. FORMULATING THE EFFECTIVE MASTER EQUATION
The system density matrix ρ for the open driven JC model obeys the Lindblad master equation (ME) where a and a † are the annihilation and creation operators for the cavity photons, σ + and σ − are the raising and lowering operators for the two-level atom, g is the dipole coupling strength, 2κ is the photon loss rate from the cavity, and γ is the spontaneous emission rate for the atom to modes other than the privileged cavity mode which is coherently driven with amplitude ε d and frequency ω d defining the detuning ∆ω d ≡ ω d −ω 0 . Hereinafter, we consider the special case γ = 2κ which considerably simplifies the analysis, while we take g/κ = 1000 after [14,15]. We now employ an effective four-level model to produce an analytical expression for the matrix elements of the system density matrix when the two-photon transition is resonantly excited in the weak excitation regime. Such a perturbative treatment, devised in [15] to study the first and second-order coherence of the forwards scattered light for the two-photon resonance in one-atom cavity QED, uses essentially an expansion in powers of (ε d /g) 2 (without loss of generality for this derivation we assume that the driving field amplitude ε d is real). The first four JC dressed states between which transitions are constrained to occur, starting from the ground state, are denoted by where |n, ± ≡ |n ⊗ |± , |n is the Fock state of the cavity field, while |+ , |− are the upper and lower states of the two-level atom, respectively. The secular approximation in the limit of strong nonperturbative coupling (g κ, γ/2) leads to the following effective ME in the laboratory frame [see [16] for further details, and Eq. (18) of [15]] where the effective Hamiltonian of the four-state model is with the following energies (shifted by δ 0 -δ 3 ) for the four states dressed by the drive, and Ω = 2 √ 2ε 2 d /g the two-photon Rabi frequency to dominant order [16]. In the effective ME of Eq. (3), we define as usual are the transition rates between the dressed states. The shifted dressed-state levels and the transition rates between them comprise the effective four-level model depicted in Fig. 1. From Eqs. (5a) and (5d) we find that the two-photon resonance, including the level shifts, must be excited with a drive frequency ω d given by including second-order corrections with respect to the drive amplitude.

III. FROM THE CAVITY-FIELD DISTRIBUTION TO THE COHERENCE OF THE FORWARDS-SCATTERED PHOTONS
The effective ME (3) can be now used to evolve the system density matrix from a given initial state. For a general time τ , the reduced density matrix for the cavity field, defined as ρ c ≡ +|ρ|+ + −|ρ|− , can be written in the form where the matrix elements are taken with respect to the dressed states |ξ 0 -|ξ 3 , while we omit their time arguments for brevity. The time dependence of the matrix elements featuring in the above expression can be calculated as follows, starting from featuring two constants to be determined from the initial conditions with Σ ≡ ρ 33 (0) − ρ 00 (0) and C expressed in terms of Σ via ρ 33 (0). We then go on to determine the sum of occupation probabilities comprising yet another constant to be determined from ρ 11 (0) + ρ 22 (0). Now, ρ 00 follows from normalization, while the off-diagonal elements involved read The two complex conjugate matrix elements in Eq. (11a) generate a quantum beat at the frequency ν = 2g+δ 2 −δ 1 [see Eqs. (5b) and (5c), as well as Eq. (19)] which is significantly higher than all the other rates involved in the dynamics, and is revealed by the intensity correlation function of the transmitted light, as we will see later on. With these results in hand, we can now propagate to the steady state. Our final expressions for the steady state are given in terms of the excitation probability of level |ξ 3 -the upper level of the two-photon transition -denoted by The Wigner distribution W (α, α * ) (in the complex variable α = x + iy and its conjugate, α * = x − iy) is the Fourier transform of the symmetrically-ordered characteristic function (see Ch. 4 of [17]) and is defined as (13) A central result to guide our discussion is the normalized Wigner distribution of the (intra)cavity field in the steady state, which we derive in the form The two terms featuring in the last line on the RHS are responsible for breaking the symmetry arising from a linear combination of the vacuum-state and Fock-state (|l , l = 0, 1, 2) phase-space distributions (see Ch. 3 of [17]). Symmetry breaking originates from the offdiagonal element |0 2| and its Hermitian conjugate in ρ c , coming directly from the driving terms in the Hamiltonian (4). As the two-photon Rabi frequency increases, the steady-state distribution of Eq. (14) develops two peaks along the line y = −x and two dips along the line y = x in phase space (clearly forming for p 3 0.15). The steady-state Wigner function remains everywhere positive for 0 ≤ p 3 ≤ 1/4.
In Fig. 2, we trace the onset of bimodality for increasing p 3 with reference to our effective four-level model, depicted in Fig. 1, to which the dynamical evolution is constrained. The analytical expression for the Wigner function suggests that azimuthal symmetry is restored asymptotically as the two-photon transition saturates and the occupation probability p 3 tends to 1/4 -its maximal value. In fact, this picture reflects the limitations of the four-level model: as we can observe in Panel I of Fig. 2, which depicts numerical results for the steadystate intracavity field distribution, the two peaks already deviate from the line y = −x for p 3 = 0.24 (being located at α m1 ≈ 0.35 − 0.53i and α m2 ≈ −0.56 + 0.45i in contrast to the analytical distribution placing the peaks at α m1,2 ≈ ±0.48 ∓ 0.48i). Agreement between exact numerical results and the analytical distributions shown in Panel II is also inherently compromised by the approximations involved in deriving Eq. (3) from (1). Upon a further increase of the drive amplitude, when the perturbative expansion breaks down and higher levels along the JC ladders enter into play, the two peaks approach each other, as we can observe in frame (d) of Panel I, until they eventually merge into a single broad squeezed distribution. Qualitatively, we can expect that the n-photon resonance will bring about a multimodality of order n through the term proportional to [α n − (α * ) n ]e −2|α| 2 . For ∆ω d ≈ −g/ √ 3, numerical simulations confirm the existence of three peaks in phase space with neighboring position vectors forming an angle of 2π/3 between them. This is in agreement with the steady-state Wigner functions given in Fig. 8 of [18] for the n-photon blockade encountered in the driven dissipative Kerr oscillator; the distribution exhibits n peaks and n dips.
The Wigner function of Eq. (14) can be used for the calculation of various moments of the cavity field in the steady state, respecting symmetric operator ordering. For example, we compute the steady-state average photon number as whence we verify that a † a ss = 5p 3 /2. Only the first three terms in Eq. (14), those originating from the vacuum state and the Fock states with l = 1, 2, contribute to the final result -the symmetry-breaking terms multiplied by |α| 2 drop out. Phase-space bimodality is also present in the transient regime. Following [15], we select the initial density matrix as the conditional state occurring at the emission of one photon from the cavity after the attainment of steady state, with the superposition generating a quantum beat of partial visibility, This particular choice determines the set of constants and ultimately the time-varying density matrix from which we derive the quasiprobability distribution of the intracavity field. It also determines the second-order correlation function for the forwards-scattered light, given in [15] as with coefficients depending solely on p 3 through the ratio γ/Ω: The quantum beat -the term with coefficient a 4 in Eq. (19) -is averaged out in Fig. 3(a), leaving behind an oscillation whose period is determined by the two-photon Rabi frequency here scaled to Ω/γ ≈ 4.54 (corresponding to ε d /γ ≈ 28.32). In Figs. 3(b)-(d), we observe the evolution of the Wigner function from zero time delay to the time when g (2) (τ ) attains its maximum value. We encounter transient symmetry breaking, with the peaks of the cavity distribution placed at different positions along the lines y = ±x in phase space, always at equal distances from the origin. The conditional state of Eq. (16) has a circularly symmetric Wigner distribution, as we can see in Fig. 3(b), which is subsequently broken in the early stages of the evolution. It corresponds to photon antibunching, since g (2) (0) ≈ 0.65 < 1 and the subsequent evolution has a positive slope. The depicted curve satisfies the criterion g (2) (0) = 0, g (2) (0) > 0 [9,10]. Interestingly, the form of Eq. (19) predicts photon bunching instead at low drive amplitudes (whence for low p 3 ), as noted in [15]. The maximum value of g (2) (τ ) is associated with a bimodal profile comprising two peaks separated by a region around the origin where the Wigner function is negative [ Fig. 3(c)]-a remaining sign of the single-photon distribution originating from the dressed states |ξ 1 and |ξ 2 involved in the quantum beat. We note that the region separating the two peaks is directly affected by the period of the beat during the evolution (2π/ν). For instance, in the Wigner function plotted for γτ = 0.3580, half a beating period away from the location of the maximum of g (2) (τ ), the two peaks are joined by a pronounced ridge, and the distribution is everywhere positive. Finally, we find that the bimodality along the line y = x [ Fig. 3(d)] occurs during the time intervals when Im(ρ 03 ) < 0. The existence of such intervals indicates a saturated two-photon transition (p 3 > 0.24).

IV. JC INTERACTION DRIVEN ON RESONANCE: FROM THE SECULAR APPROXIMATION TO QUANTUM-STATE REVIVALS
Let us now move to the weak-excitation limit on resonance, where the evolution of the system density involves primarily the first three levels of the JC ladder |ξ 0 , |ξ 1 and |ξ 2 , according to the perturbative analysis of Sec. 16.3.4 in [19] within the secular approximation relying once more on what has been termed "dressing of the dressed states". In this weak-excitation regime, to dominant order, the Wigner function of the intracavity field has azimuthal symmetry with a leading contribu- tion from the vacuum state and subsequently from the the Fock states |l l| with l = 1, 2. The corresponding coefficients of these three terms are expanded in powers of the squeezing parameter r ≡ (|ε d |/g) 2 1 to dominant order as 1 + r + r 2 , r 2 and r 2 /2, respectively. We also find that the reduced density matrix for the cavity field contains the term (1/ √ 2)[(ε * d /g) 2 |0 2| + (ε d /g) 2 |2 0|] which produces a contribution proportional to [(ε * d /g) 2 α 2 + (ε d /g) 2 (α * ) 2 ]e −2|α| 2 in the Wigner function, breaking once again the symmetry in phase space in a direction determined by the phase of ε d . In contrast to the two-photon resonance, however, these terms yield a negligible contribution since they are compared against the unity coefficient coming from the vacuum state; higher-order resonances have virtually no effect in the perturbative expansion. Cavity distributions obtained from the numerical solution of Eq. (1) confirm the prevalence of the vacuum state for the set of ratios |ε d |/g used in Fig. 2.
To obtain a non-perturbative analytical result involving transient bimodality at the few-photon level, we make a detour to briefly visit the topic of collapse and revival in the JC model, and consider a closely related nondissipative variant where now the two-level atom is coherently driven instead of the cavity mode (for the open version of the problem, the equivalence between atom and cavity driving is elucidated in [22]). The interactionpicture Hamiltonian for this system is [20] where ε d is the amplitude of the external field, assumed to be real. Introducing the displacement with D(α 0 ) ≡ exp[α 0 (a − a † )], we obtaiñ The system state can be written as a sum of two orthogonal contributions, representing excitations of almost orthogonal anharmonic oscillator modes, with the two orthogonal states expressed as where Ξ = (U, L). The coefficients c n ≡ e −α 2 0 /2 (α n 0 / √ n!) come from the familiar expansion of a coherent state into the Fock-state basis. Here, though [for n ≥ 1 and adopting a slightly different notation to Eqs. (2)], are the excited-state doublets of the JC ladder (with |0 ⊗ |− the ground state).
To visualize the evolution in the phase space we derive the Q function for the intracavity field in the transformed frame, defined as Q(z,z * , t) = (1/π)[ z|ψ (t) ψ (t)|z ] and normalized to unity, from the pure state of Eq. (25). The resulting expression reads (z =x + iỹ) where d n ≡ e −|z| 2 /2 (z n / √ n!). The evolving state of Eq. (25) identifies two counter-propagating wavepackets depicted in the two insets of Fig. 4, which shows the evolution of the intracavity photon number a † a(t) (in the original reference frame) from the initial state |0 ⊗ |+ . The two peaks observed in inset A are located atz A =x A + iỹ A ≈ 0.28 ± i1.70, with |z A | 2 being very close to α 2 0 = 3. Remnants of these peaks, depicted in the inset B, are still visible at the "revival" time when the peak is located atz B ≈ − √ 3i. We note that for high excitation (α 2 0 1), where the time evolution of the system state is described by the asymptotic analysis of [23], the Q function-depending on the scaled time gt-can be simplified as showing explicitly the two counter-propagating contributions in interference with each other, through the Euler formula applied to the cosine.

V. CONCLUDING REMARKS
To summarize, in this article we have revealed the distinct character of multi-modality accompanying a multiphoton resonance in a regime where quantum dynamics produces a visible departure from the mean-field predic-tions. For the special case of the driven two-photon resonance, we have shown that the two peaks of the bimodal distribution are bound by the unit circle for every drive amplitude for which the perturbative treatment remains valid, and are not in correspondence with the two-photon amplitude, i.e., a peak of the distribution at a position vector of magnitude equal to √ 2. If bistability is to be regarded as a feature originating from a semiclassical nonlinearity, then it ties well with the Rabi oscillation developing in the intensity correlation function for the forwards scattered photons as the two-photon transition saturates. At the same time, bimodality is present when we operate well within the strong-coupling regime; as the system-size parameter [∼ (γ/g) 2 ] of absorptive bistability -characteristic of a weak-coupling "thermodynamic limit" -tends to zero for given ratios ε d /g and γ/(2κ), the analytical expression for the Wigner function suggests an increasingly symmetric cavity-field distribution (p 3 → 1/4). Moreover, numerical simulations show a bimodal distribution for the parameters of Fig. 2(c) (Panel I) but in the absence of spontaneous emission [γ/(2κ) = 0]. While bimodality persists in this so-called "zero system size" where the relevance of a system-size parameter [∼ (g/κ) 2 ] attributed to a strong-coupling "thermodynamic limit" emerges, we encounter an asymmetric distribution where most of the excitation probability gathers in the second quadrant (with the main peak located at α m1 ≈ −0.64 + 0.39i).
On the experimental front, the phase-space profile of the cavity field in the strong-coupling limit will assist the realization of light sources producing bound states of multiple photons [24] in conjunction with autocorrelation functions of order higher than second, to provide a clearer indication of the nonclassical character of the emitted light at weak excitation [25]. Moreover, the direct measurements of the Wigner function close to the origin of the phase space, in the spirit of [26], can prove a powerful tool for assessing the multiphoton quantum nonlinearity associated with photon blockade. The proposal has been implemented in the experiment of [27] for a single-photon field. The latter is relevant for detecting the alternation between positive and negative values of the Wigner function close to the phase-space origin, correlated with the quantum beat. Finally, the scheme of [26] allows one to select the region of phase space to be explored, while being insensitive to a low detection efficiency. Such a direct measurement would readily reveal the distinct features of a multiphoton resonance in the JC oscillator, and uncover a critical behavior which is markedly different than the operation on resonance.