Quantum features of entropy production in driven-dissipative transitions

The physics of driven-dissipative transitions is currently a topic of great interest, particularly in quantum optical systems. These transitions occur in systems kept out of equilibrium and are therefore characterized by a finite entropy production rate. However, very little is known about how the entropy production behaves around criticality and all of it is restricted to classical systems. Using quantum phase-space methods, we put forth a framework that allows for the complete characterization of the entropy production in driven-dissipative transitions. Our framework is tailored specifically to describe photon loss dissipation, which is effectively a zero temperature process for which the standard theory of entropy production breaks down. As an application, we study the open Dicke and Kerr models, which present continuous and discontinuous transitions, respectively. We find that the entropy production naturally splits into two contributions. One matches the behavior observed in classical systems. The other diverges at the critical point.


I. INTRODUCTION
The entropy of an open system is not conserved in time, but instead evolves according to where 0 is the irreversible entropy production rate and is the entropy flow rate from the system to the environment. Thermal equilibrium is characterized by dS/dt = = = 0. However, if the system is connected to multiple sources, it may instead reach a nonequilibrium steady state (NESS) where dS/dt = 0 but = 0. NESSs are therefore characterized by the continuous production of entropy, which continuously flows to the environments.
In certain systems a NESS can also undergo a phase transition. These so-called dissipative transitions [1][2][3] represent the open-system analog of quantum phase transitions. Similarly to the latter, they are characterized by an order parameter and may be either continuous or discontinuous [4][5][6]. They are also associated with the closing of a gap, although the gap in question is not of a Hamiltonian, but of the Liouvillian generating the open dynamics [7,8]. The features emerging from the competition between dissipation and quantum fluctuations have led to a burst of interest in these systems in the last few years , including several experimental realizations [39][40][41][42][43][44][45].
Given that the fundamental quantity characterizing the NESS is the entropy production rate , it becomes natural to ask how behaves as one crosses such a transition, such as, i.e., what are its critical exponents? Is it analytic? Does it diverge? Surprisingly, very little is known about this and almost all is restricted to classical systems.
In Ref. [46] the authors studied a continuous transition in a two-dimensional (2D) classical Ising model subject to two baths acting on even and odd sites. They showed that the entropy production rate was always finite, but had a kink at the critical point, with its derivative presenting a logarithmic divergence. A similar behavior was also observed in a Brownian system undergoing an order-disorder transition [47], the majority vote model [48], and a 2D Ising model subject to an oscillating field [49]. In the system of Ref. [49], the transition could also become discontinuous depending on the parameters. In this case they found that the entropy production has a discontinuity at the phase coexistence region. Similar results have been obtained in Ref. [50] for the dissipated work (a proxy for entropy production) in a synchronization transition.
All these results therefore indicate that the entropy production is finite across a dissipative transition, presenting either a kink or a discontinuity. This general behavior was recently shown by some of us to be universal for systems described by classical Pauli master equations and breaking a Z 2 symmetry [51]. An indication that it extends beyond Z 2 was given in Ref. [52] which studied a q-state Potts model.
Whether or not this general trend carries over to the quantum domain remains an open question. Two results, however, seem to indicate that it does not. The first refers to the driven-dissipative Dicke model, studied experimentally in Refs. [39,40]. In this system, the part of the entropy production stemming from quantum fluctuations was found to diverge at the critical point [45]. Second, in Ref. [53] the authors studied the irreversible work produced during a unitary quench evolution of the transverse field Ising model. Although being a different scenario, they also found a divergence in the limit of zero temperature (which is when the model becomes critical). Both results therefore indicate that FIG. 1. Typical driven-dissipative scenario portraying an optical cavity with a nonlinear medium subject to an external pump E and photon losses occurring at a rate κ. quantum fluctuations may lead to divergences of the entropy production in the quantum regime. Whether these divergences are universal, and what minimal ingredients they require, remains a fundamental open question in the field.
The reason why this issue has so far not been properly addressed is actually technical: most models explored so far fall under the category of a driven-dissipative process, where dissipation stems from the loss of photons in an optical cavity [54] (see Fig. 1). The problem is that photon losses are modeled effectively as a zero-temperature bath, for which the standard theory of entropy production yields unphysical results (it is infinite regardless of the state or the process) [55,56].
This "zero-temperature catastrophe" [57,58] occurs because the theory relies on the existence of fluctuations which, in classical systems, seize completely as T → 0. In quantum systems, however, vacuum fluctuations remain. This was the motivation for an alternative formulation introduced by some of us in Ref. [59] and recently assessed experimentally in [45], which uses the Wigner function and its associated Shannon entropy as a starting point to formulate the entropy production problem. This has the advantage of accounting for the vacuum fluctuations, thus leading to a framework that remains useful even when T → 0. This paper builds on Ref. [59] to formulate a theory which is suited for describing driven-dissipative transitions. Since these transitions are seldom Gaussian, we use here instead the Husimi Q function and its associated Wehrl entropy [56,60]. Our focus is on defining a consistent thermodynamic limit where criticality emerges. This allow us to separate into a deterministic term, related to the external laser drive, plus a term related to quantum fluctuations. The latter is also additionally split into two terms, one related to the nontrivial unitary dynamics and the other to photon loss dissipation. We apply our results to the Dicke and Kerr models, two paradigmatic examples of dissipative transitions having a continuous and discontinuous transition respectively. In both cases, we find that unitary part of behaves exactly like in classical systems. The dissipative part, on other hand, is proportional to the variance of the order parameter and thus diverges at the critical point.

II. DRIVEN-DISSIPATIVE SYSTEMS
We consider a system described by a set of bosonic modes a i evolving according to the Lindblad master where H 0 is the Hamiltonian, E i are external pumps, and κ i are the loss rates for each mode (see Fig. 1). The second term in Eq. (2) is the typical Lindblad dissipator describing onephoton losses of a cavity. The results below hold for arbitrary times in the dynamics, although most of our interest will be in the NESS, defined as the fixed point dρ/dt = 0. We work in phase space by defining the Husimi function Q(μ,μ) = 1 π μ|ρ|μ , where |μ = i |μ i are coherent states andμ denotes complex conjugation. The master Eq. (2) is then converted into a quantum Fokker-Planck (QFP) equation [61], where U (Q) is a differential operator related to the unitary part (see Appendixes A-C for examples) and J i (Q) = κ i (μ i Q + ∂μ i Q) are irreversible quasiprobability currents associated with the photon loss dissipators. As our basic entropic quantifier, we use the Shannon entropy of Q, known as Wehrl's entropy [60], This quantity can be attributed an operational interpretation by viewing Q(μ,μ) as the probability distribution for the outcomes of a heterodyne measurement. S(Q) then quantifies the entropy of the system convoluted with the additional noise introduced by the heterodyning [62,63]. As a consequence, S(Q) S(ρ), with both converging in the semiclassical limit. We also mention that the Wehrl entropy has the unique advantage of being well defined for any quantum state, since Q 0. This is in contrast with the Wigner entropy, which can become imaginary if the Wigner function is negative. Next, we differentiate Eq. (4) with respect to time and use Eq. (3). Employing a standard procedure developed for classical systems [64], we can separate dS/dt as in Eq. (1), with an entropy flux rate given by and an entropy production rate The entropy flux is seen to be always non-negative, which is a consequence of the fact that the dissipator is at zero temperature, so that entropy cannot flow from the bath to the system, only the other way around. As for in Eq. (6), the last term is the typical dissipative contribution, related to the photon loss channels and also found in [59]. The extension to a finite temperature dissipator is straightforward and requires only a small modification of the currents J i [59].
The new feature in Eq. (6) is the first term, which is related to the unitary contribution U (Q). Unlike the von Neumann entropy, the unitary dynamics can affect the Wehrl entropy. This is due to the fact that the unitary dynamics can already lead to diffusionlike terms in the Fokker-Planck Eq. (3), as discussed, e.g., in Ref. [65].

III. THERMODYNAMIC LIMIT
The results in Eqs. (5) and (6) hold for a generic master equation of the form (2), irrespective of whether or not the system is critical. We now reach the key part of our paper, which is to specialize the previous results to the scenario of driven-dissipative critical systems. The first ingredient that is needed is the notion of a thermodynamic limit. For drivendissipative systems, criticality emerges when the pump(s) E i become sufficiently large. It is therefore convenient to parametrize E i = i √ N and define the thermodynamic limit as N → ∞, with i finite. In driven systems a i always scales proportionally to E i , so that we can also define a i = α i √ N, where the α i are finite and represent the order parameters of the system.
The parameter N, representing the thermodynamic limit, can therefore be thought of as being proportional to the number of photons in the pump which, in turn, is roughly the number of photons in the cavity. Thus, criticality in drivendissipative models occur when the number of photons becomes very large. From a theoretical point of view, however, N is to be viewed as knob allowing one to tune the model towards a critical behavior.
This combination of scalings implies that at the mean-field level (a i → a i ) the pump term E i (a † i − a i ) in (2) will be O(N ); i.e., extensive. We shall henceforth assume that the parameters in the model are such that this is also true for H 0 in Eq. (2) (see below for examples).
Introducing displaced operators δa i = a i − α i √ N, the entropy flux (5) is naturally split as The first term is extensive in N and depends solely on the mean-field values |α i |. It is thus independent of fluctuations. The second term, on the other hand, is intensive in N. In fact, it is proportional to the variance of the order parameter δa † i δa i (the susceptibility) and thus captures the contributions from quantum fluctuations.
We can also arrive at a similar splitting for the entropy production (6). Defining displaced phase-space variables Substituting in (6) then yields This is the main result in this paper. It offers a splitting of the total entropy production rate into three contributions with distinct physical interpretations. The first, ext , is extensive and depends solely on the mean-field values α i . It therefore corresponds to a fully deterministic contribution, independent of fluctuations. Comparing with Eq. (7), we see that a balance which holds irrespective of whether the system is in the NESS. Hence, this contribution does not affect the system entropy: At the mean-field level, all entropy produced flows to the environment. The second and third terms in Eq. (8) represent, respectively, the unitary and dissipative contributions to . These two terms account for the contributions to the entropy production stemming from quantum fluctuations. This becomes more evident in the NESS (dS/dt = 0), where combining Eqs. (1) and (9) leads to The two terms u and d therefore represent two sources for the quantum entropy q in Eq. (7). We also note in passing that while d 0, the same is not necessarily true for u , although this turns out to be the case in the examples treated below.

IV. KERR BISTABILITY
To illustrate how the different contributions to the entropy production in Eq. (8) behave across a dissipative transition, we now apply our formalism to two prototypical models. The first is the Kerr bistability model [11,26,54], described by Eq. (2) with a the single mode a and Hamiltonian where is the detuning and u is the nonlinearity strength. This model has a discontinuous transition.
The NESS of this model and the terms in Eq. (8) were computed using numerically exact methods. Details on the numerical calculations are provided Appendix B and the main results are shown in Fig. 2. In Figs. 2(a) and 2(b) we plot u and d for different sizes N. As can be seen, u has a discontinuity at the critical point when N → ∞. Conversely, d diverges. The critical behavior in the thermodynamic limit (N → ∞) can be better understood by performing a finite size analysis [Figs. 2(c) and 2(d)], where we plot u and d /N vs N ( / c − 1) for multiple values of N. Surprisingly, we find that the behavior of u matches exactly that of the classical entropy production in a discontinuous transition [49,51,52]. We also see from Fig. 2 that u is negligible compared to d . As a consequence, in view of Eq. (10) the dissipative contribution d will behave like the variance of the order parameter δa † δa , which diverges at the critical point. This is clearly visible in Fig. 2(d), which plots d /N.

V. DRIVEN-DISSIPATIVE DICKE MODEL
The second model we study is the driven-dissipative Dicke model [39,40]. It is described by Eq. (2) with a mode a, subject to photon loss dissipation κ, as well as a macrospin where J i are macrospin operators. This model does not need a drive E since the last term can already be interpreted as a kind of "operator valued pump" (as it is linear in a + a † ). In fact, this is precisely how this model was experimentally implemented in a cold-atom setup [39]. The model can also be pictured as purely bosonic by introducing an additional mode b according to the Holstein-Primakoff map It hence falls under the category of Eq. (2), with two modes a and b.
Since this is a two-mode model, numerically exact results are more difficult. Instead, we follow Refs. [39,40,45] and consider a Gaussianization of the model valid in the limit of N large. Details are provided in Appendix C and the results are shown in Fig. 3. Once again, the unitary part u of the entropy production [Figs. 3(a) and 3(b)] is found to behave like the mean-field predictions for classical transitions [46][47][48][49][50][51][52]. It is continuous and finite, but presents a kink (the first derivative is discontinuous) at the critical point λ c = ω 0 (κ 2 + ω 2 )/ω.
The dissipative part d , on the other hand, diverges at λ c . This was indeed already shown experimentally in Ref. [45]. In fact, the behavior of d at the vicinity of λ c is of the form as confirmed by the analysis in Fig. 3(d). Similarly to the Kerr model, u is much smaller than d so that the latter essentially coincides with 2κ δa † δa in the NESS [cf. Eq. (10)]. The divergence in (13) thus mimics the behavior of δa † δa .

VI. DISCUSSION
Understanding the behavior of the entropy production across a nonequilibrium transition is both a timely and important question, specially concerning driven-dissipative quantum models, which have found renewed interest in recent years. This paper provides a framework for computing the entropy production for the zero-temperature dissipation appearing in driven-dissipative models.
We then applied our formalism to two widely used models. In both cases we found that one contribution u behaved qualitatively similar to that of the entropy production in classical dissipative transitions. The other, d , behaved like a susceptibility, diverging at the critical point. Why u behaves in this way, remains an open question. Driven-dissipative systems have one fundamental difference when compared to classical systems. In the latter, energy input and output both take place incoherently, through the transition rates in a master equation. In driven-dissipative systems, on the other hand, the energy output is incoherent (Lindblad-like) but the input is coherent (the pump). A classical analog of this is an electrical circuit coupled to an external battery E. For instance, the entropy production in a simple RL circuit is RL = E 2 /RT [66] where R is the resistance and T is the temperature. If we consider an empty cavity with a single mode a and H 0 = 0, Eq. (6) predicts cavity = 2E 2 /κ. Notwithstanding the similarity between the two results, one must bear in mind that the RL circuit still contains incoherent energy input. Indeed, RL diverges as T → 0. The cavity, on the other hand, relies solely on vacuum fluctuations. This interplay between thermal and quantum fluctuations highlights the need for extending the present analysis to additional models of driven-dissipative transitions. In particular, it would be valuable to explore models which can be tuned between classical (e.g., for large T ) and quantum (T = 0) transitions.

APPENDIX A: PROPERTIES OF u
The entropy production rate in Eq. (6) of the main text has a term proportional to the unitary dynamics, The corresponding phase-space contribution U (Q) can be found using standard correspondence tables [61] to convert the master equation term −i[H 0 , ρ] into a corresponding differential operator for Q. The result is Normal ordering is convenient as it pushes all derivatives to the right. We now change variables to ν = μ − α √ N and expand the result in a power series in N.
This yields, to leading order, The remaining terms in the expansion are at least O(1/ √ N ) and thus vanish in the limit N → ∞. This expression may be further simplified by introducing the constants Then, since h rs = h * sr , we can write (A5) as This is the leading-order contributions of the unitary dynamics to the Fokker-Planck equation. The important point to notice is the existence of diffusive terms (proportional to the second derivative ∂ 2 ν and ∂ 2 ν ). This is a known feature of the Husimi function.
We now plug Eq. (A9) into Eq. (A1). Integrating by parts multiple times and using the fact that the Husimi function always vanishes at infinity, we find that the only surviving terms are which provides the leading-order contribution to u . In the limit N → ∞ this is the only contribution which survives.

APPENDIX B: SOLUTION OF THE KERR BISTABILITY PROBLEM
In this section we provide additional details on the solution methods used to study the entropy production in the Kerr model [Eq. (11) of the main text]. The NESS of this model can be found analytically using the generalized P function [54]. This includes all moments of the form (a † ) r a s , as well as the Wigner function [67]. While the Husimi function can in principle be found numerically from the Wigner function, we have found that this is quite numerically unstable due to the highly irregular nature of the latter. Instead, it is easier to simply find the steady-state density matrix ρ numerically using standard vectorization techniques (as done, e.g., in Ref. [26]).

Numerical procedure
The numerical calculations were performed as follows. We define the Liouvillian corresponding to the master equation (2) as The steady-state equation, is then interpreted as an eigenvalue/eigenvector equation: ρ is the eigenvector of L with eigenvalue 0. To carry out the calculation, we decompose ρ in the Fock basis, using a sufficiently large number of states n max to ensure convergence. From ρ we then compute the Husimi function and the corresponding integrals numerically using standard integration techniques. The Husimi function is obtained by constructing approximate coherent states A grid of the Husimi function Q(μ,μ) can then be built to be subsequently integrated numerically. Derivatives of Q do not need to be computed using finite differences. Instead, one may notice that, for instance, with similar expressions for other derivatives. Finally, convergence of the numerical integration can be verified by computing moments (a † ) r a s of arbitrary order from the Husimi function and comparing with the exact results of Ref. [54].

Bistable behavior
For fixed κ, U , and < 0, the NESS of Eq. (B1) presents a discontinuous transition at a certain critical value c . This transition is related to a bistable behavior of the model at the mean-field level. For finite N the steady state of (B1) is unique [54]. However, as shown recently in Ref. [26], in the limit N → ∞ the Liouvillian gap between the steady state and the first excited state closes asymptotically in the region between ± = n ± [κ 2 + ( + n ± u) 2 ] and From a numerical point of view, however, this causes no interference since all computations are done for finite N, where the NESS is unique.

Unitary contribution to the quantum Fokker-Planck equation
The unitary contribution U (Q) appearing in Eq. (4) of the main text can be obtained using standard correspondence tables [61] and reads When plugged into Eq. (A1), the terms proportional to and E vanish. The only surviving terms are Substituting μ = α √ N + ν yields a leading contribution of O(1) which, of course, is the same as that which would be obtained using Eq. (A10) with r = s = 2.

APPENDIX C: SOLUTION OF THE DRIVEN-DISSIPATIVE DICKE MODEL
Here we describe the calculations for the driven-dissipative Dicke model [Eq. (12) of the main text]. We consider only a single source of drive and dissipation (E, κ ) acting on the optical cavity mode a. The full master equation is then with Since this system involves two modes, direct solution by vectorization becomes computationally too costly. Instead, we tackle the problem using Gaussianization. The calculations are done in detail in Refs. [39,40,45]. Here we simply cite the main results and adapt the notation to our present interests.

Mean-field solution
We start by looking at the mean-field level by introducing a = α √ N, J − = βN, and J z = wN. For large N we then get which are independent of N, as expected. Angular momentum conservation also imposes w 2 + |β| 2 = 1/4, which leads to two choices, w = ± 1 2 1 − 4β 2 . At the steady state this implies that β * = β, and where λ c = 1 2 √ ω 0 ω (κ 2 + ω 2 ) is the critical interaction in the absence of any external drives. The ± sign in Eq. (C7) stems from the two choices w = ± 1 2 1 − 4β 2 respectively. The minus solution in Eq. (C7) always yields the trivial result β = 0. The plus solution, on the other hand, can be nontrivial when λ > λ c . For this reason, we henceforth focus on the solution of which yields either β = 0 or β ∈ [0, 1/2]. Moreover, this solution corresponds to w = − 1 2 1 − 4β 2 , so that the spin is pointing downwards.

Stabilization of the solution
The Gaussianization procedure above explicitly already takes the limit N → ∞. Because of this, it turns out that on order to obtain a stable steady state, it is also necessary to add a small dissipation to δb. Here we do so in the simplest way possible, as a zero-temperature dissipator. We therefore consider the evolution of the Gaussianized master equation where D[L] = LρL † − 1 2 {L † L, ρ}. The value of γ was actually determined experimentally in [45] and is more than six orders of magnitude smaller than κ. One must therefore use a nonzero value, but the value itself can be arbitrarily small. In Fig. 3 of the main text, we have used γ = 10 −3 κ simply to ensure numerical stability.

Lyapunov equation
Once Gaussianized, we can study the steady state by solving for the second moments of δa and δb. Define quadrature operators with identical definitions for δq a and δ p a . The Hamiltonian (C16) then transforms to Next define the covariance matrix (CM) Since both the Hamiltonian and the dissipator are Gaussian preserving, the dynamics of σ is closed and described by a Lyapunov equation, where and D = diag(γ , γ , κ, κ ).
The assumption that the state of the system can be Gaussianized allows us to write down the Husimi function of the NESS, which has the form Q = 1 π √ |σ + I 4 /2| exp − 1 2 r T (σ + I 4 /2) −1 r , where r = (x b , y b , x a , y a ) are the phase-space variables corresponding to the quadrature operators R in Eq. (C23) and I 4 is the identity matrix of dimension 4. All integrals appearing in Eq. (8) will then be Gaussian and can thus be trivially computed.