Dynamics of strongly coupled disordered dissipative spin-boson systems
Eliana Fiorelli, Pietro Rotondo, Federico Carollo, Matteo Marcuzzi, Igor Lesanovsky
DDynamics of strongly coupled disordered dissipative spin-boson systems
Eliana Fiorelli, Pietro Rotondo, Federico Carollo, Matteo Marcuzzi and Igor Lesanovsky
School of Physics and Astronomy, University of Nottingham, Nottingham, NG7 2RD, UK andCentre for the Mathematics and Theoretical Physics of Quantum Non-equilibrium Systems,University of Nottingham, Nottingham NG7 2RD, UK (Dated: May 1, 2019)Spin-boson Hamiltonians are an effective description for numerous quantum many-body systemssuch as atoms coupled to cavity modes, quantum electrodynamics in circuits and trapped ion sys-tems. While reaching the limit of strong coupling is possible in current experiments, the under-standing of the physics in this parameter regime remains a challenge, especially when disorder anddissipation are taken into account. Here we investigate a regime where the many-body spin dynam-ics can be related to a Ising energy function defined in terms of the spin-boson couplings. While inthe coherent weak coupling regime it is known that an effective description in terms of spin Hamil-tonian is possible, we show that a similar viewpoint can be adopted in the presence of dissipationand strong couplings. The resulting many-body dynamics features approximately thermal regimes,separated by out-of-equilibrium ones in which detailed balance is broken. Moreover, we show thatunder appropriately chosen conditions one can even achieve cooling of the spin degrees of freedom.This points towards the possibility of using strongly coupled dissipative spin-boson systems forengineering complex energy landscapes together with an appropriate cooling dynamics.
Introduction—
Prominent platforms for quantum sim-ulation, such as cavity, circuit [1] or waveguide quantumelectrodynamics [2] as well as trapped ions [3, 4] can bemodeled by ensembles of two-level systems interactingvia bosonic degrees of freedom (electromagnetic modes orphonons). While the weak coupling regime is relativelywell understood and can be treated by a perturbative in-tegration of the bosonic degrees of freedom, the strongcoupling limit is far more challenging [5].An additional layer of complexity is added by the pres-ence of disorder, i.e. when individual spins couple to thebosonic “environment” at different strengths. Such a set-ting is relevant for at least two reasons. First, some de-gree of quenched disorder may always be present in real-istic systems and, second, one may engineer non-uniformcouplings for practical applications: systems with tun-able quasi-random couplings often form the basis for aphysical implementation of complex optimization prob-lems, which may for instance be solved via quantum an-nealing protocols [6, 7].Disordered spin-boson systems have only recentlymoved into the focus of theoretical investigations. Ref-erences [8, 9] explore the emergence of glassiness whenmany electromagnetic modes interact with an ensembleof qubits. In Refs. [10, 11], instead, spin-glass tech-niques are employed to show that the same system effec-tively realizes an associative memory. Most of these tech-niques, however, cannot be straightforwardly generalizedto study open quantum dynamics in the strong couplingregime, and only a few studies deal with disordered openquantum systems [12–14]. This topic acquires furtherrelevance in the light of recent experimental progress inmultimodal cavity QED, which realize tunable range [15]and sign-changing [16] photon-mediated atomic interac-tions.In this work we investigate a disordered and dissipa-
FIG. 1.
Dissipative spin-boson system. (a) N weaklydriven (at strength Ω) two-level systems (spins) are stronglycoupled to a single bosonic mode with couplings g k ( k = , . . . , N ). Gain and loss of the bosons occur at rates κ and γ , respectively. (b) The resulting effective dissipative dynam-ics of the spins is related to a fully-connected Ising energyfunction ( σ k = ± E (⃗ σ ) [Eq. (4)], in which the interactionstrength between spins i and j is proportional to g i g j . Theeffective dynamics features regimes which permit cooling ofthe many-body spin state, i.e. significant population of thelow-energy configurations. tive system in which weakly driven spins are stronglycoupled to a bosonic mode (see Fig. 1a). We employ aperturbative approach which relies on the weakness ofthe driving rather than of the spin-boson coupling. Wefind that the effective spin dynamics is governed by a rateequation that depends on a fully-connected Ising energyfunction as sketched in Fig. 1b. Depending on the rates ofbosonic loss and gain we identify several distinct dynam- a r X i v : . [ c ond - m a t . s t a t - m ec h ] A p r ical regimes: two of them are high-temperature ones, inwhich the stationary state of the system is fully mixed. Afurther one mimics an effective low-temperature dynam-ics, which permits cooling of the spin system. Outsidethese the dynamics is generally non-thermal and detailedbalance is broken. This link between an open, stronglycoupled spin-boson system and the physics of disorderedIsing spin systems opens up the possibility of engineeringcomplex classical energy landscapes — with importancein the context of optimization problems [17] or associa-tive memories [18] — together with a cooling protocol. Model —
We consider an ensemble of N two-level sys-tems interacting with a single bosonic mode described bythe following Dicke Hamiltonian [19–22]:ˆ H = ω ˆ a † ˆ a + N ∑ i = g i ˆ σ zi ( ˆ a † + ˆ a ) + Ω N ∑ i = ˆ σ xi . (1)Here, ˆ σ x,y,zi are the Pauli operators and ˆ a and ˆ a † thebosonic annihilation and creation operators. The param-eters ω and Ω denote the fundamental frequency of thebosons and the coherent coupling strength between thetwo spin states, respectively. The spin-boson couplings g i are assumed to be independent and randomly distributedwith zero mean and variance g .We include dissipation on the boson in the form ofMarkovian gain and loss processes. The density matrix ρ of this open quantum system therefore evolves under aLindblad equation˙ ρ = L ρ = − i [ ˆ H, ρ ] + ∑ n = l,g ˆ L n ρ ˆ L † n − { ˆ L † n ˆ L n , ρ } (2)with the jump operators ˆ L l = √ γ ˆ a , ˆ L g = √ κ ˆ a † where γ ( κ ) is the loss (gain) rate and γ > κ ≥ N ions coupling to the centre-of-massphonon mode. As it has been shown for the quantumRabi model [25] and eventually generalized to the Dickemodel [24], the application of multiple laser fields on theions yields both the spin dependent coupling g k ˆ σ zk ( ˆ a + ˆ a † ) and the weak driving term, Ωˆ σ xk entering Eq. (1). Finally,as illustrated in Fig. 1a, the gain and loss dynamics canbe achieved by applying lasers on the ions on the edge ofthe chain, which is discussed in Ref. [26]. Since this ionis coupled to the same phonon mode as the other ionsthis effectively implements jump operators of the formintroduced in Eq. (2). Spin dynamics at strong coupling —
We explore thedynamics (2) in the strong coupling regime, i.e. whenthe driving acting on the spins is much weaker thanthe spin-boson interaction (Ω ≪ g ). In the following,we sketch the perturbative technique we employ for thispurpose. First, we split the Lindblad superoperator ac-cording to L = L + L , where L (⋅) = − i Ω [∑ i ˆ σ xi , ⋅] can be regarded as a small perturbation. Focusing now on L , we notice that each ˆ σ zi commutes with all jump op-erators and Hamiltonian terms in it, implying that the z -components of the spins constitute N independent con-served quantities [27]. Hence, the dynamics can be sepa-rated in 2 N independent sectors labeled by the classicalspin configurations ⃗ σ = ( σ , . . . , σ N ) ( σ i ∈ { − , } ), whereˆ σ zi ∣⃗ σ ⟩ = σ i ∣⃗ σ ⟩ ; in other words, states belonging to differ-ent sectors never mix under the action of L . In each sec-tor, the bosonic mode evolves according to a Lindbladian L ( ˆ σ zi → σ i ) which describes a damped quantum har-monic oscillator with a (spin-configuration-dependent)spatial displacement. This admits a single (bosonic) sta-tionary state, denoted by ρ ⃗ σ . We assume that, due tothe random and independent nature of the couplings g i ,no additional symmetries are present which could protectmore complex subspaces. Hence, for any initial state ρ ofthe spin-boson system the corresponding stationary stateunder L is of the form ρ stat = ∑ ⃗ σ p ⃗ σ ρ ⃗ σ ⊗ ∣⃗ σ ⟩ ⟨⃗ σ ∣ , wherethe coefficients p ⃗ σ form a set of classical probabilities.The perturbation L couples sectors corresponding todifferent classical spin configurations ⃗ σ . Its action can beincorporated perturbatively [28, 29] as long as Ω is smallcompared to the typical rate at which coherences betweensectors decay (estimated further below). We proceed byprojecting onto the stationary manifold of L via P ρ ( t ) =∑ ⃗ σ Tr B {⟨⃗ σ ∣ ρ ( t ) ∣⃗ σ ⟩} ρ ⃗ σ ⊗ ∣⃗ σ ⟩ ⟨⃗ σ ∣ , where Tr B denotes thepartial trace over the bosonic mode. This reduces thedynamics to the evolution of the classical probabilities p ⃗ σ ( t ) = Tr B {⟨⃗ σ ∣ ρ ( t ) ∣⃗ σ ⟩} according to a master equation˙ p ⃗ σ = ∑ ⃗ σ ′ W ⃗ σ ′ →⃗ σ p ⃗ σ ′ − W ⃗ σ →⃗ σ ′ p ⃗ σ . Here W ⃗ σ ′ →⃗ σ is the rate forswitching from configuration ⃗ σ ′ to ⃗ σ . Note, that up tosecond order in Ω, the corresponding stochastic processincludes only single spin flips (i.e., W ⃗ σ ′ →⃗ σ ≠ ⃗ σ and ⃗ σ ′ differ by a single spin). The rates read W ⃗ σ →⃗ σ ′ = ω ∫ ∞ dτ e − g i νω ( f ( τ )+ τ ) cos [ ∆ E i τ − g i s ( τ ) ω ( η + ) ] ,f ( τ ) = − η η ( η + ) [ − e − η τ cos ( τ )] − e − η τ η + ( τ ) ,s ( τ ) = η [ e − η τ cos ( τ ) − ] + [ η − ] e − η τ sin ( τ ) η + , (3)where the index i denotes which spin is being flipped andchanges sign between configurations ⃗ σ and ⃗ σ ′ .In Eq. (3) we have introduced the (scaled) differencebetween loss and gain rates η = ( γ − κ )/ ω ≡ γ / ω ( − θ ) ,the ratio θ = κ / γ ∈ [ , ) and the parameter ν = ( + θ ) η /[( η + )( − θ )] . Importantly, the sole depen-dence on the spin configuration is through the quantity∆ E i = g i σ i ∑ l ≠ i g l σ l , which can be interpreted as an en-ergy difference (see further below). Note, that there isa characteristic scale of exponential suppression of theintegrand of W ⃗ σ →⃗ σ ′ . This corresponds to the typicaltimescale involved in the loss of coherence between sec-tors belonging to different classical spin configurations. FIG. 2.
Regimes of effectively thermal spin dynamics. (a) Statistical energy distribution versus time for N =
10 spinsfor large gain-loss difference η , starting from an initial state where all spins point down. The superposed red solid line displaysthe evolution of the average energy, ⟨ E ⟩ ( t ) as a function of time, whereas the shading represents the probability of being ina configuration with energy E at time t . We set Ω = . κ = ω = η =
10 ( θ = /
11) and we select the couplings g i froma uniform distribution in [− g , g ] with g =
3. Inset: evolution of the probabilities p σ for a system of two spins with thesame parameters. We compare the effective model (solid lines) with the numerically-exact diagonalization of the full openquantum problem (dashed lines), highlighting good agreement. (b) Statistical energy distribution versus time for N =
10 spinsand corresponding evolution of the average energy, ⟨ E ⟩ ( t ) (red solid line) for η = . θ = / R i ( ∆ E ) = W i ( ∆ E )/ W i (− ∆ E ) versus energydifference ∣ ∆ E ∣ . For small values of η the rate W i ( ∆ E ) is strongly asymmetric with respect to ∆ E , whereas for large η we have W i ( ∆ E ) ≈ W i (− ∆ E ) . (d) The ratio R i ( ∆ E ) is shown for three different values of ∆ E = . , , g k are drawnfrom a uniform distribution [− g , g ] with g =
6. At small η , the ratio R i ( ∆ E ) approaches one, indicating that configurationsare visited with equal probability as for large η . (e) The three curves in (d) are rescaled according to log R i ( ∆ E )/ ∆ E . Theirasymptotic collapse in the limit η → β eff which governs the dynamicswhen η is sufficiently small. Since the function f ( τ ) is bounded, we can estimate thistimescale as t L ≈ ω /( g ν ) . Our perturbative expansionthus holds as long as Ω ≪ / t L . In the following weperform a detailed investigation of the effective spin dy-namics. It turns out that the loss-gain parameter η iscentral in determining the qualitative dynamical behav-ior: we will identify an effective high-temperature regimein the asymptotic limits η → ∞ and η → + . Furthermore,we find an effective low-temperature (cooling) dynamicswhen η < Large η : infinite temperature dynamics — As remarkedabove, the quantity ∆ E i can be interpreted as the changein the energy function E (⃗ σ ) = − ∑ i ≠ j g i g j σ i σ j . (4)occurring when the i -th spin is flipped, i.e. ∆ E i = E (− σ i ) − E ( σ i ) . In passing, we remark that the energylevels defined by Eq. (4) are (at least) doubly degener-ate, since E (⃗ σ ) = E (−⃗ σ ) . For a large gain-loss difference, η ≫
1, we find that in Eq. (3) f ( τ ) ∼ s ( τ ) ∼ O ( / η ) .Therefore both functions are approximately zero and theparameter ν ≈ + θ − θ η − determines the leading behav-ior of the timescale t L . The validity of the perturbativerequirement thus imposes an upper bound to loss-gaindifference, which must satisfy 1 ≪ η ≪ g ( + θ ) ω Ω ( − θ ) .With the above approximations the rate W ⃗ σ →⃗ σ ′ ac-quires a considerably simpler form: having neglected s ( τ ) , it no longer depends on the sign of ∆ E i , implying that the rates for inverse processes σ → σ ′ and σ ′ → σ areequal. This gives rise to an infinite-temperature dynam-ics which populates all configurations uniformly. Thisbehavior is highlighted in Fig. 2(a): we show that theaverage energy ⟨ E ⟩ ( t ) approaches (up to finite size cor-rections) zero, indicating a equal population of all spinsstates at stationarity.Interestingly, for large η and up to second order inperturbation theory, the rate W ⃗ σ →⃗ σ ′ is formally equiva-lent to the dissipative dynamics of a fictitious transversefield Ising model. The corresponding Hamiltonian isˆ H eff = Ω eff ∑ i ˆ σ xi + ξE ( ˆ σ z ) [with Ω eff = Ω λ , ξ = λ /( ωη ) ]and the spins are subject to strong dephasing at a (site-dependent) rate γ eff , i = g i λ ( + θ ) ωη ( − θ ) [30]. Here λ is an arbi-trary factor that should be chosen consistently with the(perturbative) requirement Ω eff / γ eff , i ≪
1. Therefore, inthis limit, the bosons can be interpreted as forming aninfinite temperature bath causing dephasing of the spindegrees of freedom.
Small η : approximate low-temperature dynamics — For η < κ (gain) and ω , while γ (loss) is varied.Accordingly, the limit η → + has to be interpreted as γ → κ + , so that θ = ( + ηω / κ ) − → FIG. 3.
Non-equilibrium regime.
Kolmogorov criterionin a two-spin subspace; the displayed quantity is the ratio ofthe product ∏ W c of the four rates encountered when per-forming the loop clockwise (blue arrows) divided by the anal-ogous counter-clockwise (red arrows) product ∏ W ac , plottedas a function of η . This ratio is 1 (signalling detailed balanceconditions) only for large and very small values of η . Theintermediate point η ≈ ∏ W c = ∏ W ac can be safelyignored for the following reason: for systems with more thantwo spins there are multiple loops in configuration space andKolmogorov’s criterion is never satisfied simultaneously for allof them (except in the extremal limits η → η → ∞ ) andthe dynamics does not obey detailed balance. The parametersare: ω =
1, Ω = . κ = g = To obtain approximate expressions for the transitionrate we treat it as a function of the energy difference∆ E , which we now regard as a continuous parameter.Furthermore, we note that, for sufficiently small η , theintegrand defining W ⃗ σ →⃗ σ ′ is rapidly suppressed for τ > f ( τ ) ≈ ( − cos τ )/ η .Thus, the integral is dominated by the contribution closeto τ =
0. Hence, one can expand all arguments in powersof τ (see Appendix). Setting for simplicity ω = κ = η →
0, oneobtains τ + f ( τ ) ≈ τ / η − τ / ν ≈ s ( τ ) ≈ − τ + τ / τ ∼ √ η /( g i ) , whereas the cosine termoscillates with a frequency which is approximately Γ i = ( ∆ E + g i ) . We thereby identify (i) a regime of “smallenergy jumps”, where Γ i ≪ g i / η and (ii) a “large energyjumps” one with Γ i ≫ g i / η . In case (i), we obtain (seeAppendix) W i ( ∆ E ) ≈ Ω √ πηg i e − ηg i ( ∆ E + g i ) , (5)where the index i reminds us of which spin is beingflipped [31]. The rate reaches its maximum when ∆ E =− g i ≤ W i (∣ ∆ E ∣) < W i (− ∣ ∆ E ∣) . Thismeans that spin flips which lower the energy are favored,suggesting that the dynamics enacts a form of cooling.Case (ii) can be analyzed using the asymptotic expansionof Fourier integrals [32] (see also Appendix). To leadingorder this yields a power-law decay W i ( ∆ E ) ≈ g i Γ − i ,which shares the same “cooling” properties. A numerical integration suggests that W i (∣ ∆ E ∣) < W i (− ∣ ∆ E ∣) holdsalso in between the asymptotic cases (i) and (ii).To shed further light on the cooling dynamics weanalyze this asymmetry of the rates through the ra-tio R i ( ∆ E ) = W i ( ∆ E )/ W i (− ∆ E ) . This is depicted inFig. 2(c,d) as a function of ∆ E and η , respectively. Inregime (i) we have R i ( ∆ E ) ≈ e − η ∆ E , which implies athermal dynamics with an effective inverse temperature β eff = η . Note that the r.h.s. has no dependence on theindex i , implying the existence of a unique, well-definedtemperature for all spin-flip processes. In Fig. 2(e) wedisplay the ratio log [ R i ( ∆ E )]/ ∆ E for different values of∆ E and show that different curves collapse to a single(negative) inverse temperature − β eff up to the edge ofcase (i). At η = β eff approaches zero, leadingto an infinite-temperature dynamics. This is reasonable,since in this limit the bosonic gain rate approaches theloss rate. This implies the population of arbitrarily highFock states, effectively heating the bosons. The latterthen act as a high-temperature bath on the spins. If, onthe other hand, 1 / β eff remains small or comparable withthe energy gap from the ground states of (4) — which onaverage is of order g , meaning 4 ηg ≥
1) — an effectivelow-temperature dynamics is realized.In case (ii), the ratio R i ( ∆ E ) ≈ ( ∆ E − g i ) /( ∆ E + g i ) tends to increase towards 1 as ∆ E grows. Typically, theavailable ∆ E i s populate both range (i) and (ii), implyingthe presence of type (ii) processes which do not follow thesame low-temperature rules obeyed by the “small-jump”ones. Provided the number of spins N is not too large,these non-thermal processes constitute, however, a smallperturbation for the following two reasons: first, the dis-tribution of energy jumps is peaked around 0, imply-ing that, if the parameters are adequately chosen, mostjumps lie in regime (i). Second, since the rates are de-creasing functions of ∣ ∆ E + g i ∣ , type (ii) processes occurat smaller rates than the type (i), thermal ones. A nu-merically exact analysis of the classical master equationfor N =
10, displayed in Fig. 2(b), indeed shows that theeffect of the non-thermal processes is sufficiently weakto avoid having a significant population of high-energystates in the long-time limit. The statistical energy dis-tribution tends instead to become concentrated on low-energy configurations, highlighting a clear bias of the dy-namics towards cooling, as compared e.g., to the η ≫ Breakdown of detailed balance —
Outside the thermalregimes the dynamics is not an equilibrium one, i.e. itdoes not obey detailed balance. This can be proved viaKolmogorov’s criterion [33] which we analyze for the loopformed in the configuration space of a two-spin system(see Fig. 3): (↑↑) → (↑↓) → (↓↓) → (↓↑) → (↑↑) . Tothis end we investigate the ratio between the product ofthe rates for the clockwise (blue arrows) cycle and thecorresponding product for the counter-clockwise (red ar-rows) one. This ratio goes to 1 when η → ∞ and alsowhen η →
0, signalling the emergence of the infinite-temperature dynamics. For different values η the ratio istypically different from one, which indicates the persis-tence of probability currents in the stationary state andthe absence of detailed balance. Conclusions —
We have studied a disordered dissipa-tive spin-boson system in the limit of strong coupling andweak driving, which can for example be implemented ontrapped ion quantum simulators. Many aspects of theemerging physics can be understood in terms of a disor-dered fully-connected Ising model whose state evolves ac-cording to a rate equation. In general the dynamics vio-lates detailed balance, and the system is thus out of equi-librium. However, we could identify parameter regimes inwhich the evolution is effectively thermal. Among themis one where predominantly low-energy configurations arepopulated,which mimics the action of a low-temperaturedynamics. In the future it would be interesting to seewhether this effective cooling mechanism permits to ac-cess low-energy states or even ground states of complexspin networks. This might open an elegant way for en-coding and solving computationally hard problems [17] orassociative memories [18] through Ising energy functionsand an appropriate (thermal) dynamics on quantum sim-ulators.
Acknowledgments—
The research leading to theseresults has received funding from the European Re-search Council under the European Unions SeventhFramework Programme (FP/2007-2013)/ERC GrantAgreement No. 335266 (ESCQUMA), the EPSRCGrant No.EP/R04340X/1 via the QuantERA project“ERyQSenS and from the University of Nottinghamthrough a Nottingham Research Fellowship (M.M.). P.R.acknowledges funding by the European Union throughthe H2020 - MCIF No. 766442. I.L. gratefully acknowl-edges funding through the Royal Society Wolfson Re-search Merit Award. [1] A. Blais, R.-S. Huang, A. Wallraff, S. M. Girvin, andR. J. Schoelkopf, Phys. Rev. A , 062320 (2004).[2] H. Zheng, D. J. Gauthier, and H. U. Baranger, Phys.Rev. Lett. , 090502 (2013).[3] J. I. Cirac and P. Zoller, Phys. Rev. Lett. , 4091 (1995).[4] D. Leibfried, R. Blatt, C. Monroe, and D. Wineland,Rev. Mod. Phys. , 281 (2003).[5] A. F. Kockum, A. Miranowicz, S. De Liberato,S. Savasta, and F. Nori, Nature Reviews Physics , 19(2019).[6] E. Farhi, J. Goldstone, S. Gutmann, J. Lapan, A. Lund-gren, and D. Preda, Science , 472 (2001).[7] V. S. Denchev, S. Boixo, S. V. Isakov, N. Ding, R. Bab-bush, V. Smelyanskiy, J. Martinis, and H. Neven, Phys.Rev. X , 031015 (2016).[8] P. Strack and S. Sachdev, Phys. Rev. Lett. , 277202(2011). [9] S. Gopalakrishnan, B. Lev, and P. Goldbart, Phys. Rev.Lett. , 277201 (2011).[10] P. Rotondo, E. Tesio, and S. Caracciolo, Phys. Rev. B , 014415 (2015).[11] P. Rotondo, M. Cosentino Lagomarsino, and G. Viola,Phys. Rev. Lett. , 143601 (2015).[12] E. G. D. Torre, S. Diehl, M. D. Lukin, S. Sachdev, andP. Strack, Phys. Rev. A , 023831 (2013).[13] E. Fiorelli, P. Rotondo, M. Marcuzzi, J. P. Garrahan,and I. Lesanovsky, Phys. Rev. A , 032126 (2019).[14] P. Rotondo, M. Marcuzzi, J. P. Garrahan, I. Lesanovsky,and M. M¨uller, Journal of Physics A: Mathematical andTheoretical , 115301 (2018).[15] V. D. Vaidya, Y. Guo, R. M. Kroeze, K. E. Ballantine,A. J. Koll´ar, J. Keeling, and B. L. Lev, Phys. Rev. X ,011002 (2018).[16] Y. Guo, R. M. Kroeze, V. V. Vaidya, J. Keeling, andB. L. Lev, arXiv preprint arXiv:1810.11086 (2018).[17] F. Barahona, Journal of Physics A: Mathematical andGeneral , 3241 (1982).[18] J. J. Hopfield, PNAS , 2554 (1982).[19] B. M. Garraway, Philosophical Transactions of the RoyalSociety A: Mathematical, Physical and Engineering Sci-ences , 1137 (2011).[20] P. Kirton, M. M. Roses, J. Keeling, and E. G.Dalla Torre, Advanced Quantum Technologies , 1800043(2019).[21] K. Hepp and E. H. Lieb, Annals of Physics , 360(1973).[22] K. Hepp and E. H. Lieb, Phys. Rev. A , 2517 (1973).[23] D. Porras and J. I. Cirac, Phys. Rev. Lett. , 207901(2004).[24] I. Aedo and L. Lamata, Physical Review A , 042317(2018).[25] A. Mezzacapo, U. Las Heras, J. Pedernales, L. DiCarlo,E. Solano, and L. Lamata, Scientific Reports , 7482(2014).[26] G.-D. Lin, S.-L. Zhu, R. Islam, K. Kim, M.-S. Chang,S. Korenblit, C. Monroe, and L.-M. Duan, Europhys.Lett. , 60004 (2009).[27] V. V. Albert and L. Jiang, Phys. Rev. A , 022118(2014).[28] P. Degenfeld-Schonburg and M. J. Hartmann, Phys. Rev.B , 245108 (2014).[29] M. Marcuzzi, J. Schick, B. Olmos, and I. Lesanovsky,Journal of Physics A: Mathematical and Theoretical ,482001 (2014).[30] B. Everest, M. Marcuzzi, J. P. Garrahan, andI. Lesanovsky, Phys. Rev. E , 052108 (2016).[31] Note, that in this notation if W ⃗ σ →⃗ σ ′ ≡ W i ( ∆ E ) then W ⃗ σ ′ →⃗ σ ≡ W i (− ∆ E ) ).[32] H. Dai and D. Naylor, Proceedings of the Royal Societyof London. Series A: Mathematical and Physical Sciences , 109 (1992).[33] R. K. P. Zia and B. Schmittmann, Journal of Statisti-cal Mechanics: Theory and Experiment , P07012(2007). Appendices for ”Dynamics of strongly coupled disordered dissipative spin-bosonsystems”
DERIVATION OF THE RATES
In this section we provide details on the derivation of Eq.(3) in the main text. We firstly consider the evolution ofthe state ρ ( t ) as ˙ ρ = L ( ρ ) + L ( ρ ) , with L = L − L and L (⋅) = − i Ω [∑ i ˆ σ xi , ⋅] . Secondly, we assume the stationarystate of L of the form ρ stat = ∑ ⃗ σ p ⃗ σ ρ ⃗ σ ∣⃗ σ ⟩ ⟨⃗ σ ∣ , where ∣⃗ σ ⟩ = { σ , ..., σ N } , with σ i = ± σ zi ∣⃗ σ ⟩ = σ i ∣⃗ σ ⟩ , p ⃗ σ are a set ofclassical probabilities and ρ σ is the corresponding bosonic state, that we assume to be a gaussian state. Considering L perturbatively with respect to L , and projecting the dynamics onto the stationary manifold of L , we exploit theNakajima-Zwanzig formalism to write the evolution of the spin as P ˙ ρ spin ( t ) = Tr B ∫ +∞ dt ′ P L e L T ′ L ρ stat ( t ) == Ω ∑ {⃗ σ } p ⃗ σ ( t ) ∑ i ∫ +∞ dt ′ ∑ j =± Tr B [ e V j ⃗ σ,i t ′ ( ρ ⃗ σ )] ×× ( ˆ σ xi ∣⃗ σ ⟩ ⟨⃗ σ ∣ ˆ σ xi − ∣⃗ σ ⟩ ⟨⃗ σ ∣) , (S1)where P ˙ ρ spin ( t ) = Tr B [ ˙ ρ stat ( t )] = ∑ ⃗ σ ˙ p ⃗ σ ( t ) ∣⃗ σ ⟩ ⟨⃗ σ ∣ , Tr B is the partial trace over the boson, and we have defined thespin-configuration dependent superoperators V ±⃗ σ,i (⋅) = − iω [ ˆ a † ˆ a, ⋅] + D γ (⋅) + D κ (⋅) − ig i M i [( ˆ a † + ˆ a ) , ⋅] ± ig i σ i {( ˆ a † + ˆ a ) , ⋅} , with M i = ∑ l ≠ i g l σ l , and D γ , D κ the dissipative terms representing cooling and heating, respectively. By projectingEq.(S1) on a state ∣⃗ σ ′ ⟩ , the dynamics reduces to the evolution of the classical probabilities ruled by a master equationwhose general form is the following ˙ p ⃗ σ = ∑ ⃗ σ ′ ( W ⃗ σ ′ →⃗ σ p ⃗ σ ′ − W ⃗ σ →⃗ σ ′ p ⃗ σ ) , (S2)where W ⃗ σ →⃗ σ ′ = Ω ∫ +∞ dt ∑ j =± Tr B [ e V j ⃗ σ,i τ ( ρ ⃗ σ )] , is the transition rate for the switching ⃗ σ → ⃗ σ ′ and it allows onlysingle spin-flip processes.We can now go ahead in evaluating explicitly the expression for the rates. Exploiting the superoperator’s properties,we can write e V ±⃗ σ,i ( ρ ⃗ σ ) = ( e V ± , ∗⃗ σ,i ) ρ ⃗ σ , with V ± , ∗⃗ σ,i the adjoint superoperator of V ±⃗ σ,i and I the identity operator. It isworth noticing that the identity operator can be expressed in terms of a generalised displacement operator of fieldcoherent states as I = ˆ D ( ) , where ˆ D ( τ ) = e α ( τ ) ˆ a † − β ( τ ) ∗ ˆ a e − γ ( τ ) with α ( ) = β ( ) = γ ( ) =
0. We then verify thatthe displacement operator ˆ D ( ) is mapped into the generalised one ˆ D ( τ ) by applying the adjoint superoperator V ±⃗ σ,i :indeed, by considering the differential equation ddτ [ ˆ D ±⃗ σ,i ( τ )] = V ±∗⃗ σ,i [ ˆ D ±⃗ σ,i ( τ )] we obtain the solutions for the functions α + σ i ( τ ) = α −− σ i ( τ ) , β + σ i ( τ ) = β −− σ i ( τ ) , γ + σ i ( τ ) = γ −− σ i ( τ ) . For initial conditions α ± σ i ( ) = β ± σ i ( ) = γ ± σ i ( ) =
0, we get α + σ i ( τ ) = [ β + σ i ( t )] ∗ = i g i σ i ω ( η − i ) [ − e ( i − η ) τ ] ,γ + σ i ( τ ) = g i νω η [ f ( τ ) + τ ] + ig i σ i M i ω ( η + ) [ s ( τ ) + τ ] ,f ( τ ) = − e ητ η − η [ − e − η τ cos ( τ )] − e − η τ sin ( τ ) η + ,s ( τ ) = η [ e − η τ cos ( τ ) − ] + ( η − ) e − η τ sin ( τ ) η + , (S3)where we have defined the dimensionless time τ = tω , and η = ( γ − κ )/ ω , θ = κ / γ ∈ [ , ) , and ν = ( + θ ) η ( η + )( − θ ) = ( κω + η ) η + . (S4)The previous steps allow us to write the partial trace over the boson as Tr B [ e V ±⃗ σ,i ( ρ ⃗ σ )] = e − γ ± σi ( τ ) Tr B [ ˆ D ±⃗ σ,i ( τ ) ρ ⃗ σ ] .We recall that the bosonic state ρ ⃗ σ has been assumed to be a gaussian state. In this case, we recognise the quantityTr B [ ˆ D ±⃗ σ,i ( τ ) ρ ⃗ σ ] as the characteristic function χ ±⃗ σ,i ( τ ) of the state ρ ⃗ σ . The expression of the characteristic functionfor a generic gaussian state ρ G reads χ ρ G [ α ( τ )] = e − ⃗ α T ( τ ) Σ ⃗ α ( τ )+ α ( τ )⟨ ˆ a † ⟩ G − α ∗ ( τ )⟨ ˆ a ⟩ G , (S5)where ⃗ α T = ( α ( τ ) , α ∗ ( τ )) , ⟨⋅⟩ G is the expectation value performed over the state ρ G , and Σ represent the covariancematrix which reads 2 ( −(⟨ ˆ a ⟩ G − ⟨ ˆ a ⟩ G ) (⟨ ˆ a † ˆ a ⟩ G + ⟨ ˆ a ˆ a † ⟩ G ) − ⟨ ˆ a † ⟩ G ⟨ ˆ a ⟩ G (⟨ ˆ a † ˆ a ⟩ c + ⟨ ˆ a ˆ a † ⟩ G ) − ⟨ ˆ a † ⟩ G ⟨ ˆ a ⟩ G −(⟨ ˆ a † ⟩ G − ⟨ ˆ a † ⟩ G ) ) . (S6)By applying the definition (S5), with expectation values of the operators obtained considering the lindblad operator L , we get χ ±⃗ σ,i ( τ ) = exp {− + θ ( − θ ) ∣ α ± σ i ( τ )∣ + ig i M ω ( α ± σ i ( τ ) η − i + α ±∗ σ i ( τ ) η + i )} , (S7)where M = ∑ l g l σ l . Thus, the expression of the rate reads W ⃗ σ →⃗ σ ′ = Ω ω ∫ +∞ dτ ∑ j =± e − γ jσi ( τ ) χ j ⃗ σ,i ( τ ) = ω ∫ ∞ dτ e − g i νω [ f ( τ )+ τ ] cos [ ∆ E i τ − g i s ( τ ) ω ( η + ) ] ,f ( τ ) = − η + η ( η + ) [ − e − η τ cos ( τ )] − e − η τ η + ( τ ) , (S8)where ∆ E i = g i σ i ∑ l ≠ i g l σ l retains the dependence on the spin configuration. APPROXIMATE EXPRESSION OF THE RATE FOR η SMALL
For sufficiently small η ≲
1, the exponent appearing in W ⃗ σ →⃗ σ ′ is dominated by f ( τ ) ≈ η ( − cos ( τ )) , (S9)implying that the integrand is quickly suppressed as τ grows. We therefore perform an expansion around τ = τ + f ( τ ) = ( η + ) ( η τ − τ − − η η τ ) + O ( τ ) (S10)and s ( τ ) = − τ + η + τ − η η + τ + O ( τ ) . (S11)From the first term in Eq. (S10) we see that the integrand is strongly suppressed on scales τ ≳ √ η . Noticing that inthe Taylor expansion of τ + f ( τ ) odd coefficients are finite, whereas even ones are O ( / η ) for η → z = τ √ η we see that higher orders are perturbations of order O ( η n + / z n + , η n − z n ) and can be neglected. Similar considerations can be applied to s ( τ ) , which can be therefore also approximated withits leading order − τ . In the following, for simplicity we set ω =
1, remembering that our “energy” ∆ E i is actuallymeasured by construction in units of ω . Additionally, we introduce the shorthandΓ i = η + ( ∆ E i + g i ) , (S12)so that, by keeping only the lowest orders of the expansion in τ , we can approximate our rate as W ⃗ σ →⃗ σ ′ ≈ ∫ ∞ dτ e − g i ( κ + η ) η τ cos ( Γ i τ ) , (S13)which can be integrated to give the closed expression W ⃗ σ →⃗ σ ′ ≈ Ω √ πη g i ( κ + η ) e − η Γ2 i g i ( κ + η ) . (S14)It is worth remarking that the exponent can be rewritten as − η Γ i g i ( κ + η ) = η ( κ + η )( η + ) [− ( ∆ E i + g i ) g i ] == η ( κ + η )( η + ) [ E (⃗ σ ) − ∑ j g j ] = β eff E (⃗ σ ) − const. , (S15)highlighting the “thermal” structure of the rates. Note that, in order to obtain the approximation (S13), we haveassumed that we can resum the Taylor expansion of the original cosine to the function cos ( Γ i τ ) , whose series onlycoincides with the former up to O ( τ ) . This is only valid as long as the cosine does not oscillate significantly beforethe other Gaussian term suppresses the integrand; in other words, Eq. (S14) should be valid up to values of Γ i ofthe order of ∼ /√ η . Since we wish to understand the behavior of the rates as functions of the energy difference∆ E i without restrictions imposed by the other parameters (like η ), we need to account for energies which exceed thisrange. To do this, we extract the asymptotic behavior of the rate for Γ i → ∞ . We start by rewriting the integrand in W as I ( Γ i , τ ) ≡ Re { e − g i νω ( f ( τ )+ τ ) e i Γ i τ − i g iη + ( s ( τ )+ τ ) } = (S16) = Re { A ( τ ) e i Γ i τ } . (S17)We now use the result that, if the function A admits a small τ expansion A ( τ ) = ∞ ∑ n = a n τ n , (S18)then asymptotically in the limit Γ i → ∞ one finds ∫ ∞ dτ A ( τ ) e i Γ i τ = ∞ ∑ n = i n n ! a n Γ − n − i . (S19)The leading term in this expansion corresponds to the lowest n for which one finds a non-vanishing real part. Inparticular, we note that Re [ a n i n + ] equals (− ) l + Re [ a l + ] if n = l + (− ) l + Im [ a l ] if n = l is even.For our function we find a = , (S20a) a = , (S20b) a = − g i κ + ηη , (S20c) a = g i [ κ + η − i ] . (S20d)The leading behavior in the large Γ i limit is therefore determined by Re [ a ] , implying W ⃗ σ →⃗ σ ′ = ∫ ∞ dτ I ( Γ i , τ ) ≈ g i ( κ + η ) [ η + ( ∆ E i + g i ) ] . (S21)To provide a very crude estimate of where the change from the two regimes characterized by Eqs. (S14) (“small Γ i ”)and (S21) (“large Γ i ”) occurs, we evaluate the point where the two asymptotic expressions cross (for η sufficientlysmall): setting 4 g i ( κ + η ) Γ − i = √ πη g i ( κ + η ) e − η Γ2 i g i ( κ + η ) (S22)we find Γ i ≈ g i ( κ + η ) η [ log A + (
12 log A )] , (S23)where A = g i π ( κ + η ) η −3