Gaussian trajectory approach to dissipative phase transitions: the case of quadratically driven photonic lattices

We apply the Gaussian trajectories approach to the study of the critical behavior of two-dimensional dissipative arrays of nonlinear photonic cavities, in presence of two-photon driving and in regimes of sizable loss rates. In spite of the highly mixed character of the density matrix of this system, the numerical approach is able to provide precise estimations of the steady-state expectation values, even for large lattices made of more than 100 sites. By performing a finite-size scaling of the relevant properties of the steady state, we extrapolate the behavior of the system in the thermodynamic limit and we show the emergence of a second-order dissipative phase transition, belonging to the universality class of thermal Ising model. This result indicates the occurrence of a crossover when the loss rate is increased from the weak-loss limit, in which the phase transition belongs to the universality class of the quantum Ising model

We apply the Gaussian trajectories approach to the study of the critical behavior of twodimensional dissipative arrays of nonlinear photonic cavities, in presence of two-photon driving and in regimes of sizable loss rates. In spite of the highly mixed character of the density matrix of this system, the numerical approach is able to provide precise estimations of the steady-state expectation values, even for large lattices made of more than 100 sites. By performing a finite-size scaling of the relevant properties of the steady state, we extrapolate the behavior of the system in the thermodynamic limit and we show the emergence of a second-order dissipative phase transition, belonging to the universality class of thermal Ising model. This result indicates the occurrence of a crossover when the loss rate is increased from the weak-loss limit, in which the phase transition belongs to the universality class of the quantum Ising model Dissipative phase transitions are critical phenomena emerging in the non-equilibrium steady state of open quantum systems, due to the competition between their coherent Hamiltonian dynamics and the incoherent processes [1,2]. In the last years, the possibility of realizing strongly correlated states in photonic cavity arrays [3][4][5] has stimulated a deep investigation of these phenomena, which have been discussed theoretically in photonic systems [6][7][8][9][10][11][12][13][14][15][16][17][18], lossy polariton condensates [19][20][21] and spin models [22][23][24][25][26][27][28][29][30]. From the experimental point of view, some remarkable results have been recently obtained with driven circuit quantum electrodynamics systems [31] and semiconductor microcavities [32,33], showing thus the possibility to observe dissipative phase transitions in real open quantum many-body systems.
The non-trivial interplay between the Hamiltonian evolution, driving and dissipative processes in open quantum many-body systems has given rise to several fundamental questions about the respective roles of quantum and classical fluctuations across a dissipative phase transition. In this regard, the universality classes of these critical phenomena has been matter of an intense debate [19,20,26,27,[34][35][36][37][38]. In many cases, dissipative phase transitions emerge in regimes where the dissipation rates are comparable with the typical energy scale of the Hamiltonian and their universality classes belong to those of classical thermal phase transition [27,36,39,40], such that the critical behavior can be described in term of an effective temperature emerging as a results of dissipations.
In this debate, the case of arrays of nonlinear photonic resonators in presence of two-photon -i.e., quadratic in the field -driving has recently attracted a certain attention, as these systems may undergo a dissipative phase transition which belongs to a quantum universality class in a suitable regime of parameters [18]. Quadratically driven resonators have been realized experimentally with superconducting circuits [41,42] and, in the last years, several works have considered the possibility of exploiting the Z 2 symmetry of this system in order to simulate the behavior of quantum spin systems, proposing also noise-resilient quantum codes based on this platform [43][44][45][46][47]. In the context of critical phenomena, the spontaneous breaking of the Z 2 symmetry in a lattice of coupled quadratically driven resonators has been investigated with both a mean-field [16] and a many-body approach [18], showing the emergence of a phase transition in a regime where the loss rate is small compared to the Hamiltonian energy scales, with critical exponents equal to those of the quantum transverse Ising model.
In this letter, we consider the behavior of a quadratically driven photonic lattice when the loss rate becomes comparable with the Hamiltonian parameters, showing that the universality class of the transition changes in this regime and one recovers the description in terms of a classical criticality. We perform this study using the Gaussian trajectories approach, a novel numerical technique introduced in the study of the single Kerr cavity with a one-photon pump [48]. This method turns out to be a suitable approach to investigate dissipative phase transitions, as it allows to estimate efficiently the dissipative dynamics of large two-dimensional lattices (made of more than 100 sites) and, at the same time, to include the relevant many-body correlations which arise in vicinity of the critical point. Analyzing the finite-size scaling of the steady-state properties of the system, we estimate its behavior in the thermodynamic limit of an infinitely large lattice and we highlight the emergence of a secondorder phase transition, with critical exponents equal to those of the classical Ising model.
A lattice of N coupled photonic cavities can be described by a Bose-Hubbard model, whose Hamiltonian contains a local termĥ i acting on the i-th cavity, and a non-local term modelling the photon hopping between different cavities (we set = 1): Here,â i (â † i ) is the annihilation (creation) operator acting on the i-th site and the last sum runs over the nearestneighbor pairs ij , J and z being, respectively, the hopping strength and the coordination number. The local term can be written aŝ where ∆ is the detuning between half of the two-photon driving field frequency and the resonant cavity frequency, U is the photon-photon interaction energy associated to the Kerr nonlinearity and G is the two-photon driving field amplitude.
Assuming that the dissipative processes are Markovian, the dynamics of the open quantum system is recovered in terms of the density matrixρ(t) which solves the Lindblad master Equation: (3) Here, L is the Liouvillian superoperator andΓ j,k are the jump operators which describe the coupling of the system with the external environment: in general, the system exhibits local one-and two-photon losses, which are modelled respectively with the jump operatorsΓ 1,j = √ γâ j andΓ 2,j = √ ηâ 2 j . In the limit of long time, the system evolves towards a nonequilibrium steady stateρ SS , defined by the condition Lρ SS = 0.
The Liouvillian superoperator L defined in Eq. (3) presents a Z 2 symmetry coming from its invariance under a global change of sign of the annihilation operators,â i → −â i ∀i. In the thermodynamic limit of an infinite lattice, the Z 2 symmetry can be spontaneously broken for large values of G/γ and J/γ, as indicated by a mean-field analysis of this system [16]. This leads to the emergence of a transition between a phase with a Z 2 -symmetric steady state (i.e. Tr(ρ SS iâ i ) = 0) and a coherent phase with non-zero expectation value of the Bose field (i.e. Tr(ρ SS iâ i ) = 0).
The nature of the dissipative phase transition can be investigated in terms of an approximate spin model, where two Schrödinger-cat states with opposite parity play the role of the two s = 1/2-spin states with opposite magnetization [18]. When projecting the bosonic Hamiltonian in Eq. (1) onto the basis spanned by the Schrödinger-cat states, we recover the Hamiltonian of a quantum transverse XY model, where the photon hopping term plays the role of an anisotropic spin coupling in the xy plane (the coupling can be either ferromagnetic [18] or antiferromagnetic [49], according to the sign of J) and the detuning ∆ plays the role of a transverse magnetic field in the z-direction. The quantum transverse XY model is known to present a phase transition belonging to the universality class of the quantum transverse Ising model [50].
From a computational point of view, the numerical solution of the master equation (3) becomes intractable even when considering lattices made of few sites. The corner-space renormalization method [51], which provides a wise procedure to target the relevant subset of the Hilbert space where the steady-state density matrix lives, has been used to perform a finite-size scaling analysis of the system for large values of the nonlinearity and small values of the loss rates. However, it fails in the description of regimes characterized by a large photon occupancy per cavity and a highly mixed steady state, that are the regimes we aim to study in this work.
An alternative method which can overcome these difficulties is provided by the Gaussian Trajctories approach (GTA), which has been developed in the numerical simulation of polaritonic systems [48] and applied to study the temporal coherence of a dye-microcavity photon condensate [52]. As in an exact Wave Function Monte Carlo approach [53], the GTA recovers the density matrixρ(t) of the open quantum system from the average of a set of N T pure states |ψ n (t) (usually called quantum trajectories) obtained independently by integrating a stochastic differential equation: The physical picture underlying this formalism is to consider the environment as a measurement apparatus continuously monitoring the open system, and hence the stochastic evolution of each quantum trajectory may be interpreted as the result of the different random outcomes of these measurements [54][55][56][57]. The peculiarity of the GTA approach is the Gaussian ansatz for the quantum trajectories, which makes the pure state |ψ n completely determined by its first and second central moments [58,59]. This assumption reduces notably the computational cost needed for the numerical integration of the stochastic differential equation, since the complexity of the problem within GTA scales with the square of the number N of lattice sites, rather than with its exponential.
Contrarily to the exact Wave Function Monte Carlo method, where there is no univocal choice of the stochastic differential equation for the trajectories |ψ n (t) in order to reproduce the same master equation for the density matrixρ(t), the assumption of a Gaussian ansatz for |ψ n (t) makes the accuracy of the method dependent on the choice of the stochastic unravelings used to factorize the coupling of the quantum system with the external environment. For the single site problem with quadratic driving and losses, it was shown in [60], that in the homodyne unraveling [61], exact trajectories remain nearly coherent states, a subclass of Gaussian states. Here, we will therefore use the heterodyne unraveling, a symmetrized version of homodyne detection. Heterodyne measurements measure the field quadratures with equal weights and are equivalent to two complementary homodyne detection schemes. It is a common procedure in quantum optics, where the photon losses are interfered with an out-of-resonance reference beam [61] and is mathematically closely connected to the stochastic collapse model known as Quantum State Diffusion(QSD) [62]. Interestingly, it has recently been suggested that such unraveling naturally captures macroscopic phenomena as phase transitions [63,64]. The derivation of the full stochastic differential equations for the dynamics of the quadratically driven-dissipative photonic lattice within the GTA can be found in Supplemental material [65].
Another interesting aspect of the GTA is that the Z 2 symmetry of the Liouvillian is not preserved in the equation of motion of a single quantum trajectory, but it is restored only after the average of a large set of trajectories. Therefore, in spite of the finite size of the simulated lattice, it is possible to obtain non-trivial results in the calculation of a suitable order parameter along an individual trajectory at long times. The study of the distribution of the latter expectation values over the whole set of sampled trajectories is helpful to understand the emergence of collective phases in different regimes of the physical parameters, and hence provides an important insight into the critical behavior in the thermodynamic limit [29].
We apply the GTA formalism to the study of the quadratically-driven dissipative Bose-Hubbard model in 2D square lattice of different size, with periodic boundary conditions. In Supplementary material [65], we show that the method is able to replicate the finding of quantum-critical behavior for parameters similar to [18]. Since our main objective is to investigate the dissipative phase transition in regimes where the loss rates are comparable with the Hamiltonian parameters, we here set U/γ = J/γ = 1, ∆ = −J. The last condition assures that the two-photon driving is resonant with the k = 0 mode of the single-particle energy band of the closed system, and is useful to avoid any bistable behavior. Moreover, since two-photon losses are not expected to play an important role in the emergent criticality [18], we set η = 0: this choice allows for a relevant speed-up of the numerical calculations ( [65]).
In Fig. 1, the steady-state expectation value for the relative occupation of the k = 0 mode is plotted as a function of the driving amplitude G for different lattice sizes. At large values of G/γ, we notice that n k=0 1, independently of the size of the lattice. This result indicates that the steady state of the dissipative system in this limit is a coherent state characterized by long-range order, since the local Bose fields on each cavity assume the same value. This is in agreement with the picture of the spontaneous symmetry breaking obtained the mean-field calculation [16]. At smaller values of G/γ, the relative occupation n k0 decreases when the lattice size increases, suggesting that n k=0 → 0 in the thermodynamic limit. This behavior supports the hypothesis of the presence of a disordered Z 2 -symmetric phase for small values of the driving amplitude. Given this last results, it is convenient to describe the behavior of the system across the phase transition in terms of an order parameter, which is related to the average value of the Bose fields on the different cavities. We thus define α = Im ψ| 1 N iâ i |ψ , i.e. the expectation value of the average Bose field on a single quantum trajectory (we consider its imaginary part in order to have a real-valued order parameter): this quantity has a clear physical meaning, as it corresponds to the outcome of a heterodyne unraveling in real experiments. In Fig. 2-(a-d), we plot the time evolution of α on single trajectories in a 6 × 6 lattice starting, choosing the vacuum as the initial state at t = 0, for different values of the driving amplitude. For G = 0.7γ [ Fig. 2-(a)], we notice a behavior supporting the emergence of a disordered phase, as the order parameter exhibits Gaussian fluctuations around the value α = 0. For G = 0.8γ [ Fig. 2-(b)], the situation is similar, but we can see the presence of  Fig. 2-(d)] according to the particular realization of the noise terms in the stochastic differential equation. This behavior suggest the emergence of an ordered phase with broken symmetry in this regime of parameters. However, when we consider an ensemble of many trajectories at a given G, we see that half of them will stabilize at long time around a positive value for α and the other half around the opposite value −α, retrieving thus a Z 2 -symmetric steady-state density matrix, as expected for a system of finite size. This behavior becomes clear from the calculation of the distribution p(α) of the probability to measure a given value α at long  Fig. 2-(e). Even if p(α) is an even function for all values of G, we notice that, when increasing G, the character of the distribution changes from a monomodal to a bimodal behavior.
A more quantitative analysis can be performed from the calculation of the Binder cumulant, defined as where µ m = dα α m p(α) denotes the m-th moment of the probability distribution p(α) [67,68]. This quantity has been deeply used in the finite-size scaling analysis of spin models in presence of a paramagnetic-toferromagnetic phase transition, thanks to its peculiarity of being a universal function of the ratio ξ/L, with ξ and L being respectively the correlation length and the finite size of the lattice [67]. This means that at a critical point, where the correlation length diverges in the thermodynamic limit, the Binder cumulant becomes independent of the lattice size L: therefore, the finite-size scaling analysis of U L is a useful approach to track the emergence of the criticality in our system. In Fig. 3-(a), we plot U L as a function of the driving amplitude G for lattices of different size. We find that all the different curves cross at the common point G c = 0.86γ (the common crossing is more appreciable in the inset, where we plot the difference 2/3−U L vs. G on a logarithmic scale), confirming the presence of a dissipative phase transition in our model. Some insight about the universality class of the phase transition can be obtained plotting the same data for the Binder cumulant as a function of the quantity (G − G c )L 1/ν , where ν is the critical exponent of the correlation length. From the results shown in Fig. 3-(b), we notice that setting ν = 1 makes all the curves collapse on top of each other, confirming the expected universal behavior of U L = f (ξ/L). The value ν = 1 corresponds to the critical exponent of the correlation length for a 2D classical Ising model. This latter results indicate that the dissipative phase transition of a quadratically-driven Bose-Hubbard model in the regime of large loss rates results from a classical criticality, contrarily to what happens in regimes of small loss rate where the phase transition has a quantum nature [18].
In this work, we have applied the Gaussian trajectory approach in the study of quadratically driven photonic lattices across a dissipative phase transition. This method turns out to be particularly suitable for theoretical analysis of critical phenomena in open quantum systems, as it provides accurate estimates of the physical observables even in regimes where the steady state is highly mixed.
In the case under consideration, we have been able to simulate lattices made of up to 144 sites, which has allowed to perform a precise finite-size scaling of the relevant properties of the system and to extract the critical exponent ν for the correlation length. The result ν = 1 is found, which differs from the value predicted by meanfield theory and corresponds to the classical Ising model, which reveals the ability of the method to describe the many-body correlations arising among the photons. Furthermore, in the low-loss regime ( [65]) we replicate the finding of the value ν = ν q of the quantum Ising model which reveals that the method is able to capture the entanglement leading to relevant quantum correlations. This suggests a scenario in which, similarly to the case at thermal equilibrium, a crossover between quantum and classical criticality in the non-equilibrium steady state occurs, depending on the scale of the loss rate relative to the Hamiltonian energy scale.
The possibility to use the Gaussian ansatz in the description of strongly correlated photons can open up several ways to the application of our method in quantum many-body physics with light. An interesting perspective could be to use Gaussian trajectories to investigate the emergence of collective phenomena in photonic lattices in presence of geometric frustration [49,69] or disorder [70]. Further, it is expected that the method is readily applicable to the study of the dynamics is these systems.
Our case study to the quadratically driven kerr lattice may also contribute to the general understanding of the crossover from quantum to classical critical behavior in dissipative phase transitions.
We acknowledge stimulating discussions with F. Minganti and D. Huybrechts. W.V. acknowledges financial support from the FWO project 41190, the FWO travel grant V413119N and the hospitality of LTPN.

QUANTUM PHASE TRANSITION AT SMALL LOSS RATE
In this section, we apply the GTA to the study the lattice of quadratically driven Kerr resonators in regimes where the one-photon losses is small respect to the Hamiltonian and we benchmark it with the results obtained in Ref. [2]. The purpose of this study is two-fold. At first, we want to demonstrate that the GTA is not limited to the description of classically correlated states in open quantum systems, but can be applied also in the study of dissipative phase transitions with a quantum criticality: indeed, the Gaussian ansatz may exhibit multimode entanglement, as shown also in several studies of continuous-variable quantum information [4]. Secondly, we want to show that two-photon losses do not play a relevant role in the emergence of the dissipative phase transition, considering the critical behavior of the system for η = 0.
In Fig. 1-(b), we show the results for the steady-state expectation value of the parity Π = exp iπâ †â , as obtained with the GTA, for the following set of parameters: U = 40γ, J = 20γ, ∆ = −20γ, η = 0 (i.e. the same regime of parameter considered in Ref. [2], except for the value of the two-photon loss rate). The results show the emergence of a critical point at G c 1.2γ, with the critical exponents of the quantum transverse Ising model.
Note that the parity cannot be used to study the phase transition with classical criticality in the regimes of large loss rates, that we considered in the main text. In the latter case, indeed, the strong classical fluctuations mix the even and odd eigenstates of the steady-state density matrix at all values of G, such that Π = 0 on both sides of the phase transitions.