A Equations of motion in the U\rightarrow\infty limit

# Collective phases of strongly interacting cavity photons

## Abstract

We study a coupled array of coherently driven photonic cavities, which maps onto a driven-dissipative XY spin- model with ferromagnetic couplings in the limit of strong optical nonlinearities. Using a site-decoupled mean-field approximation, we identify steady state phases with canted antiferromagnetic order, in addition to limit cycle phases, where oscillatory dynamics persist indefinitely. We also identify collective bistable phases, where the system supports two steady states among spatially uniform, antiferromagnetic, and limit cycle phases. We compare these mean-field results to exact quantum trajectories simulations for finite one-dimensional arrays. The exact results exhibit short-range antiferromagnetic order for parameters that have significant overlap with the mean-field phase diagram. In the mean-field bistable regime, the exact quantum dynamics exhibits real-time collective switching between macroscopically distinguishable states. We present a clear physical picture for this dynamics, and establish a simple relationship between the switching times and properties of the quantum Liouvillian.

Despite numerous outstanding questions, the study of quantum many-body systems in thermal equilibrium is on relatively solid ground. In particular, very general guiding principles help to categorize the possible equilibrium phases of matter, and predict in what situations they can occur Fisher (1974, 1998); Sachdev (2007). In comparison, quantum many-body systems that are far from equilibrium are less thoroughly understood, motivating a large scale effort to explore non-equilibrium dynamics experimentally, in particular using atoms, molecules, and photons Kasprzak et al. (2006); Hofferberth et al. (2007); Amo et al. (2008, 2009); Deng et al. (2010); Polkovnikov et al. (2011). At the same time, it has become clear that studying non-equilibrium physics in these systems is often more natural than studying equilibrium physics; they are, in general, intrinsically non-equilibrium. For example, thermal equilibrium is essentially never a reasonable assumption in photonic systems, where dissipation must be countered by active pumping Carusotto and Ciuti (2013). Indeed, the inadequacy of equilibrium descriptions for photonic systems has long been recognized Scully and Lamb (1967), even though close analogies to thermal systems sometimes exist DeGiorgio and Scully (1970); Torre et al. (2013); Kirton and Keeling (2013); Hafezi et al. (2015); Maghrebi and Gorshkov (2016).

Until recently, photonic systems have been restricted to a weakly interacting regime. With notable progress towards generating strong optical nonlinearities at the few-photon level, for example with atoms coupled to small-mode-volume optical devices Spillane et al. (2008); Vetsch et al. (2010); Hennessy et al. (2007); Tiecke et al. (2014); Hafezi et al. (2012); Thompson et al. (2013), Rydberg polaritons Gorshkov et al. (2011); Peyronel et al. (2012), and circuit-QED devices Houck et al. (2012); Nissen et al. (2012); Reftery et al. (2014); Barends et al. (2015), this situation is rapidly changing. The production of strongly interacting, driven and dissipative gases of photons appears to be feasible Chang et al. (2008, 2014), and affords exciting opportunities to explore the properties of open quantum systems in unique contexts, while studying the applicability of theoretical treatments designed with more weakly interacting systems in mind. For example, it is not fully understood how the steady states of these systems relate to the equilibrium states of their “closed” counterparts, or how conventional optical phenomena, such as bistability, manifests in the presence of strong optical nonlinearities and spatial degrees of freedom.

We consider an array of coupled, single-mode photonic cavities described by a driven-dissipative Bose-Hubbard model Hartmann et al. (2006, 2008); Carusotto et al. (2009); Le Boité et al. (2013, 2014), which maps onto a driven-dissipative XY spin- model in the limit of strong optical nonlinearity. We perform a comprehensive mean-field (MF) study, and identify a variety of interesting steady states including spin density waves and limit cycles, which break the discrete translational symmetry of the system. The spin density waves possess canted antiferromagnetic order for a range of drive strengths, despite the ferromagnetic nature of the spin couplings. Interestingly, the exact quantum solutions exhibit short-range antiferromagnetic correlations for parameters that have notable overlap with the MF results. The system also supports collective bistable phases, which manifest in the exact quantum dynamics as fluctuation-induced collective switching between MF-like states. We present a simple relationship between this dynamics and properties of the quantum Liouvillian.

## I Model

For a system weakly coupled to a Markovian environment, the dynamics of its density matrix is governed by a master equation , where is the Liouvillian, is the system Hamiltonian, and is a dissipator in the Lindblad form Breuer and Petruccione (2002); Carmichael (1993). Here, we consider an array of coherently coupled, nonlinear, single-mode photonic cavities driven by a spatially uniform laser field with frequency , which leak photons into the environment at a rate . In the frame rotating at the driving frequency, the Hamiltonian for the coherent drive is . The system Hamiltonian is then , where

 ^HBH =−Jd∑⟨i,j⟩^a†i^aj−μ∑i^ni+U2∑i^ni(^ni−1) (1)

is the Bose-Hubbard (BH) Hamiltonian Fisher et al. (1989). Here, () creates (annihilates) a photon in cavity , is the local number operator, and the sums run from to , the number of cavities. The notation implies a further summation over all cavities that are nearest-neighbors to cavity , the number of which is given by the coordination number . The nonlinearity of the cavities is quantified by effective two-photon interactions of strength . The first term in Eq. (1) describes the hopping of photons between cavities; is the hopping rate, and is the dimensionality of the system ( for hypercubic arrays). The laser drives the system with strength , and is detuned from the cavity resonance frequency by . The coupling of the system to the environment is described by the dissipator .

For , the system evolves into a trivial vacuum at long times, with for all . The laser drive provides a photon source, and can stabilize non-trivial steady states , which satisfy . Generally, these non-equilibrium steady states are qualitatively distinct from the equilibrium states of the closed () BH model described by , which are characterized by a superfluid order parameter that spontaneously breaks the symmetry associated with particle number conservation. This symmetry is explicitly broken by the coherent laser drive, and the driven-dissipative BH (DDBH) model conserves neither energy nor particle number. Therefore superfluidity cannot emerge in the DDBH model; this is in contrast to similar models with incoherent pumps, which can support superfluid phases Greentree et al. (2006); Angelakis et al. (2007); Hartmann and Plenio (2007); Rossini and Fazio (2007); Koch and Hur (2009); Diehl et al. (2010); de Leeuw et al. (2015). The DDBH model does, however, possess a spatial symmetry generated by discrete translations along the Bravais vectors of the cavity array, which is broken spontaneously if spatial structure develops in the steady state.

We study the steady states of the DDBH model using both a site-decoupled Gutzwiller mean-field (MF) method Rokhsar and Kotliar (1991) and exact quantum trajectory simulations of finite systems Dalibard et al. (1992); Dum et al. (1992); Plenio and Knight (1998); Daley (2014). In this MF approximation, the density matrix is decomposed as a matrix product , where is the local density matrix at cavity . Further, we restrict our study to the “hard-core” limit, where strong optical nonlinearities produce a perfect photon blockade, by taking  Carusotto et al. (2009). In this limit, the photons can be mapped onto spins by an inverse Holstein-Primakoff transformation Holstein and Primakoff (1940), resulting in an effective driven-dissipative XY spin- model Angelakis et al. (2007); Joshi et al. (2013) with symmetric, ferromagnetic spin couplings, described by the Hamiltonian , where are the Pauli matrices. We derive equations of motion for the spin components ; these are given by Eqs. (A) in the appendix. For clarity, we specialize to the case of though many of the qualitative features discussed below are valid more generally.

## Ii Mean-field phase diagram

By solving the MF equations for a variety of parameters, we find heuristically that all steady states are either spatially uniform, where all spins point in the same direction ( for all ), or have antiferromagnetic (AF) spin density wave order, where neighboring spins point in different directions, but next nearest neighbors point in the same direction ( and for all ). Because the neighboring spins are not antiparallel, this AF order is “canted.” This motivates the use of a two-sublattice ansatz, which we solve by evolving the MF equations of motion on two sites. We find a variety of interesting steady state phases, shown by the colored regions in Fig. 1. In the blue region, the system exhibits bistability between spatially uniform darker () and brighter () steady states. In the red region, there is a unique steady state with AF order. In the green region, the system is bistable between and steady states. All phase boundaries in Fig. 1 correspond to continuous transitions except those at the threshold of bistability (dark lines), where the additional steady state appears discontinuously.

The bistability in this system is inherently collective, in that it does not exist for a single cavity in the hard-core limit Drummond and Walls (1980). We note that collective bistability exists in a variety of other driven and dissipative systems Bonfacio and Lugiato (1978); Bowden and Sung (1979); Carmichael et al. (1986); Sarkar and Satchell (1987); Nagorny et al. (2003); Elsässer et al. (2004); Armen and Mabuchi (2006); Lee et al. (2012); Marcuzzi et al. (2014), and was recently observed in an gaseous ensemble of Rydberg atoms Carr et al. (2013). Gases of Rydberg atoms are also predicted to exhibit AF order Lee et al. (2011); Hu et al. (2013); Hoening et al. (2014), though unlike the model we consider here, their interactions (due to the Rydberg blockade) are effectively antiferromagnetic in nature. Other works studying the hard-core DDBH model also predict AF order, though they consider variants of the model that include spatially varying drive fields Hartmann (2010), two-cavity pumping Joshi et al. (2013), and cross-Kerr terms Jin et al. (2013, 2014). Our system exhibits AF order in the absence of these features, despite the ferromagnetic nature of the couplings.

In the region enclosed by the dashed line in Fig. 1, the AF states acquire a limit cycle (LC) character Nicolis and Prigogine (1977); Lee et al. (2011); Qian et al. (2012), exhibiting periodic oscillatory dynamics at long times, and thus breaking the continuous time-translational symmetry of the system. The inset in Fig. 1 shows an example limit cycle projected onto a Bloch sphere. Interestingly, the limit cycle exists both as a unique steady state, and as part of a bistable pair (green region). We note that similar limit cycles were recently predicted for a driven-dissipative XYZ spin- model Chan et al. (2015).

To determine the validity of the two-sublattice ansatz, we perform a linear stability analysis on the spatially uniform steady states Le Boité et al. (2013); Joshi et al. (2013); Le Boité et al. (2014) (see Sec. A.2 in the appendix). In the light blue and light green regions in Fig. 1, the steady state is dynamically unstable to the formation of incommensurate () spin density waves, which cannot be captured within the two-sublattice ansatz. To understand the effects of this instability, we solve the inhomogeneous MF equations for a 1D chain, with randomly seeded initial states. In the light blue region, the instability results in the disappearance of the steady state. In the light green region, we find steady states that exhibit AF order over finite size domains, which are interrupted by domains of the state; we refer to this as frustrated AF order (f-AF). While the AF phase is the only true periodically ordered steady state in this region, domains remain stable if they are sufficiently small, and do not sample unstable wave vectors that often exist over a very narrow range of . The limit cycles remain mostly AF ordered, with some low amplitude, small- features, and the state of the pair ceases to exist in the light green LC region.

## Iii Beyond mean-field

It is important to understand how the rich physics predicted by the Gutzwiller MF theory survives in the presence of quantum and classical fluctuations, which exist in the true steady state of the DDBH model, and play particularly important roles in low dimensions. Toward this end, we employ a quantum trajectories algorithm to study finite 1D systems Dalibard et al. (1992); Dum et al. (1992); Plenio and Knight (1998); Daley (2014) (see Sec. B in the appendix). This method provides exact results for physical observables in the steady state under the ensemble averaging of trajectories.

We present results for a system with cavities and periodic boundary conditions in Fig. 2. In panel (a), we show the nearest-neighbor part of the connected correlation function , where , which is independent of . We choose to study correlations because the -components of the spins exhibit the strongest AF order in the MF results, and can be measured via correlated homodyne detection in an experimental setup. The blue region shows where is positive, and the red region shows where is negative, corresponding to AF nearest-neighbor correlations. The inset shows for and , represented by green diamonds (red squares). While both correlation functions exhibit nearest-neighbor AF order, the correlations for have incommensurate spin density wave character while the correlations for have a true AF character, switching from negative to positive values with each successive cavity. The incommensurate spin density wave correlations are present where the MF linear stability analysis predicts an incommensurate spin density wave instability of the steady state, corresponding to the f-AF region in Fig. 1.

We performed finite size scaling calculations for these and a number of other parameters, and find that accurately captures AF order of the DDBH model in the thermodynamic limit for a large region of the phase diagram; this is due to the short-range, exponentially decaying nature of . The black lines overlaid in Fig. 2 show the MF phase diagram boundaries. Interestingly, there is reasonably good agreement between the region of short-range AF correlations in the 1D quantum system and the MF results.

In Fig. 2(b), we show the normalized fluctuations in the total photon number , given by where and . Interestingly becomes anomalously large in the darker blue region, which has significant overlap with the MF bistability, indicating that photon number fluctuations become strongly correlated when the MF theory predicts collective bistability. The origin of the enhancement is revealed upon inspection of the trajectories themselves. In Fig. 3(a), we show the photon number as a function of time for cavities with and , in the region of enhanced (shown by the white circle in Fig. 2(b)). Interestingly, the trajectory exhibits collective switching between macroscopically distinguishable states, which resemble the MF steady states for these parameters. We plot the total photon number obtained via MF calculations as a function of for in panel (b) of Fig. 3; the values at are indicated by the grey markers. We show for the quantum trajectory as a function of time in panel (c). Here, fluctuates about two mean values for extended periods of time, which are interrupted by switching events that drive the system from one MF-like state to the other.

## Iv Collective bistability & switching

A recent work Mendoza-Arenas et al. (2015) showed that collective bistability in the driven-dissipative XY spin- model vanishes as the Gutzwiller approximation is systematically improved. Our results demonstrate that while is indeed unique Spohn (1977); Schirmer and Wang (2010), it retains clear signatures of the MF bistability, namely that the system dynamically switches between two macroscopically distinguishable configurations that resemble the MF steady states and ; in fact, this is a generic feature of bistable systems Armen and Mabuchi (2006); Kerckhoff et al. (2011); Lee et al. (2012); Hu et al. (2013); Dombi et al. (2013). The approach to in the bistable regime is then characterized by switching between these two states, and the rate of convergence is directly related to the average time spent in and between switching events; we refer to these times as and , respectively. In the master equation formalism, the asymptotic rate of convergence is set by the Liouvillian gap , where is the eigenvalue of with the smallest magnitude non-zero real part Verstraete et al. (2009); Cai and Barthel (2013); this suggests that and are intimately related Kinsler and Drummond (1991). The presence of true bistability in the MF theory reflects the fact that the Liouvillian gap vanishes at this level of approximation. The collective switching in the exact dynamics, and the concomitant opening of the Liouvillian gap, can be understood to result from both quantum and dissipation-induced (classical) fluctuations that are not included in the MF theory.

It is natural to expect that MF behavior can be recovered in the quantum system as its coordination number is increased. We explore this possibility by taking the spin couplings (or photon hopping) to be infinite-range. This corresponds to modifying the hopping term in Eq. (1) to , where . While this limit may seem unnatural for arrays of photonic cavities, it could be achieved using an external mirror, and it is in fact quite natural for other open quantum systems, such as ensembles of Rydberg atoms Lee et al. (2011); Hoening et al. (2014); Weimer (2015) or trapped ions Porras and Cirac (2004); Kim et al. (2010). We calculate the Liouvillian gap exactly by taking advantage of an efficient parameterization of the accessible space of density matrices (see Sarkar and Satchell (1987); Xu et al. (2013) and Sec. C in the appendix); we show the Liouvillian gap for as a function of and in Fig. 4(a). There is striking agreement between the exact quantum results and the MF results; where MF theory predicts collective bistability, the Liouvillian gap decreases to  Sarkar and Satchell (1987).

We plot for (12) cavities and as a function of in Fig. 4(b), shown by the light (dark) blue solid lines. Savage and Carmichael proposed a two-state toy model to describe a bistable system with a small Liouvillian gap Savage and Carmichael (1988), which has a gap . We extract values for and heuristically by measuring the time spent in the dark (1) and bright (2) states of the quantum trajectory simulations, and plot in Fig. 4(b), shown by the squares (circles) for (). Already at , there is excellent quantitative agreement between the exact Liouvillian gap and the results of the simple two state model as extracted from the quantum trajectories with infinite-range couplings. This provides a clear connection between the MF and quantum solutions in the bistable regime. We note that Kinsler and Drummond performed a similar analysis for the single-mode quantum parametric oscillator, and also found good quantitative agreement for large photon numbers Kinsler and Drummond (1991).

## V Discussion

The thermodynamic limit of the 1D DDBH model with nearest-neighbor interactions is challenging to study numerically, but we expect its dynamics to exhibit collective switching over finite-size domains. This behavior is reminiscent of equilibrium systems that, while exhibiting a first-order phase transition in higher dimensions, fail to do so in 1D. Whether or not mean-field bistability is associated with a true first-order phase transition in higher spatial dimensions is an interesting question that warrants further study. Finally, we note that we have studied the soft-core DDBH model (with finite ), and identified features that are directly analogous to those discussed here. A comprehensive study of the soft-core DDBH model is the subject of future work.

###### Acknowledgements.
We thank Mohammad Maghrebi, Sarang Gopalakrishnan, and Jeremy Young for insightful discussions. RW, KM, AG, and MH thank the KITP for hospitality. We acknowledge partial support from the NSF, ONR, ARO, ARL, AFOSR, NSF PIF, NSF PFC at the JQI, and the Sloan Foundation. This research was supported in part by the NSF under Grant No. NSF PHY11-25915.

## Appendix A Equations of motion in the U→∞ limit

In this appendix, we provide technical details to support the theory and numerics in the main text. Though many of our results are valid more generally for soft-core bosons, here we specialize to the limit of hard-core bosons, valid when , where is the local interaction strength in Eq. (2) of the main text. Hard-core bosons can be conveniently mapped onto spins via an inverse Holstein-Primakoff transformation,

 ^a†i→^σ+i  and  ^ai→^σ−i, (2)

where and , is the identity operator, and are Pauli matrices that act on cavity . Following this transformation, the master equation becomes

 Missing or unrecognized delimiter for \big (3)

where is the system Hamiltonian,

 ^H=−J2z∑⟨i,j⟩(^σxi^σxj+^σyi^σyj)+Ω∑i^σxi−μ2∑i^σzi. (4)

Here the summations run over , and the notation indicates an additional sum over all cavities that are coupled to cavity ; the number of cavities are quantified by the coordination number . Equations (3) and (4) describe a driven-dissipative spin- XY model, with isotropic interactions and in the presence of a homogeneous applied field.

To study the steady states of Eq. (3), it is convenient to use equations of motion for the spin components ; these equations are readily derived by taking . Using the Pauli matrix commutation relations where is the Levi-Civita symbol, we find

 ∂tσxi =−2Jz∑⟨i,j⟩⟨^σzi^σyj⟩+μσyi−γ2σxi ∂tσyi =2Jz∑⟨i,j⟩⟨^σzi^σxj⟩−μσxi−2Ωσzi−γ2σyi ∂tσzi =2Jz∑⟨i,j⟩⟨^σxi^σyj−^σyi^σxj⟩+2Ωσyi−γ(σzi+1). (5)

### a.1 Gutzwiller mean-field approximation

In the Gutzwiller mean-field (MF) approximation, the density matrix is assumed to factorize over all cavities, . In the spin formalism, this corresponds to the factorization of all non-local two-spin expectation values. Thus, in the MF approximation, Eqs. (A) become

 ∂tσxi =−2Jzσzi∑⟨i,j⟩σyj+μσyi−γ2σxi ∂tσyi =2Jzσzi∑⟨i,j⟩σxj−μσxi−2Ωσzi−γ2σyi ∂tσzi =2Jzσxi∑⟨i,j⟩σyj−2Jzσyi∑⟨i,j⟩σxj+2Ωσyi−γ(σzi+1). (6)

In the majority of the phase diagram in Fig. 1 of the main text, the steady states are obtained by evolving these equations numerically using a order Runge-Kutta algorithm.

### a.2 Linear stability analysis

The phase diagram in Fig. 1 of the main text also shows results from a linear stability analysis on spatially uniform steady states within the MF approximation, which is carried out as follows. Spatially uniform states have the property for all . Using this as an ansatz for Eqs. (A.1), we find the following equations for the uniform steady states,

 0 =−2Jσzσy+μσy−γ2σx 0 =2Jσzσx−μσx−2Ωσz−γ2σy 0 =2Ωσy−γ(σz+1), (7)

which can be solved analytically to find the spin configurations of the uniform steady states. For example, we plot solutions for and as a function of in Fig. 5(a); these solutions exhibit collective bistability, where the dynamically stable steady state solutions are shown by the solid black lines. The dashed black line in this figure indicates a dynamically unstable solution, which we discuss in more detail below.

The ansatz of spatially uniform steady states is not appropriate under all circumstances, for example, if the steady state develops a spatial order that breaks the discrete translational symmetry of the DDBH model. The presence of such spatially ordered states can sometimes be captured by an instability of the uniform steady state. To explore this possibility, we perform a linear stability analysis on the uniform steady state solutions of Eqs. (A.2), specializing to one-dimensional systems with . In a system with infinite spatial extent, this is accomplished by adding a small plane wave perturbation to the uniform steady state of the form , where and are the solutions of Eqs. (A.2), are cavity indices, and is the wave number of the perturbation. Linearizing in , we find the equation

 ∂tδ=\mathdsMδ, (8)

where

 Unknown environment '% (9)

The matrix has eigenvalues that depend on the wave number . When the real part of is negative for all , the uniform steady state is dynamically stable. When some acquires a positive real part, this signifies an instability of the uniform steady state, and the wave number at which is maximum corresponds to the mode that is maximally unstable, and dominates the dynamics of the instability. For example, in Fig. 5(b) we plot as a function of for three sets of parameters, all with . The blue lines correspond to one of three solutions with , shown by the blue circle in panel (b) of this figure. These modes are unstable at , indicating a global instability of the uniform steady state solution. The green and red lines correspond to parameters exhibiting spin density wave instabilities. The green line exhibits an incommensurate spin density wave instability, as the modes are stable at but unstable in a small region of . The red line exhibits an antiferromagnetic (AF) spin density instability, as the most unstable eigenvalue occurs for .

and antiferromagnetic (AF) instabilities, respectively.

## Appendix B Exact numerical solution of Eq. (3)

We employ a quantum trajectories algorithm to solve the master equation (3), which is a powerful, exact method for studying open quantum systems that relies on treating the system-environment coupling (the dissipator in the master equation) stochastically. Here, we present the algorithm for generating a quantum trajectory using quantum jumps. We begin by defining the effective non-Hermetian Hamiltonian,

 ^Heff=^H−iγ2∑i^ni, (10)

so Eq. (3) can be written as

 ∂t^ρ=−i(^Heff^ρ−^ρ^H†eff)+γ∑i^ai^ρ^a†i. (11)

The quantum trajectories method amounts to evolving a wave function, as opposed to a density matrix, under , while treating the right-most “recycling” term in Eq. (11) stochastically.

Consider a wave function at time that evolves under . For a small time step , the wave function at can be written as (using Euler integration)

 |ψ′(t+dt)⟩=(1−i^Heffdt)|ψ(t)⟩. (12)

The norm of is not conserved in real-time evolution under due to its being non-Hermetian, so in general we have

 ⟨ψ′(t+dt)|ψ′(t+dt)⟩⟨ψ(t)|ψ(t)⟩=1−δp. (13)

In the quantum trajectories formalism, we choose the wave function at stochastically. With probability , we choose

 |ψ(t+δt)⟩=|ψ′(t+δt)⟩√1−δp, (14)

and with probability we take our next state to be one that emitted a photon from cavity ,

 |ψ(t+δt)⟩=^ai|ψ(t)⟩√δpi/δt, (15)

where the jump-site is chosen with probability

 Πi=δpiδp=δt⟨ψ(t)|^ni|ψ(t)⟩δp, (16)

and .

## Appendix C Infinite range interactions

In general, long-range interactions introduce considerable difficulties when numerically studying the dynamics of many-body systems. However, when all pairs of spins interact with the same strength, the model becomes symmetric under the exchange of any two spins, which allows for an extremely efficient parametrization of the accessible Hilbert space. In the absence of spontaneous emission, permutation symmetry of the Hamiltonian restricts the dynamics of initially permutation symmetric pure states to the subspace of collective spin-states, often referred to as Dicke states, and the dimension of the accessible Hilbert space is therefore . For an open quantum system governed by a permutation symmetric Liouvillian, dynamics is restricted to the space of permutation symmetric density matrices, which can be parametrized in terms of symmetrized direct products of Pauli matrices. The space of permutation symmetric, Hermitian matrices over the product Liouville space of spins is spanned by basis states of the following form,

 ^M(n) =12N∑χ(^σxχ(1)⊗⋯⊗^σxχ(nx))⊗(^σyχ(nx+1)⊗⋯⊗^σyχ(nx+ny)) ⊗(^σzχ(nx+ny+1)⊗⋯⊗^σzχ(nx+ny+nz))⊗(^σ0χ(nx+ny+nz+1)⊗⋯⊗^σ0χ(N)). (17)

Here the parameters are positive integers that are arbitrary up to the constraint , is the identity matrix on the Hilbert space of spin , and is a permutation of integers . The most general permutation symmetric Hermitian matrix can be written

 ^ρ=∑nc(n)^M(n). (18)

Note that has unity trace, while all the other basis states are traceless; hence can only be a valid density matrix if . The number of unconstrained coefficients, and hence the dimensionality of the space that must be considered in the dynamics, is the number of ways of choosing such that . It is straightforward to check that there are such , where

 D=(N+3)(N+2)(N+1)3!−1, (19)

hence the dimension of the accessible space of density matrices is . Even after the restriction , a matrix of the form in Eq. (18) is not guaranteed to be a valid density matrix, as it may have negative eigenvalues (equivalently, it need not satisfy ). However, the master equation provides a positive map, and preserves the positivity of the density matrix eigenvalues dynamically. Therefore if we start with an initial state that is a valid density matrix, it will remain restricted to the space of valid density matrices during the time evolution.

The permutation symmetry of the Liouvillian endows it with a block diagonal structure, with one of the blocks spanned by the states in Eq. (18). Dynamics of an initial density matrix within this block can therefore be calculated by determining the action on all states :

 L(^M(n))=∑mLm,n^M(m). (20)

The coefficients are straightforward to compute, and using this Eq. (20) the master equation can be written

 ∂tc(n)=Ln,mc(m). (21)

This set of first-order, linear, ordinary differential equations is straightforward to integrate for .

### References

1. M. E. Fisher, Rev. Mod. Phys. 46, 597 (1974).
2. M. E. Fisher, Rev. Mod. Phys. 70, 653 (1998).
3. S. Sachdev, Quantum phase transitions (Wiley Online Library, 2007).
4. J. Kasprzak, M. Richard, S. Kundermann, A. Baas, P. Jeambrun, J. M. J. Keeling, F. M. Marchetti, M. H. Szymańska, R. André, J. L. Staehli, V. Savona, P. B. Littlewood, B. Deveaud,  and L. S. Dang, Nature 443, 409 (2006).
5. S. Hofferberth, I. Lesanovsky, B. Fischer, T. Schumm,  and J. Schmiedmayer, Nature 449, 324 (2007).
6. A. Amo, D. Sanvitto, F. P. Laussy, D. Ballarini, E. del Valle, M. D. Martin, A. Lemaître, J. Bloch, D. N. Krizhanovskii, M. S. Skolnick, C. Tejedor,  and L. Vina, Nature 457, 291 (2008).
7. A. Amo, J. Lefrére, S. Pigeon, C. Adrados, C. Ciuti, I. Carusotto, R. Houdré, E. Giacobino,  and A. Bramati, Nat. Phys. 5, 805 (2009).
8. H. Deng, H. Haug,  and Y. Yamamoto, Rev. Mod. Phys. 82, 1489 (2010).
9. A. Polkovnikov, K. Sengupta, A. Silva,  and M. Vengalattore, Rev. Mod. Phys. 83, 863 (2011).
10. I. Carusotto and C. Ciuti, Rev. Mod. Phys. 85, 299 (2013).
11. M. O. Scully and W. E. Lamb, Phys. Rev. 159, 208 (1967).
12. V. DeGiorgio and M. O. Scully, Phys. Rev. A 2, 1170 (1970).
13. E. G. D. Torre, S. Diehl, M. D. Lukin, S. Sachdev,  and P. Strack, Phys. Rev. A 87, 023831 (2013).
14. P. Kirton and J. Keeling, Phys. Rev. Lett. 111, 100404 (2013).
15. M. Hafezi, P. Adhikari,  and J. M. Taylor, Phys. Rev. B 92, 174305 (2015).
16. M. F. Maghrebi and A. V. Gorshkov, Phys. Rev. B 93, 014307 (2016).
17. S. M. Spillane, G. S. Pati, K. Salit, M. Hall, P. Kumar, R. G. Beausoleil,  and M. S. Shahriar, Phys. Rev. Lett. 100, 233602 (2008).
18. E. Vetsch, D. Reitz, G. Sagué, R. Schmidt, S. T. Dawkins,  and A. Rauschenbeutel, Phys. Rev. Lett. 104, 203603 (2010).
19. K. Hennessy, A. Badolato, M. Winger, D. Gerace, M. Atatüre, S. Gulde, S. Fält, E. L. Hu,  and A. Imamoğlu, Nature 445, 896 (2007).
20. T. G. Tiecke, J. D. Thompson, N. P. d. Leon, L. R. Liu, V. Vuletic,  and M. D. Lukin, Nature 508, 241 (2014).
21. M. Hafezi, D. E. Chang, V. Gritsev, E. Demler,  and M. D. Lukin, Phys. Rev. A 85, 013822 (2012).
22. J. D. Thompson, T. G. Tiecke, N. P. de Leon, J. Feist, A. V. Akimov, M. Gullans, A. S. Zibrov, V. Vuletic,  and M. D. Lukin, Science 340, 1202 (2013).
23. A. V. Gorshkov, J. Otterbach, M. Fleischhauer, T. Pohl,  and M. D. Lukin, Phys. Rev. Lett. 107, 133602 (2011).
24. T. Peyronel, O. Firstenberg, Q.-Y. Liang, S. Hofferberth, A. V. Gorshkov, T. Pohl, M. D. Lukin,  and V. Vuletic, Nature (London) 488, 57 (2012).
25. A. A. Houck, H. E. Türeci,  and J. Koch, Nat. Phys. 8, 292 (2012).
26. F. Nissen, S. Schmidt, M. Biondi, G. Blatter, H. E. Türeci,  and J. Keeling, Phys. Rev. Lett. 108, 233603 (2012).
27. J. Reftery, D. Sadri, S. Schmidt, H. E. Türeci,  and A. A. Houck, Phys. Rev. X 4, 031043 (2014).
28. R. Barends, L. Lamata, J. Kelly, L. Garcia-Alvarez, A. G. Fowler, A. Megrant, E. Jeffrey, T. C. White, D. Sank, J. Y. Mutus, B. Campbell, Y. Chen, Z. Chen, B. Chiaro, A. Dunsworth, I. C. Hoi, C. Neill, P. J. J. O’Malley, C. Quintana, P. Roushan, A. Vainsencher, J. Wenner, E. Solano,  and J. M. Martinis, Nat. Comm. 6, 1 (2015).
29. D. E. Chang, V. Gritsev, G. Morigi, V. Vuletić, M. D. Lukin,  and E. A. Demler, Nat. Phys. 4, 884 (2008).
30. D. E. Chang, V. Vuletić,  and M. D. Lukin, Nat. Phot. 8, 685 (2014).
31. M. J. Hartmann, F. G. S. L. Brandao,  and M. B. Plenio, Nat. Phys. 2, 849 (2006).
32. M. J. Hartmann, F. G. S. L. Brandao,  and M. B. Plenio, Las. & Photon. Rev. 2, 527 (2008).
33. I. Carusotto, D. Gerace, H. E. Tureci, S. D. Liberato, C. Ciuti,  and A. Imamoğlu, Phys. Rev. Lett. 103, 033601 (2009).
34. A. Le Boité, G. Orso,  and C. Ciuti, Phys. Rev. Lett. 110, 233601 (2013).
35. A. Le Boité, G. Orso,  and C. Ciuti, Phys. Rev. A 90, 063821 (2014).
36. H.-P. Breuer and F. Petruccione, The theory of open quantum systems (Oxford University Press, 2002).
37. H. J. Carmichael, An Open Systems Approach to Quantum Optics (Springer, Berlin, 1993).
38. M. P. A. Fisher, P. B. Weichman, G. Grinstein,  and D. S. Fisher, Phys. Rev. B 40, 546 (1989).
39. A. D. Greentree, C. Tahan, J. H. Cole,  and L. C. L. Hollenberg, Nat. Phys. 2, 856 (2006).
40. D. G. Angelakis, M. F. Santos,  and S. Bose, Phys. Rev. A 76, 031805(R) (2007).
41. M. J. Hartmann and M. B. Plenio, Phys. Rev. Lett. 99, 103601 (2007).
42. D. Rossini and R. Fazio, Phys. Rev. Lett. 99, 186401 (2007).
43. J. Koch and K. L. Hur, Phys. Rev. A 80, 023811 (2009).
44. S. Diehl, A. Tomadin, A. Micheli, R. Fazio,  and P. Zoller, Phys. Rev. Lett. 105, 015702 (2010).
45. A.-W. de Leeuw, O. Onishchenko, R. A. Duine, ,  and H. T. C. Stoof, Phys. Rev. A 91, 033609 (2015).
46. D. S. Rokhsar and B. G. Kotliar, Phys. Rev. B 44, 10328 (1991).
47. J. Dalibard, Y. Castin,  and K. Mølmer, Phys. Rev. Lett. 68, 580 (1992).
48. R. Dum, P. Zoller,  and H. Ritsch, Phys. Rev. A 45, 4879 (1992).
49. M. B. Plenio and P. L. Knight, Rev. Mod. Phys. 70, 101 (1998).
50. A. J. Daley, Adv. Phys. 63, 77 (2014).
51. T. Holstein and H. Primakoff, Phys. Rev. 58, 1098 (1940).
52. C. Joshi, F. Nissen,  and J. Keeling, Phys. Rev. A 88, 063835 (2013).
53. P. D. Drummond and D. F. Walls, J. of Phys. A 13, 725 (1980).
54. R. Bonfacio and L. A. Lugiato, Phys. Rev. A 18, 1129 (1978).
55. C. M. Bowden and C. C. Sung, Phys. Rev. A 19, 2392 (1979).
56. H. J. Carmichael, J. S. Satchell,  and S. Sarkar, Phys. Rev. A 34, 3166 (1986).
57. S. Sarkar and J. S. Satchell, Europhys. Lett. 3, 797 (1987).
58. B. Nagorny, T. Elsässer,  and A. Hemmerich, Phys. Rev. Lett. 91, 153003 (2003).
59. T. Elsässer, B. Nagorny,  and A. Hemmerich, Phys. Rev. A 69, 033403 (2004).
60. M. A. Armen and H. Mabuchi, Phys. Rev. A 73, 063801 (2006).
61. T. E. Lee, H. Häffner,  and M. C. Cross, Phys. Rev. Lett. 108, 023602 (2012).
62. M. Marcuzzi, E. Levi, S. Diehl, J. P. Garrahan,  and I. Lesanovsky, Phys. Rev. Lett. 113, 210401 (2014).
63. C. Carr, R. Ritter, C. G. Wade, C. S. Adams,  and K. J. Weatherill, Phys. Rev. Lett. 111, 113901 (2013).
64. T. E. Lee, H. Häffner,  and M. C. Cross, Phys. Rev. A 84, 031402(R) (2011).
65. A. Hu, T. E. Lee,  and C. W. Clark, Phys. Rev. A 88, 053627 (2013).
66. M. Hoening, W. Abdussalam, M. Fleischhauer,  and T. Pohl, Phys. Rev. A 90, 021603(R) (2014).
67. M. J. Hartmann, Phys. Rev. Lett. 104, 113601 (2010).
68. J. Jin, D. Rossini, R. Fazio, M. Leib,  and M. J. Hartmann, Phys. Rev. Lett. 110, 163605 (2013).
69. J. Jin, D. Rossini, M. Leib,  and M. J. H. and. R. Fazio, Phys. Rev. A 90, 023827 (2014).
70. G. Nicolis and I. Prigogine, Self-Organization in Nonequilibrium Systems (Wiley, 1977).
71. J. Qian, G. Dong, L. Zhou,  and W. Zhang, Phys. Rev. A 85, 065401 (2012).
72. C.-K. Chan, T. E. Lee,  and S. Gopalakrishnan, Phys. Rev. A 91, 051601(R) (2015).
73. J.-J. Mendoza-Arenas, S. R. Clark, S. Felicetti, G. Romero, E. Solano, D. G. Angelakis,  and D. Jaksch, Phys. Rev. A 93, 023821 (2015).
74. H. Spohn, Lett. Math. Phys. 2, 33 (1977).
75. S. G. Schirmer and X. Wang, Phys. Rev. A 81, 062306 (2010).
76. J. Kerckhoff, M. A. Armen,  and H. Mabuchi, Optics Express 19, 24468 (2011).
77. A. Dombi, A. Vukics,  and P. Domokos, J. Phys. B 46, 224010 (2013).
78. F. Verstraete, M. M. Wolf,  and J. I. Cirac, Nat. Phys. 5, 633 (2009).
79. Z. Cai and T. Barthel, Phys. Rev. Lett. 111, 150403 (2013).
80. P. Kinsler and P. D. Drummond, Phys. Rev. A 43, 6194 (1991).
81. H. Weimer, Phys. Rev. Lett. 114, 040402 (2015).
82. D. Porras and J. I. Cirac, Phys. Rev. Lett. 92, 207901 (2004).
83. K. Kim, M.-S. Chang, S. Korenblit, R. Islam, E. E. Edwards, J. K. Freericks, G.-D. Lin, L.-M. Duan,  and C. Monroe, Nature 465, 590 (2010).
84. M. Xu, D. A. Tieri,  and M. J. Holland, Phys. Rev. A 87, 062101 (2013).
85. C. M. Savage and H. J. Carmichael, IEEE J. Quantum Electron. 24, 1495 (1988).
You are adding the first comment!
How to quickly get a good reply:
• Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
• Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
• Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
The feedback must be of minumum 40 characters