Fluorescence spectrum and thermalization in a driven coupled cavity array
FFluorescence spectrum and thermalization in a driven coupled cavity array
Dainius Kilda and Jonathan Keeling
SUPA, School of Physics and Astronomy, University of St Andrews, St Andrews, KY16 9SS, United Kingdom (Dated: December 6, 2018)We calculate the fluorescence spectra of a driven lattice of coupled cavities. To do this, weextend methods of evaluating two-time correlations in infinite lattices to open quantum systems;this allows access to momentum resolved fluorescence spectrum. We illustrate this for a driven-dissipative transverse field anisotropic XY model. By studying the fluctuation dissipation theorem,we find the emergence of a quasi-thermalized steady state with a temperature dependent on systemparameters; for blue detuned driving, we show this effective temperature is negative. In the lowexcitation density limit, we compare these numerical results to analytical spin-wave theory, providingan understanding of the form of the distribution function and the origin of quasi-thermalization.
By driving a system out of equilibrium, it is possibleto stabilize states of matter that are either not knownor are hard to achieve in thermal equilibrium. Classi-cally, driven systems have been extensively studied in theframework of pattern formation and dynamics [1]. Thestudy of quantum systems driven far from equilibriumis currently very active, in fields ranging from ultracoldatoms [2–5] to optically induced superconductivity [6],and hybrid matter-light systems [7]. One such class ofsystem is driven dissipative lattices [8–10]. This is mo-tivated by a variety of experimental platforms, includ-ing photonic crystal devices with quantum dots [11], mi-cropillar structures in semiconductor microcavities [12],trapped ions [13], and microwave cavities and supercon-ducting qubits [14, 15]. Depending on the combination ofcouplings and driving used, many different models can berealized, and for many of these models, driving and dis-sipation allow one to induce a wider variety of collectivestates than occur in thermal equilibrium [16–25].Most theoretical work on driven-dissipative lattices hasfocused on using order parameters or equal-time correla-tion functions to identify the phase diagram. For co-herent driving, such observables correspond to measur-ing the elastically scattered light. Less attention hasbeen paid to the properties of the incoherent fluores-cence from such lattices. From the quantum optics per-spective, incoherent fluorescence of a coherently drivensystem can reveal interactions and coherence times, asknown for the Mollow triplet fluorescence [26], which hasbeen seen in candidate systems for coupled cavity arrayssuch as quantum dots [27] and superconducting qubitscoupled to microwave cavities [28]. In extended systems,one can also access momentum-resolved spectra, e.g. bymeasuring the interference of light emitted from differentcavities. Moreover, second order correlations distinguishbunching or antibunching of photons — as studied the-oretically for a pair of coupled cavities [29, 30]. Appliedto extended systems, such measurements can make con-tact with quantities typically seen with condensed matterprobes such as angle resolved photon emission, spectro-scopic scanning tunneling microscopy, or neutron scat-tering. i.e. they measure the excitation and fluctuation spectrum of a correlated state, revealing the nature ofcorrelated states.There are other reasons to anticipate that calculationsof two-point and two-time correlations can provide un-derstanding beyond single-time observables. Firstly, forany correct treatment of a finite size system, symme-try breaking should not occur. This can also be truefor certain numerical approaches in infinite systems: un-less one uses the non-commuting limits of symmetry-breaking fields and system size, one finds a steady statedensity matrix with equal mixtures of symmetry-brokenstates [31, 32]. Two-time correlations allow one to in-stead ask how long symmetry breaking persists in re-sponse to a probe — i.e., long time correlations cor-respond to divergences of the zero frequency responseof a system. For driven systems, similar results maybe extended to the treatment of limit cycles and ‘non-equilibrium time crystals’ [24, 33, 34]. The density ma-trix, as an ensemble averaged quantity, involves averag-ing over the phase (or equivalently origin in time) of anylimit cycle, washing out any time dependence in the den-sity matrix. In contrast two-time correlations reveal suchcycles as a diverging response at non-zero frequency.Another motivation for studying two-time correlationsis to investigate thermalization. Thermalization in drivensystems has been studied in a number of contexts, includ-ing the ‘low energy effective temperature’ in the Keldyshfield theory of driven atom-photon systems [35–41] andthe mode populations in photon [42, 43] and polari-ton condensates [44–46]. This steady state behavior ina continuously driven system can also be connected tothe emergence of a prethermalized state following a sud-den quench in an isolated system [47–49] — in such aprethermalized state, there is a flow of energy betweendegrees of freedom at different scales. For a thermal-ized state, we expect the density matrix takes the Gibbsform ρ = exp ( − H eff /T eff ) with some effective Hamilto-nian H eff . One may note however that any density matrixcan be written in the Gibbs form; to make the criterionmeaningful one thus needs a method to independently de-termine H eff . This means simultaneously measuring theoccupations and densities of available states — this is the a r X i v : . [ c ond - m a t . o t h e r] D ec essence of the fluctuation dissipation theorem, which wediscuss below.In this letter, we find the two-time correlations of adriven dissipative lattice, and see the emergence of aquasi-thermalized state. We calculate both on-site andinter-site correlations, giving access to the momentumresolved fluorescence spectrum of a driven coupled cav-ity array. In order to eliminate boundary and finite-size effects, we work always with the translationally in-variant infinite lattice. On-site calculations in a finite-size lattice have also been recently studied for the XXZmodel [50]. While the methods we present are general,our work will focus on the transverse field anisotropic XYmodel (which has both the Ising and XY models as spe-cial cases), a driven dissipative realization of which wasproposed by Bardyn and ˙Imamo˘glu [51], and the steadystate properties studied [20, 52] using matrix productstate approaches. As shown in [51] and reviewed in thesupplementary material [53], this model can be realizedby an array of coupled cavities in the photon blockaderegime, with a two-photon pump that creates pairs ofphotons in adjacent sites (see Fig. 6). J κ J J J κ κ κ
FIG. 1. Coupled cavity array with hopping J , photon loss κ and two-photon pumping (blue line). When strong nonlinear-ity (purple shading) in each cavity leads to photon blockade,this yields the transverse field anisotropic XY model [51]. Following [20, 53], working in the rotating frameof the pump, the effective Hamiltonian has the form: H = − J (cid:80) j (cid:2) gσ zj + σ xj σ xj +1 + − ∆2 σ yj σ yj +1 (cid:3) . The di-mensionless transverse field g depends on the pump-cavity detuning, and the anistropy parameter ∆, givenby the ratio of pump strength and photon hopping J .For ∆ = 1 we recover the Ising model and for ∆ = 0the isotropic XY model. In the following we will work inunits of J. For the driven system, the Hamiltonian is ac-companied by photon loss at rate κ into empty radiationmodes. We thus have the master equation: ∂ t ρ = L{ ρ } = − i [ H, ρ ]+ κ (cid:88) j (cid:0) σ − j ρσ + j − σ + j σ − j ρ − ρ σ + j σ − j (cid:1) . (1)While a non-driven system would equilibrate with thebath, the time-dependent driving breaks detailed bal-ance and leads instead to a nonequilibrium steady state(NESS).The fluctuation and response spectra discussed aboverequire evaluating two-time correlation functions which, for a Markovian system, can be found using the quantumregression theorem [26]: (cid:68) O ( j )2 ( t ) O ( i )1 (0) (cid:69) = Tr (cid:104) O ( j )2 e t L O ( i )1 ρ ss (cid:105) , (2)where i, j label two lattice sites and 1 , ρ ss of the master equation (1). Wedo this by using the infinite Time Evolving Block Dec-imation (iTEBD) algorithm [54, 55] to find the trans-lationally invariant infinite MPS such that L { ρ ss } = 0.Starting from the NESS, we then calculate two-time cor-relations using Eq. (43). Because applying local oper-ators ˆ O to ρ ss breaks translational invariance, we canno longer propagate using iTEBD. For a finite size lat-tice, TEBD could be used, but this restricts the extent ofcorrelations in both space and time, as excitations are re-flected from the boundaries [56]. Fortunately, a methodto find such correlations in an infinite lattice has beendeveloped by Ba˜nuls et al. [57] for unitary evolution.This approach [57], which we extend to open systems,writes the time evolution between applying ˆ O and ˆ O as a tensor network, and contracting this network givesthe desired correlator (see [53] for details).Using this approach, we calculate the fluctuationspectrum S O,O † ( ω ) and the response function of thesystem χ (cid:48)(cid:48) O,O † ( ω ) which are at the heart of thefluctuation-dissipation theorem [26, 41], S O,O † ( ω ) = F ( ω ) χ (cid:48)(cid:48) O,O † ( ω ) , with the distribution function F ( ω ) dis-cussed below. Both S O,O † ( ω ) and χ (cid:48)(cid:48) O,O † ( ω ) are theFourier transforms of two-time correlations˜ S O,O † ( t ) = 12 (cid:68) { ˆ O ( t ) , ˆ O † (0) } (cid:69) , (3)˜ χ O,O † ( t ) = iθ ( t ) (cid:68) [ ˆ O ( t ) , ˆ O † (0)] (cid:69) , (4)which we may evaluate using Eq. (43).Figure 2 shows the on-site ( i = j ) fluctuation and re-sponse functions in frequency domain for ˆ O = ˆ O = ˆ O ∈{ σ x , σ z } and a range of values of transverse field g . Weshow both the Ising limit, (∆ = 1, left two columns) aswell as at small ∆ (right column), where analytic resultscan be found using spin-wave theory as discussed fur-ther below. The panels (a–c) show S ( ω ) which measuresthe occupations while, panels (d–f) show response func-tion χ (cid:48)(cid:48) ( ω ), which measures the density of states (DoS).We note that while at g = 0 , S ( ω ) for σ x ispeaked at ω = 0, its value always remains finite as thereis no phase transition in this open one-dimensional sys-tem [20, 25]. As we will discuss later, the form of thedensity of states seen here can be understood from themomentum resolved correlation functions.The bottom row of Fig. 2 shows the inverse distributionfunctions F ( ω ) − = χ (cid:48)(cid:48) O,O † ( ω ) /S O,O † ( ω ) for ˆ O = σ x , σ z respectively. In an equilibrium system, the distribu-tion function F ( ω ) depends only on whether ˆ O obeys (a) S ( ω ) XX, ∆ = 1.0 (b)ZZ, ∆ = 1.0 (c)XX, ∆ = 0.05 -0.60.0 (d) χ " ( ω ) -0.8-0.7-0.6-0.5-0.4-0.3-0.2-0.1 0 (e) -0.8-0.7-0.6-0.5-0.4-0.3-0.2-0.1 0 (f) -1.00.0 0 4 8 12 (g) F ( ω ) - ω -1-0.8-0.6-0.4-0.2 0
0 7 14 21 (h) ω -1-0.8-0.6-0.4-0.2 0
0 4 8 12 (i) ω g=1.0 g=5.0g=1.0, SW g=5.0, SWg=0.0 FIG. 2. Spectrum of fluctuations S ( ω ), imaginary part ofresponse function χ (cid:48)(cid:48) ( ω ), and inverse distribution function F ( ω ) − . Left two columns: Ising limit ∆ = 1, Right columnshows ∆ = 0 .
05 where spin-wave theory (solid lines) matcheswell. Energies given in units of J . Other parameters used: κ = 0 . Fermionic or Bosonic (anti-)commutation relations; forBosons it is: F ( ω ) ≡ n B ( ω ) + 1 = coth(( ω − µ ) / T ) . Ina driven dissipative system, F ( ω ) may take a more gen-eral form. However as identified in other contexts [7, 35–41], quasi-thermalisation of low energy modes often oc-curs, leading to the identification of a low energy effectivetemperature F ( ω ) ∼ T eff /ω . Note that since all calcula-tions are performed in the rotating frame, all frequenciesare measured relative to the pump frequency — i.e. thepump frequency acts as an effective chemical potential µ that sets the frequency at which F ( ω ) diverges.As seen in Fig. 2(g,h), F ( ω ) − is linear ω → | F ( ω ) | − ≤
1. At high frequencies the distributionfunction of a fully thermalised system asymptotically ap-proaches this value. In our non-equilibrium system wesee that in some cases the inverse distribution | F ( ω ) | − approaches 1 over a range of frequencies, however in allcases it falls falls below one at higher frequencies, in-dicating higher fluctuations than for a thermal state.The results shown give some indication that, at leastfor Fig. 2(g), the F ( ω ) approaches a thermal form moreclosely at larger g .The right column of Fig. 2 compares the MPS results(points) to analytic spin-wave theory [53, 58], which isvalid if the density of excitations is small. We see a goodagreement between spin-wave theory and MPS numer- -6-4-2 0 0 2 4 (a) T eff , ∆ = 1.0 T e ff g T xx T zz -6-4-2 0 0 3 6 9 (b) T eff , ∆ = 0.05 T e ff g T xx T yy -g/2 FIG. 3. The effective temperature T eff against transverse field g . We find T eff by fitting F ( ω ) (cid:39) A coth( bω ) for low fre-quencies ( ω ≤ . T eff ≡ A/ b . (a) MPS re-sults for σ x,z fluctuations and the transverse-field Ising limit(∆ = 1 . σ x,y fluctuations at∆ = 0 .
05. Energies given in units of J . Other parametersused: κ = 0 . ics at ∆ = 0 .
05 for σ x correlations (we do not show the σ z spectra for this ∆, as these vanish in the linearisedspin-wave theory). Remarkably, the agreement for F ( ω )is better than for S ( ω ) , χ (cid:48)(cid:48) ( ω ) individually. It is notablethat despite being a linear (i.e. non-interacting) theory,the spin-wave result reproduces both the low energy ef-fective temperature and the emergent plateau F ( ω ) (cid:39) k -dependent function F ( ω, k ) = (2 T eff ,k + λ k ω ) /ω ,with weighting by the k -dependent density of states [53].This form (which follows directly from the structure ofthe relevant linearised theory) leads directly to the exis-tence of a low energy effective temperature. The plateauat F ( ω ) (cid:39)
1, seen only at larger g , results from the lo-cal spectra averaging over many momentum states [53],however the form F ( ω, k ) inevitably leads to F ( ω ) ∝ ω at high frequencies, corresponding to the breakdown ofthe plateau.As well as the deviation from the thermal F ( ω ), a sec-ond distinction from an equilibrated system is that boththe distribution and the low-energy effective temperatureextracted differ depending on the system operator con-sidered. Figure 3(a) shows how T eff of σ x and σ z correla-tors vary with transverse field g . Fig. 3(b) shows similarresults for the spin-wave theory at small ∆ for σ x and σ y correlators (as noted above, σ z correlators vanish ina linearised theory). We observe that for ∆ → g → ∞ the σ x,y excitations thermalize to the same effective tem-perature, T eff ≈ − g/
2. This can be understood as T eff ,k becomes k independent in this limit, see [53].We only show results for g > S ( ω ), χ (cid:48)(cid:48) ( ω ), F ( ω ) − for values g and − g . This duality,discussed in [20] arises because a combination of g (cid:55)→ − g and a π rotation of the spin on every second site leadsto H (cid:55)→ − H . (A more general discussion of such dual-ities can be found in [59].) This duality means that onchanging the sign of g , the state of the system should cor-respond to reversing the sign of all energies. We may then -14-7 0 7 14 ω (a) S xx ( ω ,k), g=5.0 -20-10 0 10 20 (b) S zz ( ω ,k), g=5.0 - π - π /2 0 π /2 π k -8-4 0 4 8 ω (c) S xx ( ω ,k), g=1.0 - π - π /2 0 π /2 π k -10-5 0 5 10 (d) S zz ( ω ,k), g=1.0 FIG. 4. S ( ω, k ) momentum-resolved fluctuation spectrum forexcitations of: (a) σ x at g = 5 .
0; (b) σ z at g=5.0; (c) σ x at g = 1 .
0; (d) σ z at g = 1 .
0. Energies given in units of J . Otherparameters used: κ = 0 . note that fluctuation and dissipation spectra show differ-ent parity; χ (cid:48)(cid:48) ( − ω ) = − χ (cid:48)(cid:48) ( ω ), while S ( − ω ) = S ( ω ),and so F ( − ω ) = − F ( ω ). As such, the energy sign rever-sal under g (cid:55)→ − g yields a sign change of the distributionfunction and effective temperature. We find g < g > T eff < g is proportional to the pump-cavity detuning,so that g < g > g = 0, the susceptibility χ (cid:48)(cid:48) ( ω ) vanishes, so F ( ω ) − = 0 and the effective temper-ature diverges.So far we have evaluated correlation functions at equalpositions; this corresponds to recording all light from onecavity, which implicitly integrates over momentum. Moreinformation on the structure of the correlations is avail-able if we consider the momentum-resolved spectrum.This requires evaluating correlations at non-equal sites i, j , and performing a double Fourier transform with re-spect to separation in time t and space | i − j | . The re-sulting fluctuation spectra S ( ω, k ) are displayed in Fig. 4(The response function χ (cid:48)(cid:48) ( ω, k ) shows similar featuresas S ( ω, k )). We show the case for ˆ O = σ x , σ z , and twovalues of g (we consider only g >
0, since the dualitydiscussed above allows one to understand the effects of asign change of g ). All the features visible in these spectracan be described straightforwardly using excitation spec- tra derived from the Jordan-Wigner solution of H TFI (seee.g. [60] for details).At large positive g , the NESS is known [20] to be amaximum energy state with spins pointing in the − ˆ z direction, opposing the magnetic field. The spectrumof the σ x operator corresponds to single spin flips, sofollows the single particle dispersion ω ( k ) = (cid:15) ( k ) ≡ J (cid:112) g + 2 g cos( k ) (where we consider de -excitationsof the maximum energy state). This expression is shownby the black line in Fig. 4(a). In contrast, the σ z oper-ator corresponds to two-particle excitations, which comein two varieties. The first one is a two-particle contin-uum with ω = (cid:15) ( k ) + (cid:15) ( k ) , k + k = k . The envelopeof these states is given by (cid:15) min ( k ) < ω ( k ) < (cid:15) max ( k ) with (cid:15) max / min ( k ) = 4 J (cid:112) g ± g cos( k/ q to q + k , i.e. ω ( k ) = ∆ (cid:15) ( q, k ) ≡ (cid:15) ( q + k ) − (cid:15) ( q ). Thedominant contribution comes from q = 0, since this cor-responds to the maximum energy mode, which is max-imally occupied for a negative temperature state. Theblack solid line shows ∆ (cid:15) (0 , k ) which indeed matches thedominant feature observed. Given these momentum re-solved results, the momentum integrated spectral func-tions in Fig. 2(a,c) can be easily understood, with peaksarising from van Hove singularities at the band edges.Near g = 1 . σ x spectrum, Fig. 4(c) andFig. 2(a,c), the peaks at k = ± π become dominant. Inthe σ z spectrum, Fig. 4(d) and Fig. 2(b,d), the scatter-ing band and two-particle continuum overlap. The blacklines show the same expressions as discussed above. Aground state phase transition occurs for | g | <
1, hencethe gap closing at g = 1 .
0. In contrast, the NESS at g = 1 . g = 1 .
0. As one continues todecrease g → ω = 0, as seen in Fig. 2(a–d).In conclusion, we have calculated the two-time corre-lations of a driven-dissipative coupled-cavity array, pro-viding the fluorescence and absorption spectra. Due tothe duality between red and blue detuned scenarios, wefind that a blue pump-cavity detuning produces a quasi-thermalized state with a negative temperature. We havealso shown how the structure of F ( ω ) and emergent ther-malization can be understood using a spin-wave theory,and how momentum resolved fluorescence reveals thenature of quasiparticle excitations in the quasi-thermalstate. The system we have studied here is in the photonblockade regime, with at most one excitation per site.This restricts us to study “first-order” correlation func-tions. When generalizing to problems with a larger on-site Hilbert space, second-order photon counting correla-tions may also be of interest, in revealing the coherenceand statistics of any ordered state. Our results illustratehow calculating such correlations of the fluorescence canprovide new insights into the state of many-body drivendissipative systems.DK acknowledges support from the EPSRC CM-CDT(EP/L015110/1). JK acknowledges support from EP-SRC program TOPNES (EP/I031014/1). We are grate-ful to M. Hartmann and S. H. Simon for helpful discus-sions, and to M. Hartmann for helpful comments on themanuscript. [1] M. C. Cross and P. C. Hohenberg, Rev. Mod. Phys. ,851 (1993).[2] I. Bloch, J. Dalibard, and W. Zwerger, Rev. Mod. Phys. , 885 (2008).[3] H. Ritsch, P. Domokos, F. Brennecke, and T. Esslinger,Rev. Mod. Phys. , 553 (2013).[4] A. J. Daley, Adv. Phys. , 77 (2014).[5] T. Langen, R. Geiger, and J. Schmiedmayer, Annu. Rev.Condens. Matter Phys. , 201 (2015).[6] M. Mitrano, A. Cantaluppi, D. Nicoletti, S. Kaiser,A. Perucchi, S. Lupi, P. D. Pietro, D. Pontiroli, M. Ricc,S. R. Clark, D. Jaksch, and A. Cavalleri, Nature ,461 (2016).[7] I. Carusotto and C. Ciuti, Rev. Mod. Phys. , 299(2013).[8] S. Schmidt and J. Koch, Ann. Phys. (Berlin) , 395(2013).[9] C. Noh and D. G. Angelakis, Rep. Prog. Phys. , 016401(2016).[10] M. J. Hartmann, Journal of Optics , 104005 (2016).[11] A. Majumdar, A. Rundquist, M. Bajcsy, andJ. Vuˇckovi´c, Phys. Rev. B , 045315 (2012).[12] V. G. Sala, D. D. Solnyshkov, I. Carusotto, T. Jacqmin,A. Lemaˆıtre, H. Ter¸cas, A. Nalitov, M. Abbarchi, E. Ga-lopin, I. Sagnes, J. Bloch, G. Malpuech, and A. Amo,Phys. Rev. X , 011034 (2015).[13] K. Toyoda, Y. Matsuno, A. Noguchi, S. Haze, andS. Urabe, Phys. Rev. Lett. , 160501 (2013).[14] A. A. Houck, H. E. T¨ureci, and J. Koch, Nature Physics , 292 (2012).[15] M. Fitzpatrick, N. M. Sundaresan, A. C. Y. Li, J. Koch,and A. A. Houck, Phys. Rev. X , 011016 (2017).[16] I. Carusotto, D. Gerace, H. Tureci, S. De Liberato,C. Ciuti, and A. ˙Imamo˘glu, Phys. Rev. Lett. , 033601(2009).[17] M. J. Hartmann, Phys. Rev. Lett. , 113601 (2010).[18] T. Grujic, S. Clark, D. Jaksch, and D. Angelakis, NewJ. Phys. , 103025 (2012).[19] F. Nissen, S. Schmidt, M. Biondi, G. Blatter, H. E.T¨ureci, and J. Keeling, Phys. Rev. Lett. , 233603(2012).[20] C. Joshi, F. Nissen, and J. Keeling, Phys. Rev.A ,063835 (2013).[21] J. Jin, D. Rossini, R. Fazio, M. Leib, and M. J. Hart-mann, Phys. Rev. Lett. , 163605 (2013).[22] J. Jin, D. Rossini, M. Leib, M. J. Hartmann, andR. Fazio, Phys. Rev.A , 023827 (2014). [23] A. Biella, L. Mazza, I. Carusotto, D. Rossini, andR. Fazio, Phys. Rev.A , 053815 (2015).[24] M. Schir´o, C. Joshi, M. Bordyuh, R. Fazio, J. Keeling,and H. T¨ureci, Phys. Rev. Lett. , 143603 (2016).[25] T. E. Lee, S. Gopalakrishnan, and M. D. Lukin, Phys.Rev. Lett. , 257204 (2013).[26] H.-P. Breuer and F. Petruccione, The theory of openquantum systems (Oxford University Press, Oxford,2002).[27] A. N. Vamivakas, Y. Zhao, C.-Y. Lu, and M. Atat¨ure,Nature Physics , 198 (2009).[28] C. Lang, D. Bozyigit, C. Eichler, L. Steffen, J. M. Fink,A. A. Abdumalikov, M. Baur, S. Filipp, M. P. da Silva,A. Blais, and A. Wallraff, Phys. Rev. Lett. , 243601(2011).[29] T. C. H. Liew and V. Savona, Phys. Rev. Lett. ,183601 (2010).[30] M. Bamba, A. Imamo˘glu, I. Carusotto, and C. Ciuti,Phys. Rev. A , 021802 (2011).[31] S. R. K. Rodriguez, W. Casteels, F. Storme, N. Car-lon Zambon, I. Sagnes, L. Le Gratiet, E. Galopin,A. Lemaˆıtre, A. Amo, C. Ciuti, and J. Bloch, Phys.Rev. Lett. , 247402 (2017).[32] T. Fink, A. Schade, S. H¨ofling, C. Schneider, andA. ˙Imamo˘glu, Nature Physics , 365 (2018).[33] C.-K. Chan, T. E. Lee, and S. Gopalakrishnan, Phys.Rev.A , 051601 (2015).[34] F. Iemini, A. Russomanno, J. Keeling, M. Schir`o, M. Dal-monte, and R. Fazio, Phys. Rev. Lett. , 035301(2018).[35] S. Diehl, A. Micheli, A. Kantian, B. Kraus, H. B¨uchler,and P. Zoller, Nat. Phys. , 878 (2008).[36] S. Diehl, A. Tomadin, A. Micheli, R. Fazio, and P. Zoller,Phys. Rev. Lett. , 015702 (2010).[37] B. ¨Oztop, M. Bordyuh, ¨O. E. M¨ustecaplıo˘glu, and H. E.T¨ureci, New J. Phys. , 085011 (2012).[38] E. G. D. Torre, S. Diehl, M. D. Lukin, S. Sachdev, andP. Strack, Phys. Rev.A , 023831 (2013).[39] M. Buchhold, P. Strack, S. Sachdev, and S. Diehl, Phys.Rev.A , 063622 (2013).[40] L. Sieberer, S. Huber, E. Altman, and S. Diehl, Phys.Rev. Lett. , 195301 (2013).[41] L. M. Sieberer, M. Buchhold, and S. Diehl, Rep. Prog.Phys. , 096001 (2016).[42] J. Klaers, J. Schmitt, F. Vewinger, and M. Weitz, Nature , 545 (2010).[43] P. Kirton and J. Keeling, Phys. Rev.A , 033826 (2015).[44] T. D. Doan, H. T. Cao, D. B. Tran Thoai, and H. Haug,Phys. Rev. B , 085301 (2005).[45] J. Kasprzak, D. D. Solnyshkov, R. Andr´e, L. S. Dang,and G. Malpuech, Phys. Rev. Lett. , 146404 (2008).[46] Y. Sun, P. Wen, Y. Yoon, G. Liu, M. Steger, L. N. Pfeif-fer, K. West, D. W. Snoke, and K. A. Nelson, Phys. Rev.Lett. , 016602 (2017).[47] A. Polkovnikov, K. Sengupta, A. Silva, and M. Vengalat-tore, Rev. Mod. Phys. , 863 (2011).[48] J. Eisert, M. Friesdorf, and C. Gogolin, Nat. Phys. ,124 (2015).[49] T. Langen, T. Gasenzer, and J. Schmiedmayer, J. Stat.Mech. Theor. Exp. , 064009 (2016).[50] S. Wolff, J.-S. Bernier, D. Poletti, A. Sheikhan, andC. Kollath, “Evolution of two-time correlations in dis-sipative quantum spin systems: aging and hierarchical dynamics,” (2018), preprint, 1809.10464.[51] C.-E. Bardyn and A. ˙Imamo˘glu, Phys. Rev. Lett. ,253606 (2012).[52] E. Mascarenhas, H. Flayac, and V. Savona, Phys. Rev.A , 022116 (2015).[53] Supplementary Material containing derivation of the ef-fective model, details of the spin wave calculation, anda full description of the tensor network method for two-time correlations.[54] U. Schollw¨ock, Ann. Phys. , 96 (2011).[55] R. Or´us, Ann. Phys. , 117 (2014).[56] H. N. Phien, G. Vidal, and I. P. McCulloch, Phys. Rev.B , 245107 (2012).[57] M. C. Ba˜nuls, M. B. Hastings, F. Verstraete, and J. I.Cirac, Phys. Rev. Lett. , 240603 (2009).[58] D. C. Mattis, The theory of magnetism made simple:an introduction to physical concepts and to some usefulmathematical methods (World Scientific Publishing Com-pany, 2006).[59] A. C. Y. Li and J. Koch, New Journal of Physics ,115010 (2017).[60] S. Sachdev, Quantum phase transitions (Cambridge Uni-versity Press, Cambridge, 2007).[61] M. O. Scully and M. S. Zubairy,
Quantum Optics (Cam-bridge, 1997).[62] G. Ford and R. O’Connell, Phys. Rev. Lett. , 798(1996).[63] A. Altland and B. D. Simons, Condensed matter fieldtheory (Cambridge University Press, Cambridge, 2010).[64] E. Stoudenmire and S. R. White, New J. Phys. ,055026 (2010). SUPPLEMENTARY MATERIAL FOR: FLUORESCENCE SPECTRUM AND THERMALIZATION IN ADRIVEN COUPLED CAVITY ARRAYDRIVEN-DISSIPATIVE XY MODEL
This section provides the derivation of the effectivetransverse field anisotropic XY model which we study,starting from a model of a coupled cavity array, follow-ing Refs. [20, 51].We consider a 1D lattice of optical or microwave cav-ities supporting photon modes b j with tunneling ampli-tude J between adjacent cavities, and an on-site opticalnonlinearity U which induces effective photon-photon in-teractions in each cavity. Such a coupled cavity array isthus described by the Bose-Hubbard Hamiltonian: H = (cid:88) j (cid:104) ω c b † j b j + U b † j b † j b j b j − J (cid:16) b † j b j +1 + H.c. (cid:17)(cid:105) . In addition to these elements, we consider a two-photondrive Ω cos(2 ω P t ) near two-photon resonance ω P ≈ ω c .We then work in the limit of strong optical nonlinearitywith a perfect photon blockade, which restricts occupa-tions to at most one photon in each cavity. This strongnonlinearity then implies that the two-photon drive isonly resonant with creation of photon pairs on adjacentcavities. The above considerations allow us to replaceeach cavity mode with a spin-1 /
2, equivalent to replac-ing the bosonic operators by Pauli matrices: b j → σ − j .Our model Hamiltonian then becomes: H = (cid:88) j ω c σ zj − J (cid:88) j (cid:2) σ + j σ − j +1 + H.c. (cid:3) − Ω (cid:88) j (cid:2) σ + j σ + j +1 e − iω p t + H.c. (cid:3) . (5)If we then define the dimensionless parameters g = ( ω p − ω c ) / J , ∆ = Ω /J , we can transform H to a rotatingframe (at pump frequency ω p ) to gauge away the explicittime-dependence and write: H = − J (cid:88) j (cid:20) gσ zj + 1 + ∆2 σ xj σ xj +1 + 1 − ∆2 σ yj σ yj +1 (cid:21) . (6)The Hamiltonian H of a coupled cavity array thus takesa form of XY model where g acts as the transverse mag-netic field, and ∆ is the anisotropy of spin-spin interac-tions. The limit ∆ = 0 corresponds to the isotropic XYmodel and ∆ = 1 to the Ising model. SPIN-WAVE APPROXIMATION AT SMALLEXCITATION NUMBER
In this section we present further of the spin wave the-ory [58] used to describe the behavior at small excitationnumber. In particular, we can use this to understand ei-ther the limit of large | g | or small ∆, as both lead to smallexcitation number. Such an approach was used in Joshi et al. [20] for small ∆ to calculate static correlation func-tions; here we extend this to dynamical correlation func-tions and associated spectra.At ∆ = 0 (i.e. zero pumping Ω /J = 0), or at g → −∞ ,the NESS of our model corresponds to an empty state.For a small ∆, one can thus use spin-wave approxima-tion, which ignores the constraint on double occupancyof a lattice site, and so is only valid for a low density ofexcitations. In this small excitation number regime wecan revert from spin-1 / σ − j → b j , hence recovering aspects of aweakly interacting model. (Note that for large positive g a similar argument can be made, making use of theduality under g → − g discussed in the manuscript.) Calculating correlation functions
Spin-wave approximation and equations of motion
We first follow the steps described in [20] to derivethe Hamiltonian in terms of Bosonic system operators b k and b †− k . Working in the momentum basis, b k = (cid:80) j e ikj b j / √ N , the master equation, written as Eq. (1)in the main text becomes: ∂ t ρ = − i (cid:88) k [ h k , ρ ] + κ (cid:88) k (cid:16) b k ρb † k − b † k b k ρ − ρ b † k b k (cid:17) , (7)where h k = − (cid:0) b † k b − k (cid:1) (cid:18) g + cos( k ) ∆ cos( k )∆ cos( k ) g + cos( k ) (cid:19) (cid:18) b k b †− k (cid:19) , (8)and we have set J = 1, so all energies are measured inunits of J . Note that when ∆ controls the strength ofpair creation, while max( κ, g,
1) determines the cost ofcreating these excitations, so the small excitation regimecorresponds to ∆ (cid:28) max( κ, g, b k and b †− k can be written in a matrix form: ∂ t f ( t ) = M f ( t ) + v ( t ) , (9)with the vectors: f ( t ) = (cid:18) b k ( t ) b †− k ( t ) (cid:19) , v ( t ) = √ κ (cid:18) b in k ( t ) b † in − k ( t ) (cid:19) , (10)and the matrix: M = (cid:18) − κ + 2 i ( g + cos( k )) 2 i ∆ cos( k ) − i ∆ cos( k ) − κ − i ( g + cos( k )) (cid:19) . (11)Here, coupling to Markovian bath introduces the inputnoise term b in k ( t ). Since we consider a zero temperaturebath, there is only vacuum quantum noise, and the onlynonzero correlator is (cid:68) b in k ( t ) b † in k (cid:48) ( t (cid:48) ) (cid:69) = δ k,k (cid:48) δ ( t − t (cid:48) ).The solution of (9) is: f ( t ) = e Mt f (0) + (cid:90) t dt (cid:48) e M ( t − t (cid:48) ) v ( t (cid:48) ) . In the long-time limit t → ∞ we find the expressions forsystem operators: b k ( t ) = √ κ (cid:90) t dt (cid:48) (cid:104) G ( t − t (cid:48) ) b in k ( t (cid:48) ) + G ( t − t (cid:48) ) b † in − k ( t (cid:48) ) (cid:105) ,b †− k ( t ) = √ κ (cid:90) t dt (cid:48) (cid:104) G ∗ ( t − t (cid:48) ) b † in − k ( t (cid:48) ) + G ∗ ( t − t (cid:48) ) b in k ( t (cid:48) ) (cid:105) . (12)where the propagators G , ( τ ) are matrix elements of e Mt given by: G ( τ ) = e − κτ (cid:20) cos( ξ k τ ) + i(cid:15) k sin( ξ k τ ) ξ k (cid:21) , (13) G ( τ ) = iη k e − κτ sin( ξ k τ ) ξ k , (14)with dispersions (cid:15) k = 2( g + cos( k )), η k = 2∆ cos( k ), and ξ k = (cid:112) (cid:15) k − η k . Correlations and effective temperatures for ˆ σ x After deriving the system operators, we now proceedto calculate the frequency-resolved spectra for XX cor-relations. Since σ xj → b j + b † j in the spin-wave limit, wecan express the on-site XX two-time correlator as:˜ C xx ( τ ) = (cid:104) σ x (0) σ x ( τ ) (cid:105) = (cid:10) b † (0) b † ( τ ) (cid:11) + (cid:104) b (0) b ( τ ) (cid:105) + (cid:10) b † (0) b ( τ ) (cid:11) + (cid:10) b (0) b † ( τ ) (cid:11) , (15)where the correlations are given by a Fourier transformfrom momentum to real space: (cid:10) b † (0) b ( τ ) (cid:11) = (cid:90) π − π dk e ikl (cid:68) b † k (0) b k ( τ ) (cid:69) (cid:12)(cid:12)(cid:12)(cid:12) l =0 , and similar expressions for other correlators. We thensubstitute in the solutions for operators (12), and evalu-ate the time integrals at unequal times, t + τ and t . Thisgives a two-time correlator:˜ C xx ( τ ) = e − κτ π (cid:90) π − π dk (cid:34) cos( ξ k τ ) + i ( η k − (cid:15) k ) sin( ξ k τ ) ξ k + η k ( η k − (cid:15) k ) ξ k + κ (cid:18) cos( ξ k τ ) + κ sin( ξ k τ ) ξ k (cid:19)(cid:35) . (16)The quantities of interest are the fluctuation spec-trum S ( ω ) and susceptibility χ (cid:48)(cid:48) ( ω ) given by the Fouriertransforms of ˜ S ( τ ) = ( ˜ C ( τ ) ∗ + ˜ C ( τ )) and ˜ χ ( τ ) = i Θ( τ )( ˜ C ( τ ) ∗ − ˜ C ( τ )) respectively. Plugging in (16) andtaking a Fourier transform with respect to τ , we obtain XX spectra: S xx ( ω ) = κπ (cid:90) π − π dk P k + ω + 2 η k ( η k − (cid:15) k ) Q k ( ω ) , (17) χ (cid:48)(cid:48) xx ( ω ) = 2 κωπ (cid:90) π − π dk η k − (cid:15) k Q − k ( ω ) , (18)where we have introduced auxiliary functions P k = ξ k + κ , and Q k ( ω ) = ( P k − ω ) + (2 ωκ ) .One can then substitute e ik → z , and the integrals in(17), (18) become contour integrals around a unit circle C with | z | = 1. The values of (17), (18) are then determinedby the residues of poles z = Z located inside C (i.e. with | Z | < Z = θ s ,s ± (cid:113) θ s ,s − ,θ s ,s = − g + s (cid:113) g ∆ − − ∆ [ κ − ω + s iωκ ]1 − ∆ , with s = ± s = ±
1. Evaluating the contour integralsgives fluctuation spectrum: S xx ( ω ) = 2 κ (cid:88) | Z n | < Z n α n , (19)where α n = (cid:2) (1 − ∆) Z n + 2 gZ n + (1 − ∆) (cid:3) + ( ω + κ ) Z n (1 − ∆ ) (cid:81) m =1 ,m (cid:54) = n ( Z n − Z m ) . Similarly the susceptibility is given by: χ (cid:48)(cid:48) xx ( ω ) = − κω (cid:88) | Z n | < Z n β n , (20)where β n = (1 − ∆) Z n + 2 gZ n + (1 − ∆)(1 − ∆ ) (cid:81) m =1 ,m (cid:54) = n ( Z n − Z m ) . In both expressions, the sum runs over poles Z n insidethe unit circle C with | Z n | <
1. From (19) and (20)it is straightforward to derive the distribution function F xx ( ω ) of fluctuation-dissipation theorem: F xx ( ω ) = S xx ( ω ) χ (cid:48)(cid:48) xx ( ω ) = − ω (cid:80) | Z n | < Z n α n (cid:80) | Z n | < Z n β n . (21)It is evident that in the low frequency limit ω → F xx ( ω ) is dominated by the 1 /ω diver-gence. The effective thermalization of NESS thus alreadyemerges in spin-wave theory, leading to an effective tem-perature T eff : T eff ,xx = − (cid:80) | Z n | < Z n α n (cid:80) | Z n | < Z n β n (cid:12)(cid:12)(cid:12)(cid:12) ω =0 . (22) Correlations and effective temperatures for ˆ σ y Similarly to the previous subsection, we can calculatethe fluctuation spectrum and susceptibility for
Y Y exci-tations in the spin-wave limit where σ yj → − i ( b j − b † j ).By following the same steps as for XX correlators, weobtain the Y Y spectra: S yy ( ω ) = κπ (cid:90) π − π dk Q − k ( ω ) (cid:2) P k + ω + 2 η k ( η k + (cid:15) k ) (cid:3) , (23) χ (cid:48)(cid:48) yy ( ω ) = 2 κωπ (cid:90) π − π dk Q − k ( ω )( η k + (cid:15) k ) . (24)The only differ from the expressions for XX spectra (17),(18) by ( η k − (cid:15) k ) → ( η k + (cid:15) k ). Note that as a result,the contours integrals to evaluate have the same poles asfor XX correlators, but with different residues and thusdifferent weights.Continuing along the same steps as for XX correlators,we find the Y Y fluctuation spectrum: S yy ( ω ) = 2 κ (cid:88) | Z n | < Z n γ n , (25)where γ n = (cid:2) (1 + ∆) Z n + 2 gZ n + (1 + ∆) (cid:3) + ( ω + κ ) Z n (1 − ∆ ) (cid:81) m =1 ,m (cid:54) = n ( Z n − Z m ) . The
Y Y susceptibility is given by: χ (cid:48)(cid:48) yy ( ω ) = − κω (cid:88) | Z n | < Z n δ n , (26)where δ n = (1 + ∆) Z n + 2 gZ n + (1 + ∆)(1 − ∆ ) (cid:81) m =1 ,m (cid:54) = n ( Z n − Z m ) . and the poles Z n are the same as in XX spectra. One canthen straightforwardly derive the distribution function F yy ( ω ): F yy ( ω ) = S yy ( ω ) χ (cid:48)(cid:48) yy ( ω ) = − ω (cid:80) | Z n | < Z n γ n (cid:80) | Z n | < Z n δ n . (27)and the effective temperature T eff : T eff ,yy = − (cid:80) | Z n | < Z n γ n (cid:80) | Z n | < Z n δ n (cid:12)(cid:12)(cid:12)(cid:12) ω =0 . (28) Vanishing correlations for σ z Note that while the above approach allows calculationof the XX and Y Y correlators, we cannot use this smallexcitation number approximation to find the ZZ correla-tors. To see this, note that if we express the ZZ two-timecorrelator using σ zj = b † j b j − b j b † j , then in the spin-wavelimit:˜ C zz ( τ ) = (cid:104) σ z (0) σ z ( τ ) (cid:105) == (cid:10) b † (0) b (0) b † ( τ ) b ( τ ) (cid:11) − (cid:10) b (0) b † (0) b † ( τ ) b ( τ ) (cid:11) + (cid:10) b (0) b † (0) b ( τ ) b † ( τ ) (cid:11) − (cid:10) b † (0) b (0) b ( τ ) b † ( τ ) (cid:11) . (29)Since the problem involves non-interacting bosons, thesteady state is Gaussian and we can expand the four-field correlators using Wick’s theorem, which leads to˜ C zz ( τ ) = 0. Thus, at this order of approximation all ZZ spectra are trivially zero. This is to be expected, since σ z correlations are quartic and the spin wave theory islinear. Temperature of individual bosonic modes
Building on the spin-wave theory introduced above,we next discuss thermalization and the appearance of aneffective temperature from this linear theory. To under-stand these effects, it is helpful to consider the contri-bution of individual bosonic modes. From (17), (18) wethus define the momentum-resolved fluctuation spectrumand susceptibility: S xx ( ω, k ) = κπ Q − k ( ω ) (cid:2) P k + ω + 2 η k ( η k − (cid:15) k ) (cid:3) . (30) χ (cid:48)(cid:48) xx ( ω, k ) = 2 κωπ Q − k ( ω )( η k − (cid:15) k ) . (31)Then the distribution function of an individual bosonic k -mode: F xx ( ω, k ) = S xx ( ω, k ) χ (cid:48)(cid:48) xx ( ω, k ) = ξ k + ω + κ + 2 η k ( η k − (cid:15) k )2 ω ( η k − (cid:15) k ) . (32)0Focusing on the frequency dependence of this expression,we see it can be written in the form: F xx ( ω, k ) = 2 T eff ,xx,k + λ xx,k ω ω , (33)where λ xx,k = [2( η k − (cid:15) k )] − , and the effective temper-ature of an individual bosonic mode, defined from the ω → T eff ,xx,k = κ + ( η k − (cid:15) k ) η k − (cid:15) k ) , (34)where we have used the definition of ξ k = (cid:15) k − η k .From the above, we see that despite considering a lin-earized (i.e. non-interacting) theory, a low energy ef-fective temperature emerges for each individual k mode.However, the functional form of F xx ( ω, k ) does not showa plateau around F xx ( ω, k ) (cid:39)
1, as expected for an equi-librium system, and as sometimes seen from the MPSnumerics for F xx ( ω ). Instead | F xx ( ω, k ) | − shows a peakat ω = ω ∗ xx,k ≡ (cid:112) T eff ,xx,k /λ xx,k = (cid:112) κ + ( η k − (cid:15) k ) ,with a peak height | F xx ( ω ∗ xx,k , k ) | − = 1 (cid:112) T eff ,xx,k λ xx,k = (cid:115) ( η k − (cid:15) k ) κ + ( η k − (cid:15) k ) . One may see that as required, this peak value is alwaysless than one, and approaches one if damping is weakcompared the energy difference | η k − (cid:15) k | . While this mo-mentum resolved distribution function does not show aplateau, as seen in Fig. 2, a plateau does arise for the site-local (i.e. momentum integrated) result. This site-localdistribution function can be expressed as a weighted aver-age of distributions F ( ω, k ) of individual bosonic modes: F ( ω ) = (cid:82) π − π dk F ( ω, k ) χ (cid:48)(cid:48) ( ω, k ) (cid:82) π − π dk χ (cid:48)(cid:48) ( ω, k ) . (35)Because the location of the peak for each mode k differs, this weighted average shows a plateau aris-ing from combining all these peaks, stretching overa range of frequencies, set by the range of peakfrequencies, i.e. (cid:112) κ + 4[ g − (1 − ∆)] < ω ∗ xx,k < (cid:112) κ + 4[ g − (1 + ∆)] . Note that as such the locationof the plateau moves to higher frequencies as we increase g , as seen in Fig. (2) of the ma text. Note also that theplateau is always finite, and F ( ω ) ∝ ω at large enoughfrequency.As noted earlier, if we look at Y Y correlations in placeof XX , the only change is to replace η k − (cid:15) k → η k + (cid:15) k in the above expressions, including in the density ofstates χ (cid:48)(cid:48) ( ω, k ). Because of this change, one has both that T eff ,yy,k differs from T eff ,xx,k , as well as the distributionof occupied modes changing.In the limit of large g , Eq. (34) becomes T eff ,xx,k ≈ T eff ,yy,k ≈ − g/
2, independent of momentum k and the operator being measured. Since this result is indepen-dent of momentum, the local effective temperature fromEq. (35) also approaches this value, T eff ≈ − g/ . (36)In Fig.(3) of the main text we show the extracted tem-perature of the spin-wave theory for both σ x and σ y cor-relators: at large g , where quasi-thermalization holds,we see both temperatures approach this same value.In the opposite limit, of small g , one may note that η k ± (cid:15) k = 2(∆ ±
1) cos( k ) ± g can now pass through zerofor some real k , leading to a divergence of both T eff ,k and λ k . This divergence is however integrable, giving a finiteform of F ( ω ) and the corresponding T eff , except as seenat g = 0 in Fig.(3). Correlations for the Ising limit
As noted earlier, the small excitation density limit canbe understood as resulting either from small ∆ or large g . As such, this approximation should remain valid evenfor the Ising limit, ∆ = 1, as long as g (cid:29)
1. However,at first appearance the above results are singular in thelimit ∆ = 1. This is in fact not the case as one mayreadily check. For example, considering ˆ O = ˆ σ x , we seethat as ∆ → Z are still given by θ s ,s but the values of θ s ,s become singular, specifically: θ s ,s = (cid:40) − ( κ + s iω ) / g s = +1 − g/ (1 − ∆) s = − . Inserting this into the definition of Z , we find that of theeight poles, four remain finite, while two tend to zero andtwo to infinity. Pairs of poles at zero and infinity in factcancel, as long as the residue at the remaining poles isfinite. This can be checked to be true, with the residuesbecoming α n = ( ω + κ ) / g for the finite poles. Thebehavior in the limit ∆ → (a) S ( ω ) ω TN, g=1.0SW, g=1.0 TN, g=5.0SW, g=5.0 -0.60.0 0 4 8 12 (b) χ " ( ω ) ω -1.00.0 0 4 8 12 (c) F ( ω ) - ω FIG. 5. Correlation functions for ∆ = 1, comparing MPSnumerics (points) with spin wave calculations (lines). Panels(a-c) show the spectrum of fluctuations S ( ω ), imaginary partof response function χ (cid:48)(cid:48) ( ω ), and the inverse distribution func-tion F ( ω ) − respectively. The results match well for g = 5 . g = 1 . As seen in Fig. 5, the spin wave theory indeed accu-rately captures the behavior at large g , but clearly failsin the case g = 1 .
0, where it would not be expected to1hold. As with the small ∆ results in the main text, wesee that the distribution function matches more accu-rately than the fluctuation and response functions sepa-rately. We may note that for ∆ = 1, there is never a trueplateau in the distribution function. This can be under-stood from the fact that for ∆ = 1, (cid:15) k − η k = 2 g , so both T eff ,xx,k = − ( κ +4 g ) / g and λ xx,k = − / g become in-dependent of k , meaning that F xx ( ω, k ) is also indepen-dent of k , so the integrated version F ( ω ) follows the formof Eq. (33), having a peak at ω ∗ = (cid:112) κ + 4 g as is visiblein Fig. 5. Note however that for Y Y correlations the samestatement would not be true, as (cid:15) k + η k = 2 g + 4 cos( k )is k -dependent. Fluctuation-dissipation relation in linear theories
It is notable that, as discussed above, the linearized thespin-wave result predicts a low energy quasi-thermal dis-tribution with a non-zero effective temperature for XX and Y Y excitations. This is particularly notable in thelight of papers, e.g. Ford and O’Connell [62] which sug-gest that the fluctuation dissipation theorem fails for aMarkovian dissipation of a bath, and one would expect tofind the effective F ( ω ) to be frequency independent, cor-responding approximately to zero temperature. This sec-tion discusses why the model we consider does not showsuch behavior. Specifically, as already noted, the struc-ture of distribution function is dependent on the modeconsidered. If in place of the XX and Y Y correlationswe had directly considered correlations of the anihilationoperators, ˜ C bb † ( τ ) = (cid:10) b (0) b † ( τ ) (cid:11) , then we would havefound the distribution function F bb † ( ω ) would be flat.We first discuss this point, and why it is that XX and Y Y distribution functions are not flat. Since Ref. [62]also considers the analogue of XX correlations, we thenaddress further differences between the model discussedthere and our results.We first discuss the observation that the distributionfunction for correlators of annihilation and creation op-erators generally leads to a flat distribution. A quantumharmonic oscillator (with frequency Ω, and field oper-ators b , b † ) interacting with a bath of radiation modes(with frequency ω k , field operators B † k , B k for each mode)via coupling strength is described by the Hamiltonian H = Ω b † b + (cid:88) k ω k B † k B k + (cid:88) k (cid:104) g k bB † k + H.c. (cid:105) . (37)Here g k is the system-bath coupling strength. Using thestandard input-output formalism [61] in the Markovianlimit we derive Heisenberg-Langevin equation of motion ∂ t b ( t ) = − i Ω b ( t ) − κb ( t ) + √ κb in ( t ) , (38)where the input noise operator b in ( t ) is introduced bycoupling to Markovian bath. For a zero-temperature bath there is only vacuum noise and the only non-zerocorrelator is (cid:10) b in ( t ) b † in ( t (cid:48) ) (cid:11) = δ ( t − t (cid:48) ). Next, one canobtain the steady-state solution (at t → ∞ ): b ( t ) = √ κ (cid:90) t −∞ dt (cid:48) e − i Ω( t − t (cid:48) ) − κ ( t − t (cid:48) ) b in ( t ) . (39)The only non-vanishing two-time correlator then is (cid:104) b (0) b † ( τ ) (cid:105) = e − κ | τ | + i Ω τ . (40)Taking a Fourier transform of the symmetrized cor-relator ˜ S bb † ( τ ) = (cid:10) { b (0) , b † ( τ ) } (cid:11) = e − κ | τ | + i Ω τ andresponse function ˜ χ bb † ( τ ) = iθ ( τ ) (cid:10)(cid:2) b (0) , b † ( τ ) (cid:3)(cid:11) = iθ ( τ ) e − κ | τ | + i Ω τ gives the fluctuation spectrum and sus-ceptibility: S bb † ( ω ) = χ (cid:48)(cid:48) bb † ( ω ) = 2 κ ( ω − Ω) + κ . (41)Subsequently, one obtains a flat distribution spectrum for b , b † modes, in contrast to the quasi-thermal distributionof the XX and Y Y modes: F bb † ( ω ) = S bb † ( ω ) χ (cid:48)(cid:48) bb † ( ω ) = 1 . (42)Our spin-wave equations differs from the above deriva-tion in that the spin-wave theory has anomalous termsproportional to ∆. However, the derivation above canbe extended just as well to a linear theory with anoma-lous terms since its Hamiltonian can be diagonalized eas-ily using the Bogoliubov transformation [63]. The cru-cial difference that occurs in the spin wave theory dis-cussed above is our calculation of XX and Y Y correla-tions, which mean fluctuation and dissipation terms in-volve sums and differences of correlators (cid:104) b (0) b † ( τ ) (cid:105) and (cid:104) b † (0) b ( τ ) (cid:105) . Once these are both considered, the singlemode functions results in Eq. (30–32) follow, giving afrequency dependent result.As noted earlier, the above result is notable in con-nection to the argument by Ford and O’Connell [62]that quantum regression can never give a thermal spec-trum. The problem considered there is similar to oursin that there are anomalous terms (since no rotating-wave approximation is made in the system bath cou-pling), and the correlations considered are the XX cor-relations. However, there is a crucial difference in thatRef. [62] considers Ohmic rather than Markovian dissi-pation. This difference leads to the different conclusionsof the previous section. TENSOR NETWORK APPROACH FORTWO-TIME CORRELATIONS
In this section, we describe the tensor network methodfor computing two-time correlations in open quantum2
T(O ) T Λ (2) Γ (1) Γ (2) Λ (1) ... Λ (2) ... O O ... ... δ t/2 δ t R j T(O ) L i δ t x t δ t δ t/2 FIG. 6. Tensor network used to evaluate two-time correla-tions, using boundary eigenvectors (orange) for calculatingtwo-time correlations under non-unitary Liouvillian propaga-tor. systems in the thermodynamic limit. We use quantum re-gression to calculate two-time correlations Eq. (43), start-ing from the NESS density matrix ρ SS , represented bya translationally invariant infinite MPS that was previ-ously computed using infinite TEBD algorithm [54, 55]. (cid:68) O ( j )2 ( t ) O ( i )1 (0) (cid:69) = Tr (cid:104) O ( j )2 e t L O ( i )1 ρ ss (cid:105) . (43)For a finite size lattice, one could still directly use TEBDalgorithm to perform the time evolution in (43). How-ever, such direct propagation is incompatible with infiniteTEBD: application of a local operator O to ρ SS breakstranslational invariance. A naive solution would be touse a finite size extrapolation, which is prone to bound-ary and finite-size effects. In particular, the finite latticesize would restrict the extent of correlations in both spaceand time, as excitations will be reflected back from theboundaries, and the simulation will be no longer valid atlater times [56]. Such simulation would also inccur anadditional computational cost that scales linearly withthe system size which will be inefficient for large latticesneeded to approximate the thermodynamic limit, in com-parison to only two sites required in the infinite TEBD.Nonetheless, such a method has been recently used tocalculate aging dynamics in the XXZ model [50].Fortunately, it is possible to avoid these issues thatarise due to finite size altogether. A method to computetwo-time correlations in an infinite system directly (i.e.without resorting to a finite size extrapolation) has beenproposed by Ba˜nuls et al [57] for unitary evolution inisolated systems. In our work, we extend this approach to open quantum systems whose dynamics is governed bythe quantum master equation for density matrices. Weprovide a brief description of the algorithm below.The idea is to construct a network representing the en-tire time evolution, instead of evolving the MPS in timestep by step. We start from an infinite MPS represent-ing vectorized NESS density matrix | ρ ss (cid:105) and apply thefirst operator O ( i )1 (0) at the initial time t = 0. Then, forevery time evolution step we insert a propagator MPO.After repeating this for the required number of time stepswe apply the second operator O ( j )2 ( t ) at the final time t .Taking the trace at the final time removes the danglingphysical dimension at each site of the last MPO propaga-tor. This procedure produces a 2D tensor network thatis infinite along the spatial axis but finite along the timeaxis, giving an unnormalized two-time two-point correla-tor: (cid:68) O ( j )2 ( t ) O ( i )1 (0) (cid:69) ∝ T N →∞ T O T | i − j |− T O T N →∞ , (44)where T = is a transfer matrix of the evolved densitymatrix, and T O , are transfer matrices containing an ap-plication of operators O , at the initial and final timesat lattice sites i , j , as shown in Fig. 1(b) of the maintext.Since the network is translationally invariant, T is thesame on every site (except at the sites where the op-erators are applied) and we may effectively replace thesemi-infinite lattices, to the left and right of the siteswhere ˆ O and ˆ O act, by the left and right eigenvec-tors of T corresponding to its largest eigenvalue λ sincelim N →∞ T N = λ N | R (cid:105)(cid:104) L | . In practice, we compute anMPS approximation to the eigenvectors | R (cid:105) , (cid:104) L | by us-ing the MPS-MPO power method. We multiply an ini-tial arbitrary MPS (oriented along the time axis) by T (represented as an MPO along the time axis) a suffi-cient number of times until it converges to | R (cid:105) , (cid:104) L | forthe right- and left-multiplication respectively. We trun-cate the MPS bonds after each multiplication using themethod described in [64], performing truncation alongthe time axis. Once we have calculated | R (cid:105) and (cid:104) L | ,the resulting network is finite along both space and timeaxes. It can then be easily contracted using MPO-MPSand MPS-MPS multiplications [54, 64] to give any two-time two-point correlator: (cid:68) O ( j )2 ( t ) O ( i )1 (0) (cid:69) = (cid:10) L | T O T | i − j |− T O | R (cid:11) λ | i − j | +1 , (45)normalized by trace Tr( ρ ) = (cid:10) L | T | i − j | +1 | R (cid:11) = λ | i − j | +1+1