## Abstract

We investigate the mean-field dynamics of a system of interacting photons in an array of coupled cavities in the presence of dissipation and disorder. We follow the evolution of an initially prepared Fock state, and show how the interplay between dissipation and disorder affects the coherence properties of the cavity emission, and show that these properties can be used as signatures of the many-body phase of the whole array.

## 1. Introduction

The idea of understanding the behaviour of complex quantum many-body systems using experimentally controllable *quantum* *simulators* can be traced back to a pioneering keynote speech given by Richard Feynman [1]. After 30 years, quantum simulation is now a thriving field of research [2,3], driven by the increasing ability to design and fabricate controllable quantum systems, in contexts ranging from superconducting circuits [4] to ultracold atoms [5] or trapped ions [6]. These systems allow the realization of archetypal models and the exploration of new physical regimes. An area of recently developing interest has been coupled cavity arrays (CCAs), lattices of coupled matter-light systems, providing highly tunable but dissipative quantum systems [7–11].

Physical realizations of CCAs have been proposed in a variety of systems, such as photonic crystal nanocavities [12], coupled optical waveguides [13] or lattices of superconducting resonators operating in the microwave regime [4,14,15]. While experiments have not yet probed the collective behaviour predicted in large-scale arrays, progress towards such realizations is very encouraging. For these different systems, the radiation modes involved range from microwave to optical frequency; we will nonetheless refer to these as coupled ‘matter-light’ systems in the following, taking ‘light’ to refer to the radiation modes.

A generic CCA consists of a lattice of cavities, each supporting a confined photon quasi-mode. We refer to quasi-modes [16], because the finite quality of the cavities implies there will be coupling to the outside world. As shown in figure 1, the cavities are coupled through photon hopping. A purely optical system would be entirely linear, and thus unable to simulate interacting many-body quantum systems. To introduce nonlinearity requires coupling to matter (e.g. suitable optical emitters such as semiconductor quantum dots or superconducting qubits, indicated by blue circles in figure 1). This leads to a system of photons hopping on a lattice, with an on-site nonlinearity and on-site losses.

A wide variety of different microscopic models can be realized depending on the precise design of the array [10,11,15]. We consider here the archetypal case where the CCA can be mapped [10] onto the Bose–Hubbard model [17]. This corresponds to the regime of weak matter-light coupling, i.e. strong detuning between cavity mode and matter degrees of freedom [18].

In the absence of dissipation, the ground-state phase diagram of the Bose–Hubbard model has been extensively studied [5]: at zero temperature, a quantum phase transition results from the competition between on-site nonlinearities and inter-site photon hopping. When the nonlinearities are strong and prevail over the hopping, the photons are localized by interactions, leading to an insulating Mott phase; in the opposite regime, when photon hopping dominates, a superfluid phase characterized by long-range coherence occurs. Coupled matter-light systems however are naturally studied under non-equilibrium conditions, as there are invariably photon losses. As such, a steady state in a CCA requires external pumping. The steady state of cavity arrays, resulting from the competition of external driving and dissipation, has been the subject of much recent theoretical interest [18–23].

Besides the steady state, the transient dynamics of cavity arrays may bring additional interesting information. By engineering suitable laser pulse sequences [24,25], it is possible to prepare a specific initial state and follow its subsequent evolution by analysing the properties of the light escaping from the system. This is the situation considered in [26], considering the evolution following preparation of a Fock state. Such a protocol is equivalent to studying a quantum quench in an open system. Quantum quenches in closed systems are an intensively studied paradigm of non-equilibrium dynamics of many-body systems [27], and in the closed Bose–Hubbard model have been already addressed both theoretically and experimentally [27–30]. Quantum quenches in open spin chains have also revealed the possibility of non-trivial decay of coherence [31]. Cavity arrays appear to be ideally suited to explore the transient dynamics following a quench.

Remarkably, even in the lossy system, the dynamics following such a quench clearly map out a superfluid–insulator transition [26]. In fact, in mean-field theory, these can be directly related to the equilibrium phase boundary [26,32]. In particular, by rescaling correlation functions by the decaying density one finds that the behaviour at long times is distinct for values of hopping in the superfluid and insulating phases—the rescaled coherence vanishes in the insulating phase, but attains a non-zero asymptote in the superfluid phase. As discussed further below, this can be traced back to the way in which for the closed system quench, the linear stability of an initial Mott state reproduces the equilibrium phase diagram [33]. Given the notable ability of the rescaled correlations of the open system to reproduce the equilibrium phase boundary, and the significant interest in the disordered Bose–Hubbard model [17,34–36], a natural question to ask is how the open system quench dynamics are affected by disorder. This is the question we begin to address in this paper. Building on the results of Tomadin *et al.* [26], we analyse the non-equilibrium mean-field dynamics of an array of nonlinear coupled cavities in the presence of photon leakage and disorder in the on-site cavity energies.

The paper is organized as follows: in §2, we introduce the model used to describe the dissipative and disordered cavity array; in §3, we analyse the dynamics of the correlation functions of the cavity array, first summarizing the results for the ideal *clean case* and then assuming a Gaussian distribution for the on-site cavity energies, and compare the results; in §4, we briefly discuss an anomalous behaviour of the second-order correlation function in the presence of disorder; in §5, we summarize our conclusions.

## 2. Model system and initial state

As discussed above, we consider the Bose–Hubbard model [17] described by the Hamiltonian
*i*th site, and the corresponding number operator is *i*th cavity is *ε*_{i}, *J* denotes the amplitude for inter-site photon hopping, whereas *U* represents the on-site nonlinearity resulting from the coupling to matter. Including also the presence of Markovian photon loss leads to the master equation
*κ*)^{−1}. We consider the case where *U*,*κ*,*J* are independent of site, but there may be disorder in the energies *ε*_{i}. We discuss below the ‘clean’ case in which all *ε*_{i} are the same, and the disordered case where we choose cavity energies *ε*_{i} to be drawn from Gaussian distribution having mean value *σ*^{2}, with *σ* representing the strength of the disorder. The mean value

We explore the quench dynamics of this model in the mean-field approximation *z*, where *z* is the number of nearest neighbours of each site. The dynamics of each cavity is governed by the master equation
*z* nearest neighbours of site *i*. In the clean case, all sites are equivalent and so *z*. In the disordered case, each site evolves separately, and the connectivity of the lattice does affect the dynamics.

Our goal in the following discussion is to study the non-equilibrium dynamics following preparation of the cavity array in a product of Fock states, i.e. *ρ*_{i}(0)=|*Ψ*_{0}〉〈*Ψ*_{0}|, where *η*≪1 will grow, and drive the array to a different asymptotic state.

## 3. Dynamics of correlation functions

### (a) Clean case

Before exploring the role played by disorder, we first summarize the results of the clean case [26] and present also a discussion of the initial instability, first considered in [32] and here thoroughly discussed. In the absence of dissipation, it was shown [33] that when *zJ*/*U*>*zJ*/*U*|_{cr}, where *zJ*/*U*|_{cr} is the critical value corresponding to the equilibrium superfluid–insulator phase transition, an initial Fock state is linearly unstable. As such, the existence of this instability can be used to trace the equilibrium transition between the Mott and the superfluid phase [32,33] and, within our mean-field analysis, a transition survives also in the presence of dissipation.

To see this instability for the clean case, one may write coupled equations for *ρ*_{n0−1,n0},*ρ*_{n0,n0+1} (where *n*_{0} is the occupation of the initial state, under the approximation *ρ*_{n0,n0}≃1, and that all other elements of *ρ* are negligible). Denoting *X*=(*ρ*_{n0−1,n0},*ρ*_{n0,n0+1})^{T}, one finds that *X* obeys the equation ∂_{t}*X*=*MX* with
*M* become positive, as this corresponds to exponential growth of fluctuations. In the absence of dissipation (*κ*=0), the eigenvalues *ξ* are given by
*ξ*=*iμ*, one may show that Det(*iμ*1−*M*)=0 is equivalent to
*n*_{0}th Mott lobe. As *zJ* increases, the critical values of *μ* in the equilibrium phase boundary approach each other, and at large enough *zJ*, there is no longer any real value *μ* that can satisfy equation (3.3). As such, the instability of the Fock state corresponds to the locations of the tips of the equilibrium Mott lobes, *κ* can be found in [32]. Further, this analysis allows one to study how *zJ*/*U*|_{cr} changes when the initial state is a statistical admixture, e.g. when *ρ*(0)=*α*^{2}|*n*_{0}〉〈*n*_{0}|+*β*^{2}|*n*_{0}+1〉〈*n*_{0}+1| (*α*^{2}+*β*^{2}=1). In this case, for *κ*=0, the eigenvalues *ξ* are given by the solutions of
*linear* instability of the normal state evolves continuously. This is answered by figure 2, showing that the critical hopping *zJ*/*U*|_{cr} for this instability (obtained solving equation (3.4) for different *n*_{0}) does evolve continuously as a function of *β*^{2}. In the remainder of this manuscript, we restrict to the case *β*=0.

The consequence of this instability can be seen in the time evolution of the coherence *ψ*_{i}(*t*) for a clean system with *zJ*/*U*>*zJ*/*U*|_{cr}. Figure 3 is plotted for parameters *κ*=10^{−2}*U* and *zJ*=3*U*, with an initial state with *n*_{0}=1, *η*=10^{−5}—we consider this same initial state throughout the remainder of the manuscript. In the initial time range *ψ*(*t*) grows exponentially even though the photon population *n*(*t*)=*n*(0) *e*^{−2κt}, see the bottom inset in figure 3*a*. This exponential growth leads to a regime beyond the validity of linearization, which features underdamped relaxation oscillations of *ψ*(*t*). Note that for the conservative case *κ*=0 this is replaced by undamped periodic oscillations [33]. At longer times, *tκ*>1, the oscillations are damped out, but the amplitude *ψ*(*t*) also decays to zero as *a*. In contrast, if one considered instead a case where *J*/*U*<*J*/*U*|_{cr}, then there would be no instability. One should however note that while the distinction between initial stability and instability depends only on the parameters of the model, the asymptotic value of *η* chosen—for smaller *η*, the instability takes longer to reach the nonlinear regime, and so the average photon number at this point is smaller, affecting the final state reached. Further information about the state can also be found from the correlation function *n*_{0}=1). One may in fact show that if *J*=0, *g*_{2}(*t*)=*g*_{2}(0)=1−1/*n*_{0} remains fixed at the value determined by the initial Fock state.

Because the rescaled field *g*_{2}(*t*) approach asymptotic values at late times, the behaviour for a given set of initial conditions can be characterized by these values. Formally, these are extracted by finding the time-averaged values 〈|*ψ*|〉_{t}, 〈*g*_{2}〉_{t}, averaging over a time window that neglects the initial transients as proposed in [26]. Such an approach is illustrated in figure 4, which shows the time-averaged 〈|*ψ*|〉_{t} and 〈*g*_{2}〉_{t} in the time interval 1<*tκ*<2, as a function of the hopping amplitude *J*/*U* and for three different values of the photon decay rate *κ*. One may clearly see that below a threshold value of *zJ*/*U* both 〈|*ψ*|〉_{t} and 〈*g*_{2}〉_{t} vanish. As *κ*→0, this threshold approaches the equilibrium superfluid–insulator transition which occurs at *zJ*/*U*|_{cr}∼0.17 for *n*_{0}=1. While these results occur in mean-field theory, similar results have been reported by Tomadin *et al.* [26] using a cluster mean-field approach.

### (b) Disordered case

We now explore the role played by the on-site cavity disorder *ε*_{i} in the non-equilibrium dynamics. As discussed above, we thus solve the Liouville problem equation (2.3) drawing the cavity energies *ε*_{i} from a Gaussian distribution with standard deviation *σ*. In order to characterize the properties of the ensemble, rather than those specific to a particular realization, we average the expectations |*ψ*| and *g*_{2} over different realizations of disorder, and additionally average over all sites within a given realization.

Figure 5 shows the time evolution of both the order parameter (*a*) and the second-order correlation function (*b*) for a one-dimensional array consisting of *N*=48 cavities, averaged over 10 realizations of the energies. Except for the distribution of site energies *ε*_{i}, all other parameters are as in figure 3. The site energies are drawn from Gaussian distributions with *σ*=0.1*U* (blue, dark grey line) and 0.3*U* (red, light grey line). At short times *tκ*<0.1, the dynamics is characterized by the same linear instability as is seen in the clean case (black line), and both *ψ*(*t*) and *g*_{2}(*t*) increase exponentially. At later times, however, disorder strongly modifies the behaviour. The oscillations seen previously in the clean case are washed out by averaging over different cavities, and so a quasi-steady-state value is reached earlier.

As in the clean case, the appearance of a plateau at late times suggests it is possible to characterize the evolution by its asymptotic value (see however the discussion in §4). Figure 6 shows the time-integrated |*ψ*| (*a*) and *g*_{2} (*b*) as the ratio *J*/*U* is varied at a fixed photon dissipation constant *κ*=10^{−2}*U* and for increasing values of disorder strength. As in the clean case, a threshold value of *zJ*/*U* is required before the instability occurs. This threshold value of *zJ*/*U* for the instability of the *ψ*=0 state appears to increase with increasing disorder. In equilibrium, there is a ‘Bose glass’ phase between the Mott insulator and the superfluid [17], where particles are no longer localized by interactions, but are instead localized by disorder. The critical hopping *zJ*/*U* for the equilibrium transition between the Bose glass and superfluid phase increases with hopping, and so our observation of increasing critical *zJ*/*U*. The increasing critical *zJ*/*U* we observe for the threshold in the open system is consistent with this.

## 4. Rare site statistics at late times

While in the clean case, the plateau reached by *g*_{2}(*t*) around *tκ*≃2 reflects the asymptotic time dependence, this turns out not to be the case for the disordered lattice. Despite the appearance of an apparent plateau seen in figure 5, this only indicates a temporary plateau. At later times, the values of 〈*g*_{2}(*t*)〉_{dis} starts to rise further as shown in figure 7. Intriguingly, this rise of 〈*g*_{2}(*t*)〉_{dis} at late times in fact reflects the existence of rare sites with large and exceptionally large values of *g*_{2}(*t*) that dominate the disorder average. This is illustrated in figure 8, which shows the probability density of *g*_{2} (*P*_{g2}) calculated for the same hopping and decay constants used in figure 7 using a sample of 2000 values of *g*_{2}, generated after the simulation of a one-dimensional array of 100 cavities for 20 different realizations of disorder with *σ*=0.3*U* (see the red, light grey curve in figure 7). The right-hand panel shows that, at late times (*tκ*=6), some sites exhibit large values (>5) of *g*_{2}. The presence of these large values of *g*_{2} are the reason that 〈*g*_{2}(*t*)〉_{dis} continues to evolve, and thus why figure 6*b* does not indicate a steady state. In the clean system, *g*_{2} has the same value for all sites, as the state is translationally invariant, thus the appearance of these anomalous values of *g*_{2} is a consequence of features of the disorder landscape as we discuss next. This behaviour can be better understood by looking at figure 9, which illustrates for a smaller hopping *zJ*=0.5*U* and *κ*=0.01*U*, how the disorder average 〈*g*_{2}(*t*)〉_{dis} rises when *σ*=0.3*U* (see the red, light grey curve) as a result of the occurrence of rare sites with large values of *g*_{2}, as shown by the probability density evaluated at *tκ*=1 and *tκ*=6, shown in figure 7. It is worth noting that all the anomalous sites (i.e. with *g*_{2}≫1) are local minima of the potential (energy landscape), but not all the local minima exhibit such anomalous values of *g*_{2}. Further, following the evolution to much later times becomes difficult, as the large ratio between different elements of the density matrix introduces numerical errors.

## 5. Conclusion

We have studied the time evolution of the disordered open Bose–Hubbard model following an initially prepared Mott state. As for the clean case, a dynamical transition does occur between small and large values of hopping, signalled by the asymptotic behaviour of the rescaled field *zJ*/*U*|_{cr}, i.e. does increase the hopping required for the superfluid instability to develop. However, even within mean-field theory, disorder can produce anomalous long time dynamics, where ensemble averages are dominated by the effect of rare sites.

## Funding statement

We acknowledge financial support from an ICAM-I2CAM postdoctoral fellowship (under the ICAM branches cost-sharing fund) for the execution of this work. C.C. gratefully acknowledges support from ESF (Intelbiomat Programme), R.F. from IP-SIQS and PRIN (Project 2010LLKJBX), J.K. from EPSRC programme ‘TOPNES’ (EP/I031014/1) and EPSRC grant (EP/G004714/2), and H.E.T. from NSF CAREER (DMR-1151810) and The Eric and Wendy Schmidt Transformative Technology Fund.

## Acknowledgements

We acknowledge helpful discussions with Peter Littlewood.

- Received April 22, 2014.
- Accepted June 17, 2014.

- © 2014 The Author(s) Published by the Royal Society. All rights reserved.