Critical Theory for the Breakdown of Photon Blockade

Photon blockade is the result of the interplay between the quantized nature of light and strong optical nonlinearities, whereby strong photon-photon repulsion prevents a quantum optical system from absorbing multiple photons. We theoretically study a single atom coupled to the light field, described by the resonantly driven Jaynes--Cummings model, in which case the photon blockade breaks down in a second order phase transition at a critical drive strength. We show that this transition is associated to the spontaneous breaking of an anti-unitary PT-symmetry. Within a semiclassical approximation we calculate the expectation values of observables in the steady state. We then move beyond the semiclassical approximation and approach the critical point from the disordered (blockaded) phase by reducing the Lindblad quantum master equation to a classical rate equation that we solve. The width of the steady-state distribution in Fock space is found to diverge as we approach the critical point with a simple power-law, allowing us to calculate the critical scaling of steady state observables without invoking mean-field theory. We propose a simple physical toy model for biased diffusion in the space of occupation numbers, which captures the universal properties of the steady state. We list several experimental platforms where this phenomenon may be observed.

Photon blockade is the result of the interplay between the quantized nature of light and strong optical nonlinearities, whereby strong photon-photon repulsion prevents a quantum optical system from absorbing multiple photons. We theoretically study a single atom coupled to the light field, described by the resonantly driven Jaynes-Cummings model, in which case the photon blockade breaks down in a second order phase transition at a critical drive strength. We show that this transition is associated to the spontaneous breaking of an anti-unitary PT -symmetry. Within a semiclassical approximation we calculate the expectation values of observables in the steady state. We then move beyond the semiclassical approximation and approach the critical point from the disordered (blockaded) phase by reducing the Lindblad quantum master equation to a classical rate equation that we solve. The width of the steady-state distribution in Fock space is found to diverge as we approach the critical point with a simple power-law, allowing us to calculate the critical scaling of steady state observables without invoking mean-field theory. We propose a simple physical toy model for biased diffusion in the space of occupation numbers, which captures the universal properties of the steady state. We list several experimental platforms where this phenomenon may be observed.
Deviations from the equilibrium ensemble are particularly common in the field of quantum optics [41][42][43][44][45][46][47][48][49][50][51][52][53][54]. Since the typical temperature of a black-body that emits in the optical frequency range (roughly 500 THz) is of order 4800 K, almost all experiments involving optical photons are conducted far from the equilibrium temperature scale [55]. Furthermore, processes of interest often * jcurtis1@umd.edu involve strong (coherent) driving, e.g. by lasers, in part to overcome the losses due to the fact that the system of interest is typically coupled to a highly incoherent environment. This combination of strong driving and incoherent loss processes often results in a quantum system which is far from thermal equilibrium.
Non-equilibrium effects are especially pronounced in dynamics of systems with strong interactions. For instance, strong optical nonlinearities give rise to the photon blockade [24,25,27,36,[56][57][58][59][60], whereby strong photon-photon repulsion inhibits the absorption of more than one photon, even in the presence of strong external driving. In this case, the effective single-occupancy constraint becomes readily apparent in photon transport through the device.
Photon blockade can be understood qualitatively by considering a single-mode of an electromagnetic resonator with non-linear spectrum E n = ω 0 n + 1 2 U n(n − 1) as a function of photon number n (i.e. Hubbard term in the case of a lattice). To lowest order in perturbation theory, coherently driving the photon field at frequency ω will connect a Fock state with n photons to a state with n + 1 photons, resulting in a correction to the eigenstates with an energy denominator ∼ E n+1 −E n −ω. For the Kerr-type of non-linearity described above, this can result in at most one divergent term, which occurs when the drive frequency is tuned to satisfy ω − ω 0 − U n = 0. The net result of this effect is that for sufficiently strong non-linearities U , the driven resonator effectively becomes a two-level system, indicating a strong photon anti-bunching and hence a departure from the thermal ensemble with Poissonian photon number statistics.
A convenient way of engineering a strong optical nonlinearity is to couple the resonator photon to an atom, which acts as a strongly non-linear "hard-core" boson. If we describe the system using the basic Jaynes-Cummings model [61] (see Sec. II), we find the celebrated √ nnon-linearity, resulting in the excitation spectrum E n ∼ ω 0 n ± g √ n, with g the atom-photon coupling constant.
As first observed in Ref. [46], this has important implications for the fate of the photon-blockade. Re-tooling the argument from the previous section, we see that there are divergences in the perturbative expansion when the driving frequency satisfies ω −ω 0 ±g √ n + 1 − √ n = 0.
For ω = ω 0 , the non-linearity vanishes for n → ∞ as n −1/2 . This implies that the driven oscillator can effectively tunnel off to a high-photon number state, for which the non-linearity is less important.
Remarkably, it appears that the breakdown of the photon-blockade proceeds via a continuous dissipative phase transition [46,62]. More specifically, this phase transition presents as a non-analyticity in the steady state of the Lindblad equation which governs the drivendissipative dynamics of the system. Since the system we consider is a single atom coupled to a single-mode of the electromagnetic field, it constitutes a system with no spatial extent. It is then not clear a priori what the relevant "system-size" parameter is, and what the appropriate thermodynamic limit is. Since the Fock space of the photon is formally infinite, we can see that nonanalyticities may arise if the entirety of the Fock space is accessible. As we discuss in more detail later, this identifies the thermodynamic limit with the limit of vanishing dissipation [54,63,64]. It is the purpose of this work to provide a comprehensive analysis of the nature of this critical point. The main results can be found in Tab. I, which presents the predicted critical scaling of various steady-state expectation values of observables as the critical point is approached. The approach to the critical point is controlled by the dimensionless parameter = 2E/g, which measures the drive strength E relative to the strength of the atom-photon coupling constant g. Criticality occurs at 2 = 1, with 2 < 1 corresponding to the disordered phase.
This work is structured as follows. In Sec. II, we discuss the model in more detail, focusing on the coherent portion of the evolution. In Sec. III, we present a semiclassical treatment of the problem and a brief description of the physical intuition behind the critical point. The main calculation is contained in Sec. IV, where we use a rate equation to solve for the steady-state behavior. Finally, in Sec. V, we use the results of Sec. IV to compute the critical scaling of physical observables near the critical point. We conclude in Sec. VI with a discussion of the implications of our calculation as well as some future directions of interest. More technical aspects of the analysis are presented in the appendices.

Observable
Scaling TABLE I. Scaling of observables as we approach the critical point. The parameter which controls the distance to the critical point is = 2E/g, with E the coherent drive strength and g the atom-photon coupling strength. The critical point is located at 2 = 1, with 2 < 1 the disordered phase. This applies in the thermodynamic limit, which in this case is obtained by taking the zero-dissipation limit of the steady-state ensemble (see Sec. IV).

A. Hamiltonian and master equation
The model we consider in this work is that of an atom coupled to a single electromagnetic mode of a cavity. The atom is modeled as a two-level system, with ground state |↓ and excited state |↑ . For simplicity we assume, as in Ref. [46], that the bare atomic and cavity frequencies are tuned to resonance with each other. Generically, the critical theory describing the breakdown of photon blockade exhibits a first-order phase transition when the detuning of the cavity drive is tuned away from being on resonance with the atom and cavity frequencies. Only when all three frequencies are resonant does the system exhibit a second order phase transition. We restrict our focus to the continuous phase transition and assume that all three of the bare frequencies are tuned to resonance at ω 0 . We model the atom-photon interaction by a simple dipole transition with coupling constant g and assume the rotating-wave approximation to be applicable.
The intrinsic coherent dynamics is described by the celebrated Jaynes-Cummings model [61,65] with Hamil-tonianĤ Herein,â is the photon annihilation operator andσ ± = 1 2 (σ x ±iσ y ) are the atomic transition operators written in terms of Pauli matrices. In order to model the external coherent driving, e.g. by a laser imposed onto the cavity, we add the term to the Hamiltonian. This corresponds to a monochromatic driving of the cavity photon field at resonance frequency ω 0 with strength E > 0. We incorporate dissipation through single-photon loss from the cavity to a zero temperature reservoir at rate 2κ. The dynamics of the density matrix is then given by the open-system Lindblad master equation We are particularly interested in the steady-state solutions of this equation. Let us briefly review the diagonalization of the Jaynes-Cummings HamiltonianĤ 0 . Due to the rotating-wave approximation, the "polariton" numberN =â †â +σ +σ− is conserved. Thus, the eigenspectrum ofĤ 0 separates into a direct sum of decoupled two-level systems, consisting of entangled atom-matter excitations known as polaritons. The energy eigenstates ofĤ 0 are labeled by a single signed quantum number ν = 0, ± √ n with n = 1, 2, ... and read These states have energies E = ω 0 n ± g √ n, respectively. In the following, we omit the direct product symbol.
In order to understanding the steady state of Eq. (3), it is helpful to first focus on the coherent evolution generated byĤ 0 +Ĥ d (t), neglecting dissipation. Due to the time-dependent driveĤ d (t), neither energy nor polariton number is conserved. However, we may still utilize the U (1) symmetry generated byN to analyze the system. By applying the unitary operatorÛ d (t) = exp(−iω 0 tN ), we can go to a frame co-rotating with the drive frequency. The resulting rotating-frame Hamiltonian is time-independent and, because the drive is resonant, linear in the photon operators. (When the drive is nonresonant, there is a quadratic photon term reflecting the splitting due to the finite detuning.) As first observed in Ref. [42] and later clarified in Ref. [43], the rotating frame HamiltonianĤ rf develops an instability as the drive strength E is increased whereby a discrete spectrum at weak driving gives way to a continuum at strong driving. This transition occurs for all eigenstates in the spectrum simulataneously and may be understood as arising from the competition between the termâ †σ − +âσ + , which can be diagonalized by the polaritonic eigenstates from Eq. (4), and the term a +â † ∝X, which is proportional to the position operator and has no normalizable eigenbasis (the eigenfunctions of theX operator are Dirac delta functions which cannot be properly normalized). Accordingly, the discrete, quantized spectrum prevails when g/E 1, while the non-normalizable continuum emerges for E/g 1. The critical point where the spectrum becomes continuous occurs at the critical drive strength To see this, note that for E < E c the Hamiltonian is diagonalizable and eigenstates can be found exactly [43].
(We rederive this result in Appendix A). In this regime, the eigenvalues retain their √ n-like spacing even for E = 0 and are given by indexed by the signed quantum number ν = 0, ± √ n introduce above. Here we introduce the dimensionless drive strength parameter As we approach the critical point from below, the effective coupling that controls the level spacing is g eff = g 1 − 2 3/4 . Crucially, for 2 = 1, the level spacing collapses to zero, and the discrete spectrum condenses into a continuum [43]. An alternative description of this transition can be found in Ref. [66], where the eigenvalue problem is mapped on to that of a charged Dirac particle in both electric and magnetic fields. Above the critical drive strength, the HamiltonianĤ rf is no longer diagonalizable and exhibits dynamical instability. In this case, the cavity dissipation is crucial, since it ultimately limits the photon number in the absence of detuning.

B. Symmetries
In this section, we discuss the important role of symmetries in our model. Though the cavity driving provides a preferred reference phase and thus destroys the U (1) symmetry, there is still a remnant anti-unitary Z 2 discrete symmetry associated with the rotating-frame HamiltonianĤ rf [42]. This corresponding transformation acts jointly on the photon and atom in an anti-unitary fashion and is of the form where P = e πiâ †â is the bosonic parity operator, which implements "spatial" inversion on the bosonic mode while acting trivially on the atomic degree of freedom, and T is time-reversal. Since the two-level system is not a real spin, but a pseudo-spin, T acts on the atomic degree of freedom through complex conjugation alone. Hence T 2 = +1, and there is no Kramers degeneracy. Using the quadrature representation of the bosonic ladder operator asâ =X +iP √ 2 , we have The first of these relations is nothing but the statement that C is anti-unitary. Under these transformations, the rotating-frame Hamiltonian obeys the particle-hole type symmetry This implies that to each eigenstate |E of H rf with energy E > 0 there exists an eigenstate | − E = C |E with energy −E. For E < E c , these are, of course, the eigenvalues displayed in Eq. (7), but the analysis presented here reveals that the particle-hole symmetry also exists for E > E c . Note that, for E = 0, the operator C relates the two "polariton branches" in Eq. (4). Indeed, using that C |n = (−1) n |n for the noninteracting bosonic Fock states, we see that the transformation inverts the spectrum by switching between the polariton branches for n > 0 according to C |ν 0 ∝ | − ν 0 . The noninteracting vacuum state |0 with E = 0 is automatically particle-hole symmetric. We find that, while the transformation leaves the ladder-operators intact (apart from an overall phase), it acts non-trivially on coherent states built from them. In particular, for a bosonic coherent state the conjugated state is C |α = | − α * . The critical point investigated in this work is associated with a spontaneous breaking of this Z 2 symmetry in the nonequilibrium steady state. Potential order parameters are the operators X ∝ Re â and σ y , both being odd under the transformation C . In the next section, we compute the steady-state behavior of these quantities using a semiclassical approximation and show that indeed they ought to exhibit nonanalytic behavior at the critical point.

III. MEAN-FIELD THEORY
In this section, we analyze the master equation (3) in terms of a mean-field or semi-classical theory. In this scheme, we compute the time evolution of the expectation values under the approximation that correlations factorize according to Hereˆ σ = (σ x ,σ y ,σ z ) T is the vector of the atomic Pauli matrix operators. Under this assumption, the equations of motion for the expectation values close, and we have Notably, these equations conserve the length of the atomic Bloch vector Since this quantity is conserved, we have a one-parameter family of steady states labeled by the length of the Bloch vector . This multiplicity is not expected to survive once we include fluctuations, and is most likely an artefact of the approximation in Eq. (16). The steady-state solution of equations (17) is obtained by setting the time-derivatives to zero and solving the resulting algebraic equations. For a fixed value of 2 , the solution to these equations depends on whether the drive is above or below the critical value of 2 c = 2 . In the following, we examine these solutions in more detail. A summary of our findings is presented in Tab. II.
For 2 < 2 (below the critical drive strength), the solution to the steady-state equations is In this case, we may interpret the atomic Bloch vector ˆ σ as a dipole re-radiating a coherent field which selfconsistently cancels the externally imposed driving field, producing a total photon field of â = 0. This steady state corresponds to the "disordered" phase with respect to the Z 2 symmetry induced by C . Indeed, the expectation values which transform non-trivially under C , are both zero in this regime. Importantly, the solution is independent of κ (or, more precisely, independent of the dimensionless parameter g/2κ). Furthermore, a simple analysis of the non-linear equations of motion around this steady state reveals that the solution with σ 3 < 0 is a dynamically stable fixed point, while σ 3 > 0 is dynamically unstable (see Appendix C for the analysis). For 2 > 2 (above the critical point), we find the steady-state solution to form a new pair of dynamical The two solutions labelled by "±" are related by the action of C . In these steady states, both Re â and σ y are nonzero, thereby spontaneously breaking the Z 2 symmetry. The order parameter characterizes the continuous nonequilibrium phase transition. Let us remark that, in contrast to the disordered phase, the mean-field steady state for 2 > 2 does depend on the value of κ, with the mean-field bosonic order parameter X diverging as κ → 0. Interestingly, there is precedent for the occurrence of PT -symmetry breaking transitions in non-Hermitian dynamics of classical spins [67]. The non-analyticity of the steady state in this zerodimensional system is not related to the more common thermodynamic limit in which the system size diverges. Rather, as was clarified in Ref. [46], the thermodynamic limit here corresponds to κ → 0 + , or, equivalently, diverging excitation numbers [54,63,64]. While a finite value of κ > 0 limits the population of a single mode, a system without dissipation may exhibit a phase transition through a diverging population of certain modes. As such, for finite κ, we expect the system to exhibit finitesize crossover behavior instead of a true non-analytic phase transition. We simplify matters by exclusively focusing on the κ → 0 + limit; however, studying the effects of finite-size fluctuations [63] due to a small, but finite κ would be an interesting subject for future investigations. In this arrangement the system is organized as a semi-infinite one-dimensional lattice with two parallel excitation pathways. The lattice sites are indexed by n = 0, 1, 2, ..., which plays the role of a "spatial" variable, and by s = +1 (blue) or s = −1 (red), indicating the branch. States that are symmetric under the transformation C are characterized by an equal distribution function on both branches. We map the rate equation (31) to the problem of hopping on this lattice, schematically depicted here for a generic site by the arrows emanating from the site.
Bearing this in mind, we proceed on to study the κ → 0 limit in the following section by mapping the quantum master equation on to a classical rate equation.

IV. RATE EQUATION
We now proceed to derive a rate equation governing the long-time behavior of the Lindblad equation Eq. (3) in the rotating frame. This analysis is restricted to the regime below the critical point, so we only consider < c in this section.

A. Mapping to Classical Rate Equation
Recall from the exact diagonalization ofĤ rf that, below the critical driving strength, the Hamiltonian which governs the coherent part of the evolution, has a discrete spectrum. We label the associated dressed energy levels by the single quantum numbers µ = r √ m, ν = s √ n with r, s = ± and m, n = 0, 1, 2.... We work in units where g = 1. The eigenenergy from Eq. (7) is then given by The spacing between two levels is finite and reads with equality if and only if ν = µ. The eigenstates |µ , |ν are known in closed form and given in Appendix A. It is helpful to visualize the quantum states as forming a lattice in Fock space, see Fig. 2.
Working in the energy basis simplifies the quantum master equation to the form Phase Symmetry Photon Field Atomic Bloch Vector Linear Stability stable   TABLE II. Summary of semi-classical steady-state solutions for the photon field and atomic Bloch vector from Eqs. (17). For a fixed value of , we find two phases separated by a critical point at 2 = 2 . The associated order parameter is Re â = − g 2κ σy .
where ρ µν = µ|ρ|ν and In the following, we fix 0 < < 1 to some value close to unity, below the transition point at 2 c = 1. Then we find the steady state by taking κ → 0 + for fixed by considering the limit Here ρ(t) obeys Eq. (3).
To lowest order in κ, the steady-state solution is given by a diagonal density matrix. Indeed, Eq. (30) reads (E µ − E ν )ρ µν = 0 in this case. Hence, we find the first nontrivial correction due to κ to be given by the classical master equation for the diagonal entries ρ µ ≡ ρ µµ of the density matrix. The transition rates read Note that the diagonal elements are fixed by the sum rule, which ensures the master equation to be trace preserving, Starting from Eq. (33), it is then evident that we need to solve the linear algebraic problem Up to an overall equilibration time-scale defined by κ, the transition-rate matrix Γ only depends on . Consequently, the critical behavior of the steady state in the thermodynamic limit is a function of the control parameter .
Using the rate equation, it is easy to see why κ → 0 + controls the number of excitations. The rate equation description only applies to those levels which are well separated as compared to the decay rate κ. The spacing between adjacent levels decreases with increasing n, and for a given maximal n ≤ N max , the smallest spacing is This must be much larger than κ in order to justify throwing away the off-diagonals in the density matrix. This provides us with the unitless figure of merit We interpret this N max as the effective system size, since our description will encounter finite-size quantum fluctuations (due to the reappearance of the quantum coherence) when the typical quantum numbers n are of order N max . Thus, the thermodynamic limit is taken by first fixing 2 < 1 so the numerator is finite, and then taking κ → 0 so that N max is effectively unconstrained. This is precisely the limit described above and constitutes the parameter range considered in the following sections.

B. Analysis of Transition Rates
The aim of this section is to understand those properties of the transition rates Γ νµ that govern the criticality of the steady-state solution. In particular, we perform an asymptotic expansion of the rates for large ν + µ, which gives a very accurate approximation to the full expressions Γ νµ , and, at the same time, allows us to solve the rate equation.
Using the exact expression for the eigenstates presented in Appendix A, we find that the matrix elements take the form with ν = s √ n, µ = r √ m, u = cosh η, v = sinh η, and As → 1, we have η → ∞, leading to a divergence of the Bogoliubov coefficients u, v ∼ e η /2 at the critical point.
In contrast, the matrix elements n, s|â|m, r remain finite at the transition. Therefore, we will henceforth set = 1 in computing the matrix elements and only take into account the leading divergences of u and v.
The transition rates respect the Z 2 symmetry corresponding to C (which relates the two excitation branches ± √ n to each other). Thus, Furthermore, at large excitation number n, the rate for "interbranch" transitions [sign(ν) = −sign(µ))] decays as ∼ 1/n, see Fig. 3. We can therefore restrict our attention to only those processes which induce "intrabranch" transitions [sign(ν) = sign(µ)]. Within this well-justified approximation, the steady-state density matrix will also respect the Z 2 symmetry and can be written as Recall that n = ν 2 ⇔ ν = ± √ n, so that the steady state contains an equal mixture of both positive and negative energy states.
Characterizing the steady state amounts to determining the population function ρ(n). We set r = s = +1 and write in a slight abuse of notation. The behavior of the steadystate solution ρ(n) is completely determined by the control parameter . From Eq. (36), we find and so the steady-state population ρ(m) is a rightnullvector of the transition rate matrix Γ = (Γ nm ). We now discuss the behavior of the rates close to criticality. The leading divergence for η → ∞ follows from In the second equality, we have utilized the fact that the matrix elements ofâ in the eigenbasis |n, + are purely real. We introduce the symmetric and anti-symmetric matrix elements Up to an overall prefactor, which can be absorbed into a rescaling of time, we can then write the transition rate matrix as with asymmetry parameter A nonzero value of p implies Γ nm = Γ mn . If p = 0, then Γ is symmetric and the constant state (1, 1, . . . , 1) is a nullvector. (This follows from the trace preserving property discussed in Eq. (35), which implies n Γ nm = 0.) Such a state, however, is non-normalizable in an infinite Hilbert space, where n is allowed to be arbitrarily large. Thus, the small asymmetry induced by p > 0 is crucial in determining the width of the steady state in the thermodynamics limit, and we expect the width to diverge as p → 0.
The matrix elements S nm and B nm can be computed in closed form, see Appendix B. However, they are too complicated to be useful in practice. Therefore, in order to simplify the rate equation, we consider their asymptotic behavior for n + m → ∞ while keeping n − m fixed. We then find This scaling turns out to be a good approximation even for moderate values of n + m. In the following, we sketch how these forms are obtained analytically, and corroborate them numerically using the full expressions. The analytical derivation of Eq. (49) is presented in detail in Appendix B. The central steps are the following. Define the matrix element which can be expressed in terms of associated Laguerre polynomials. Now set m = n + δ and consider the limit n → ∞ while keeping δ ∈ Z fixed. Using the asymptotic behavior of Laguerre polynomials, we find, for δ = 0, that with J k (y) the Bessel function of the first kind. This implies that the symmetric and anti-symmetric matrix elements behave as For δ → ∞, the Bessel functions can be expanded in terms of Airy functions, giving a simple power-law behavior in δ as S n,n+δ ∼ C 1 √ n|δ| −5/3 and B n,n+δ ∼ C 2 |δ| 1/3 with two constants C 1,2 . This is actually a very good approximation even for moderate values of δ of order unity. We have C 2 ≈ 1. After another rescaling of time, we eventually arrive at Eq. (49).
We now numerically test the simplifications that lead to Eq. (49) against the exact matrix elements Γ νµ . In Fig. 3, we show the exact result for the symmetric matrix elements for intrabranch and interbranch transitons, S nm ∝ (Γ νµ + Γ µν ) r=s=1 and S (inter) nm ∝ (Γ νµ + Γ µν ) r=−s=1 , respectively. This justifies our approximation to neglect the interbranch transition rates, which, as can be seen clearly in the plot, are several orders of magnitude smaller than the intrabranch rates.
In Fig. 4, we plot the exact and asymptotic forms of the intrabranch rates S nm and B nm for fixed values of n + m. We verify a clear power-law decay of S nm , which is tightly localized around the main diagonal. Obviously, the asymptotic formula becomes worse as |n − m| approaches n+m 2 , but is still surprisingly accurate. Furthermore, in this region, the magnitude |S nm | ∼ 10 −3 . This needs to be compared to values of order unity on the main diagonal. The outer regions are thus unimportant in comparison to transitions which occur near the main diagonal.
Having established the behavior of the symmetric part of the matrix, we now turn to the anti-symmetric part B nm . By construction, we have B nm ∼ sign(m − n). In Fig. 4, we plot slices of B nm for constant n + m as a function of |n − m|. Apart from the anti-symmetry, we see that all slices have the same universal behavior near the main diagonal, given by the power law from Eq. (49). Again, deviations increase as we move further away from the main diagonal, but in this region the function |S nm | 2 is small and thus these deviations are not important for the solution of the master equation.

C. Solving the Rate Equation
In the previous section, we have established, both analytically and numerically, that the transition-rate matrix Nonequilibrium steady state below criticality, parametrized through the distribution ρ(n) in Eq. (42). We numerically obtain ρ(n) by solving Eq. (44), using the approximations (54) for the transition rates. The solid line corresponds to a truncated Fock space with n ≤ Nmax = 4000, while the dashed line is for Nmax = 5000. Both solutions agree well for small values of n. The straight dotted lines indicate that the initial portion of the distribution is wellapproximated by an exponential distribution, independent of Nmax. The universal part of ρ(n) only depends on the asymmetry parameter p = 2 √ 1 − 2 , which controls the distance from the critical point. As p → 0, the distribution becomes broader, occupying more of the highly excited Fock states.
can be very accurately approximated by (for n = m) with the asymmetry parameter p = 2e −2η = 2 √ 1 − 2 controlling the approach to the critical point. [The diagonal elements Γ nn are fixed by the sum rule Eq. (35).] We numerically solve Eq. (44) using this form of Γ nm for small values of p. Specifically, we find the steady-state population ρ(n) as the right-eigenvector of the matrix Γ with eigenvalue zero. By virtue of the approximations made in the transition-rate matrix, we are able to compute the matrix elements for large values of n ≤ N max . The result of this analysis is shown in Fig. 5.
As we decrease p → 0, the distribution ρ(n) becomes wider. The decay in n is roughly exponential, which we will later support by a simple qualitative argument. The steady-state distribution can be characterized by its spread in terms of the mean quantum number n = n nρ(n).
Note that, since we work in the dressed basis, this is not the average photon number, but a closely related quantity, as we detail below. If ρ(n) decays exponentially over a length scale ξ, thenn ∼ ξ. (This may also be taken as a definition of ξ.) As p → 0, the distribution becomes wider andn diverges. In Fig. 6 of p for two different cutoffs N max . For sufficiently large system size, we observen ∼ 1/p. To summarize, we have established that, below the critical drive strength, the steady state of the master equation (3) in the thermodynamic limit of κ → 0 + can be approximated by the steady state of the classical rate equation (33). The latter problem can be solved efficiently by replacing the exact transition rate matrix with its asymptotic form (54), which features quasi-local hopping with asymmetry controlled by the distance to the critical point. The solution obtained in this manner shows universal behavior for large population that is independent of the truncation of the Fock space.

V. STEADY-STATE PROPERTIES
In this section, we discuss general properties of the steady-state solution found in the previous section. First, we propose a simple qualitative picture of the behavior of the steady state near the critical point. We then discuss the scaling of observables at the critical point.

A. Qualitative Model
In order to derive a simple qualitative picture of the physical processes at the transition, we consider a simplified model with asymmetric nearest-neighbor hopping on the semi-infinite lattice comprised of integers n = 0, 1, 2, . . . . This is meant to model the population transfer in the dressed Fock space induced by photon loss events. The corresponding classical rate equation, for n ≥ 1, reads This represents a biased diffusion process on the dressed Fock space lattice depicted in Fig. 2: hopping to the left (decreasing n) is favored over hopping to the right (increasing n) via the asymmetry parameter p > 0. This leads to a significant population of high-energy states in the steady state. We choose reflecting boundary conditions at the origin, so that only the first and third term in Eq. (56) remain for n = 0 (with appropriate adjustment of the onsite term to maintain conservation of probability). Equation (56) can be solved through the ansatz with "momentum" q. We obtain the dispersion relation This admits two steady-state solutions with γ = 0. The first one, with q = 0, is not permissible on the semiinfinite system as it is not normalizable. The second solution consists in choosing q such that cosh q−p sinh q− 1 = 0. For small p ∼ 0, we have q ∼ 2p > 0, in which case the solution decays exponentially according to This model qualitatively reproduces two important features of the full solution to the classical master equation The toy model (56) satisfies Kolmogorov's criterion for reversible dynamics [68]. That is, for any closed cycle of states in configuration space {s 1 , s 2 , ..., s M }, the product of transition rate matrix elements is the same in either sense this loop is traversed, so that A consequence of this is that detailed balance is satisfied. Hence, the steady-state distribution of this toy model represent an effective equilibrium distribution ρ n ∼ e −En/T , with E n = n and effective temperature T eff ∼ q ∼ 2p.
Given that this toy model appears to accurately reproduce qualitative features of the solution to the actual rate equation (33), we may also question whether the solution to our problem is truly out-of-equilibrium or whether it too features emergent equilibrium behavior [50][51][52][53]69]. On the one hand, one can see that the transition rates of Eq. (54) do not satisfy the Kolmogorov criterion. This may be checked for the simple three-cycle of adjacent neighbors, for which we find the asymmetry (61) Thus, Kolmogorov's criterion is violated and detailed balance is not satisfied for every finite p > 0.
Nevertheless, one can verify that the steady-state distribution approximately satisfies detailed balance for p → 0, at least in the region where its support is largest. Thus, it can be approximated by a "thermal" distribution with effective temperature T eff ∼n ∼ 1/p in the region of small n (see Fig. 5 and Fig. 6). As a result of this, the critical point describing the photon blockade breakdown may fall into a classical universality class known from equilibrium phase transitions [50,51]. We leave exploring this exciting direction for future work.

B. Scaling of Observables
Using the steady-state density matrix computed in Section IV, we now compute the scaling of several experimentally accessible observables in the regime < c . We begin by showing that the order parameter â is exactly zero in the thermodynamic limit κ → 0 + . This is due to the steady state respecting the Z 2 symmetry generated by C (which swaps the two excitation pathways ± with each other). For this, it is important that the eigenstates |ν are related to the "Fock states" |n, s [Eq.
Consequently, the physical photon operator obeys (see App. B) Thus, if ρ ss is an equal admixture of the two branches s = ± according to Eq. (42), then the two contributions cancel, implying â = 0. We next consider the photon number, which characterizes the fluctuations of the order parameter in the disordered phase. From Eq. (B12), we obtain the expectation value in the eigenstate |ν as which is independent of the sign s. This implies The exponent which governs the divergence of the photon number upon approaching the critical point is known as the "photon flux exponent." We find that this exponent is 1, which is consistent with the "classical" equilibrium Dicke model [69]. In this case, this is not entirely an obvious result as it arises from a combination of two distinct contributions: a power of 1 2 comes from the divergent number of photons in each individual eigenstate |ν and another power of 1 2 comes from the divergent width of the distribution ρ(n) over the eigenstates. These contributions combine to produce the overall exponent of one.
Finally, let us comment on the behavior of the atomic Bloch vector. In the steady state, we have the exact relation This implies iκ â = E + g σ − . Since â = 0 for < c , we arrive at The underlying physical picture is that of the Bloch vector coherently canceling out the driving field. This relation is satisfied by our mean-field theory, see Eq. (20). We can also check that the explicit steady-state solution to the rate equation (33) satisfies this condition. For this, we use the results of App. B to deduce ν|σ x |ν = − and ν|σ y |ν = 0 for all the eigenstates. Since the steady state is diagonal in the energy eigenbasis, we conclude that σ y = 0 (70) for the steady state as well, which eventually yields Eq. (68). Computing the expectation value σ z is slightly more subtle. We have ν|σ z |ν = 0 for every ν = 0 . Thus, the only possible contribution to the steady-state expectation value comes from the dressed vacuum, which has ν = 0 and produces Thus, the steady-state expectation value ofσ z depends on the population of the dressed vacuum state according to The value of ρ(0) is not universal and depends on the particular approximation scheme used to compute the steady state. However, if we model the steady-state decay in n as exponential over a length ξ ∼n (as confirmed by Fig. 6), we can estimate that, as 2 → 1, the population of the state |0 behaves as For our steady-state solution from the rate equation we have 1/n ∼ p ∼ √ 1 − 2 and so we arrive at This is different from the mean-field theory prediction σ z MF ∼ √ 1 − 2 . Although obtained by directly numerically solving the rate equation, the scalingn ∝ 1/p is supported on more general grounds by the solution to the qualitative tight-binding model, the solution of which behaves as Eq. (59).
This departure from mean-field theory is even more evident if we compute the length of the Bloch vector, 2 = | ˆ σ | 2 , which also serves as a measure of the purity of the spin (upon tracing out the photon). Recall that in mean-field theory this is a constant of motion and so does not scale with the control parameter . Using the above results, we find instead that to leading order Thus, the rate equation predicts that, as the critical point is approached, the length of the atomic Bloch vector approaches unity. This implies that, as we approach the critical point, the reduced density matrix for the atom tends towards a pure state. Knowing this may be useful for future investigations, as it motivates a natural simplifying ansatz for studying the critical behavior.

VI. EXPERIMENTAL PLATFORMS AND OUTLOOK
In this paper, we have provided a comprehensive analysis of the critical behavior of the driven-dissipative Jaynes-Cummings model upon approaching the breakdown of photon blockade. In the following, we discuss potential experimental platforms and point out future directions for extending our considerations.
Several experimental platforms qualify to observe this critical behavior. The necessary experimental ingredients are (i) a qubit degree of freedom with long coherence time, (ii) a boson with long coherence time that can be coherently driven, and (iii) the ability to realize strong coupling between the qubit and boson.
Arguably the most direct route to realize the model, as originally envisioned in Ref. [61], is by using an atom, a molecule, or a quantum dot coupled to a singlemode high-quality optical cavity. Equivalently, one could couple a superconducting qubit to a microwave cavity. Indeed, photon blockade has been demonstrated with atoms [25], quantum dots [24], and superconducting qubits [27,60]. However, a number of obstacles remain in order to observe the critical breakdown of photon blockade in these systems, including limitations on the lifetimes of both cavity and qubit, presence of additional cavity modes, and deviations from the two-level approximation for the qubit.
Another promising candidate is to realize the system using an internal-state qubit of a trapped atomic or molecular ion coupled to its motion. If multiple ions are trapped, the qubit can be coupled to a phonon mode of the corresponding Wigner crystal. Trapped-ionbased platforms offer a number of technical advantages, making such platforms well-suited for quantum simulation [33,34,70], and the analogue of photon-blockade has been demonstrated in trapped ions [35]. Since the boson in this system is the phonon mode of the ion or of the ion crystal, the diverging boson number associated with the photon blockade breakdown may actually entail loss of the the ion from the trap or destruction of the ion crystal, though it may be possible to still get to large phonon occupation number before this occurs.
A third type of system which may be promising centers around replacing the photon (i.e. the bosonic subsystem) with motional quanta of a nanomechanical resonator [21,22,39,40,71]. In order to realize the qubit, one may envision strong coupling of an atom to the nanomechanical resonator, as in Ref. [71]; however, a particularly promising emerging qubit candidate is the solid-state spin qubit [39]. In particular, NV-centers in diamond can provide qubits with extremely long coherence times, and coupling the qubit to the nanomechanical motion has been demonstrated [21]. Of special note is the coupling of the spin qubit to the motion of a levitated micromagnet [22,40], which can potentially reach the required strong-coupling regimes, while remaining wellisolated from the environment.
Another, more abstract proposal might capitalize on the analogy between the Jaynes-Cummings model and Landau levels in graphene [66,72] in order to realize the critical theory in a purely solid state setting. In this case, the bosonic degree of freedom is realized by the cyclotron motion of the electrons, while the atomic pseudo-spin degree of freedom is realized by the electron's Bloch-band degree of freedom.
In order to compare our findings to potential experiments, it is crucial to understand how deviations from the model affect the physics. In the following, we discuss several potential extensions of the model that comprise important immediate directions for future research.
First, the critical point identified in Ref. [46] and studied here focuses on the specific point in parameter space where the bare frequencies of atom, cavity, and drive are all resonant. This is motivated by the observation that the Z 2 particle-hole symmetry induced by C is no longer a symmetry of the model when the cavity frequency is changed, while keeping the drive resonant with the atom. The rotating-frame Hamiltonian is then modified to include the cavity detuning δ ph : While the last two terms, which are discussed in this work, are odd under C , the first term is even. Thus, the Z 2 symmetry is explicitly broken by the detuning δ ph = 0, which inhibits critical fluctuations. The photon blockade breakdown transition is first-order in this case, see Refs. [45,46,[49][50][51]. Quite remarkably, a detuning δ a = 0 of the atomic frequency with respect to the drive (while keeping the cavity resonant with the drive) does not necessarily imply a first-order transition. As was shown in Ref. [43], adding a finite atomic detuning still preserves the exact solubility of the driven model. In such a case, the dressed spectrum is given by We see that, even though the C -invariance seems to be spoiled by the detuning, there may still be a critical point due to the collapse of the eigenvalue spacing, which vanishes as ∼ √ 1 − 2 . In our analysis, although we accounted for photon loss, we neglected spontaneous emission from the atomic degree of freedom. The inclusion of atomic decay is expected to drastically modify the critical behavior [46]. To understand this, recall that photon annihilation operatorâ essentially has no interbranch matrix elements. However, this is not true for the atomic operatorσ − . For instance, in the absence of driving, the amplitude for switching branches due to atomic emission is n − 1, −|σ − |n, + = 1 2 , irrespective of the value of n. Contrast this to the interbranch transition amplitude for photonic emission given by n − 1, −|â|n, + = 1 2 ( √ n − √ n − 1) ∼ n −1/2 , Whereas this goes to zero for large quantum numbers n, thereby effectively decoupling the two branches when the dissipation is due to photon loss, the inclusion of atomic dissipation may qualitatively alter the steady-state properties.
Two more important effects to be considered, especially if a connection to experiment is to be made, are counter-rotating terms and the presence of multiple atomic levels. The former has been discussed to some degree in Refs. [62,66], and there is reason to believe it is not a fundamental impediment towards realizing this critical point in experiment. The latter may present a more important issue since many proposed qubit systems do not have the sufficiently strong non-linearity needed to project out higher internal excitations. This is particularly troublesome in superconducting qubits [26][27][28][29][30][31], where higher lying levels can present complications to the effective Jaynes-Cummings picture.
Even within the confines of the present model, several intriguing and experimentally relevant open problems remain. For instance, how does the transition look when approached from the "other side", i.e. for 2 > 1. Based on our semiclassical calculations, as well as our argument based on the rate equation, we should not expect this part of the phase diagram to be well defined if we insist on taking the κ → 0 + limit. Nevertheless, we can still study the behavior of the steady-state for finite κ. In this case, we expect the critical point to evolve into a smooth crossover. Studying this in a quantitative manner may be achieved in terms of a functional inte-gral description [50,51,69,73,74], by building upon the saddle-point semiclassical solutions, now including fluctuations due to finite κ. This has the additional advantage in that it allows us to study dynamical correlation functions whereas in this work we are limited to only studying static correlation functions. This may be particularly relevant given that the steady state appears to admit an emergent equilibrium description [50,51], despite the microscopic violation of detailed balance. Studying the dynamical correlations of the system can then allow for determining whether the system obeys a generalized fluctuation-dissipation relation, as is often the case when systems display emergent equilibrium [50,63,68,69,74].
This approach may also naturally lend itself towards a generalization to include a many-body version of this transition. Specifically, we might imagine a Bose-Hubbard-like system where, instead of an onsite Hubbard non-linearity of the formâ †â †ââ , we consider an onsite Jaynes-Cummings-type non-linearity. Such a model has been considered in Refs. [75][76][77][78][79], but not in the context of the breakdown of photon blockade. In principle, this could dramatically alter the ground-state phase diagram. To see why, we might imagine a quantum fluctuation in the order parameter of a neighboring site which is sufficiently strong as to induce a breakdown of the photon blockade. This would result in a destruction of the Mott state locally. In principle, these tunneling events may then proliferate and ultimately destabilize the Mott lobe entirely. We follow Ref. [43] and obtain the eigenstates and eigenvalues of the Hermitian operator (rescaled by g) the unitless control parameter. We will henceforth, without loss of generality, consider > 0. We use the convention for the displacement and squeezing operators (η ∈ R) which act on the annihilation operator as with u = cosh(η) and v = sinh(η). Let us consider the eigenvalue problem We will use the convention that ν labels the eigenstates and λ ν is the eigenvalue for eigenstate ν. First, with the knowledge of the exact solution, we perform a squeezing transformation of the form |ν ≡ S(η)|n, s . (A7) We will later determine the exact squeezing parameter η. For the moment, the notation |n, s is merely suggestive, but, as we will see, this is consistent with the quantum number labeling scheme used throughout the paper. This transformation acts onĤ to produce where we have defined the new two-by-two matrix We observe that hence, the matrixτ − is singular when This only has a solution if 0 < 2 /(2 − 2 ) < 1 ⇒ 2 < 1. If this condition is met, then we can choose η such that This leads to the following matrix representation forτ − : Let χ R , χ † L be the normalized left-and right-eigenvectors with zero eigenvalue respectively, so that Note that, by virtue ofσ xτ+σx =τ − , these are related by χ L =σ x χ R . We have the specific representation of χ R as We begin by determining the ground-state (vacuum) |0 , as the state which is annihilated by bothτ − andâ. The eigenvalue is zero and the ket is where |vac is the photon vacuum state. Next we determine the excitations above the vacuum. These are organized into two-dimensional sub-spaces of the form |n, s = Aχ R ⊗ |ψ R + Bχ L ⊗ |ψ L . (A17) We insert this ansatz and then act on the left with χ † R and χ † L . Using that these are zero left-eigenvectors ofτ + andτ − respectively, we find (A18) These may be decoupled to obtain the equation for the "left" ket as and the equation for the "right" ket as (A20) We can solve this by a displaced Fock state such that the displacement and the eigenvalue, with n = 0, 1, 2, 3.... This has two branches of solution, The kets |n, s are then determined by using equation (A19) and normalizing.
We thus obtain the final result for the spectrum as with "ladder" states where s = ± and n = 1, 2, ..... The squeezing parameter is obtained as and α ν , λ ν are and, as usual, u = cosh η, v = sinh η.

Appendix B: Matrix Elements
In this Appendix, we compute the matrix elements of various operators in the driven-eigenbasis computed in Sec. A.
(B2) Recalling that v = sinh η, we obtain the scaling upon approaching the critical point of For the expectation values of the spin operators, we have a relatively simple calculation of For the states with ν = 0, the result is a bit more complicated. We first perform the squeezing transformation to go from the eigenstates |ν to the ladder states |n, s to obtain ν|â †â |ν = n, s| uâ + vâ † uâ + vâ † |n, s , (B6) ν|ˆ σ|ν = n, s|ˆ σ|n, s .
The field expectation value is the simplest, and can be computed to be n, s|â|n, s = 3 2 α ν .
We combine these with the behavior of the Bogoliubov coefficients u = cosh η and v = sinh η to obtain The results for both the ν = 0 and ν = 0 states are summarized in Tab. III. Next, we compute the transition elements between eigenstates via single-photon emission events. Let ν and µ be the two eigenvalues connected by the transition. Then we need the off-diagonal elements for the transition µ → ν. These are expressed in terms of the ladder states of Eq. (A26) as ν|â|µ = n, s|uâ + vâ † |m, r , where ν = s √ n is the final quantum number and µ = r √ m is the initial quantum number. To compute ν|â|µ , define the ladder-state matrix elements ofâ as A n,s|m,r = n, s|â|m, r .
This quantity is regular at the transition point and so we set = 1. We have Due to orthogonality this simplifies to with In order to evaluate this expression, we employ the generating function (B20) The derivative with respect to α * may be performed to obtain 0|â p D(α)â †(q) |0 = e |α| 2 /2 ∂ ∂α q α p e −α * α . (B21) As the derivative treats α and α * independently, we can now rescale α by α * via α = r/α * which then allows us to write with associated Laguerre polynomials L (k) n (z). The second equality follows by application of the Rodrigues formula for associated Laguerre polynomials. We can now take α to be real to obtain the result We thus obtain the manifestly real expression Note that the matrix elements satisfy and so we only need to consider two out of four possible combinations of the signs r and s: one same-sign combination and one opposite-sign combination. We are interested in the asymptotic behavior of A n,s|m,r for large m, n → ∞ as a function of the difference m − n. We consider the case of equal sign and set r = s = 1. Define the relative coordinate so that α = R for our sign choice. We have We write m = n + δ with integer δ = 0 so that To study the asymptotic behavior of these expressions in the limit n ∼ ∞ while keeping δ fixed, we use We then have The asymptotic behavior of the Laguerre polynomials L (k) N for x > 0, N → ∞ and k fixed is determined by with J k (y) the Bessel function. Note that Hence and so Similarly, we obtain We conclude that the asymptotic behavior of Eq. (B27) for for n → ∞ with m = n + δ and δ ∈ Z is given by where The first line follows from Eq. (B24) in the limit α → 0, and the second line follows from J k−1 (y) + J k+1 (y) = 2k y J k (y) with k = δ − 1. We next justify neglecting the interbranch transitions by showing that they decay exponentially in n for large n, m. For this purpose, we choose s = 1 and r = −1 in Eq. (B24). Setting n = m (corresponding to δ = 0), we have α = 2 √ n and arrive at A n,1|n,−1 ∼ √ ne −2n
Note that, while f (δ) is oscillatory for real δ < 0, we only evaluate the function for integer arguments δ ∈ Z\{0}, where these oscillations are not visible. For y > 0 we use the integral representation of the Bessel functions to write For large δ → ∞, the integral will be dominated by the region of small angles. We expand sin(θ) − θ = − 1 6 θ 3 + O(θ 5 ) and introduce the variable x 3 = y 2 y 3 and extend the integration boundaries to infinity. We then arrive at In fact, these formulas are excellent approximations even for small y ∼ O(1). We further have a(y) s(y) where the prefactor is close to unity.

Appendix C: Dynamic Stability of Mean-Field Equations
In this section, we investigate the stability of the set of mean-field equations (17). We write a = â , σ = σ − , σ 3 = σ z .