Dynamics of strongly coupled disordered dissipative spin-boson systems

Spin-boson Hamiltonians are an effective description for numerous quantum many-body systems such as atoms coupled to cavity modes, quantum electrodynamics in circuits and trapped ion systems. While reaching the limit of strong coupling is possible in current experiments, the understanding of the physics in this parameter regime remains a challenge, especially when disorder and dissipation are taken into account. Here we investigate a regime where the many-body spin dynamics can be related to a Ising energy function defined in terms of the spin-boson couplings. While in the coherent weak coupling regime it is known that an effective description in terms of spin Hamiltonian is possible, we show that a similar viewpoint can be adopted in the presence of dissipation and 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 that under 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 for engineering complex energy landscapes together with an appropriate cooling dynamics.

Introduction-Prominent platforms for quantum simulation, such as cavity, circuit [1] or waveguide quantum electrodynamics [2] as well as trapped ions [3,4] can be modeled by ensembles of two-level systems interacting via bosonic degrees of freedom (electromagnetic modes or phonons). While the weak coupling regime is relatively well understood and can be treated by a perturbative integration of the bosonic degrees of freedom, the strong coupling limit is far more challenging [5].
An additional layer of complexity is added by the presence of disorder, i.e. when individual spins couple to the bosonic "environment" at different strengths. Such a setting is relevant for at least two reasons. First, some degree of quenched disorder may always be present in realistic systems and, second, one may engineer non-uniform couplings for practical applications: systems with tunable quasi-random couplings often form the basis for a physical implementation of complex optimization problems, which may for instance be solved via quantum annealing protocols [6,7].
Disordered spin-boson systems have only recently moved into the focus of theoretical investigations. References [8,9] explore the emergence of glassiness when many electromagnetic modes interact with an ensemble of qubits. In Refs. [10,11], instead, spin-glass techniques are employed to show that the same system effectively realizes an associative memory. Most of these techniques, however, cannot be straightforwardly generalized to study open quantum dynamics in the strong coupling regime, and only a few studies deal with disordered open quantum systems [12][13][14]. This topic acquires further relevance in the light of recent experimental progress in multimodal cavity QED, which realize tunable range [15] and sign-changing [16] photon-mediated atomic interactions.
In this work we investigate a disordered and dissipa- tive system in which weakly driven spins are strongly coupled to a bosonic mode (see Fig. 1a). We employ a perturbative approach which relies on the weakness of the driving rather than of the spin-boson coupling. We find that the effective spin dynamics is governed by a rate equation that depends on a fully-connected Ising energy function as sketched in Fig. 1b. Depending on the rates of bosonic loss and gain we identify several distinct dynam-ical regimes: two of them are high-temperature ones, in which the stationary state of the system is fully mixed. A further one mimics an effective low-temperature dynamics, which permits cooling of the spin system. Outside these the dynamics is generally non-thermal and detailed balance is broken. This link between an open, strongly coupled spin-boson system and the physics of disordered Ising spin systems opens up the possibility of engineering complex classical energy landscapes -with importance in the context of optimization problems [17] or associative memories [18] -together with a cooling protocol. Model -We consider an ensemble of N two-level systems interacting with a single bosonic mode described by the following Dicke Hamiltonian [19][20][21][22]: Here,σ x,y,z i are the Pauli operators andâ andâ † the bosonic annihilation and creation operators. The parameters ω and Ω denote the fundamental frequency of the bosons and the coherent coupling strength between the two spin states, respectively. The spin-boson couplings g i are assumed to be independent and randomly distributed with zero mean and variance g 2 .
We include dissipation on the boson in the form of Markovian gain and loss processes. The density matrix ρ of this open quantum system therefore evolves under a Lindblad equatioṅ with the jump operatorsL l = √ γâ,L g = √ κâ † where γ (κ) is the loss (gain) rate and γ > κ ≥ 0.
A physical realization of this model can for instance be achieved on trapped-ion quantum simulators [23,24]: Following the scheme represented in Fig. 1, such system would consist of N ions coupling to the centre-of-mass phonon mode. As it has been shown for the quantum Rabi model [25] and eventually generalized to the Dicke model [24], the application of multiple laser fields on the ions yields both the spin dependent coupling g kσ z k (â+â † ) and the weak driving term, Ωσ x k entering Eq. (1). Finally, as illustrated in Fig. 1a, the gain and loss dynamics can be achieved by applying lasers on the ions on the edge of the chain, which is discussed in Ref. [26]. Since this ion is coupled to the same phonon mode as the other ions this effectively implements jump operators of the form introduced in Eq. (2).
Spin dynamics at strong coupling -We explore the dynamics (2) in the strong coupling regime, i.e. when the driving acting on the spins is much weaker than the spin-boson interaction (Ω ≪ g). In the following, we sketch the perturbative technique we employ for this purpose. First, we split the Lindblad superoperator according to L = L 0 + L 1 , where L 1 (⋅) = −iΩ[∑ iσ x i , ⋅] can be regarded as a small perturbation. Focusing now on L 0 , we notice that eachσ z i commutes with all jump operators and Hamiltonian terms in it, implying that the z-components of the spins constitute N independent conserved quantities [27]. Hence, the dynamics can be separated in 2 N independent sectors labeled by the classical in other words, states belonging to different sectors never mix under the action of L 0 . In each sector, the bosonic mode evolves according to a Lindbladian L 0 (σ z i → σ i ) which describes a damped quantum harmonic oscillator with a (spin-configuration-dependent) spatial displacement. This admits a single (bosonic) stationary state, denoted by ρ ⃗ σ . We assume that, due to the random and independent nature of the couplings g i , no additional symmetries are present which could protect more complex subspaces. Hence, for any initial state ρ 0 of the spin-boson system the corresponding stationary state σ form a set of classical probabilities. The perturbation L 1 couples sectors corresponding to different classical spin configurations ⃗ σ. Its action can be incorporated perturbatively [28,29] as long as Ω is small compared to the typical rate at which coherences between sectors decay (estimated further below). We proceed by projecting onto the stationary manifold of where Tr B denotes the partial trace over the bosonic mode. This reduces the dynamics to the evolution of the classical probabilities Here W ⃗ σ ′ →⃗ σ is the rate for switching from configuration ⃗ σ ′ to ⃗ σ. Note, that up to second order in Ω, the corresponding stochastic process includes only single spin flips (i.e., W ⃗ σ ′ →⃗ σ ≠ 0 only if ⃗ σ and ⃗ σ ′ differ by a single spin). The rates read where the index i denotes which spin is being flipped and changes sign between configurations ⃗ σ and ⃗ σ ′ . In Eq. (3) we have introduced the (scaled) difference between loss and gain rates η = (γ − κ) ω ≡ γ ω(1 − θ), the ratio θ = κ γ ∈ [0, 1) and the parameter ν = 4(1 + θ)η [(η 2 + 4)(1 − θ)]. Importantly, the sole dependence on the spin configuration is through the quantity ∆E i = g i σ i ∑ l≠i g l σ l , which can be interpreted as an energy difference (see further below). Note, that there is a characteristic scale of exponential suppression of the integrand of W ⃗ σ→⃗ σ ′ . This corresponds to the typical timescale involved in the loss of coherence between sectors belonging to different classical spin configurations. Since the function f (τ ) is bounded, we can estimate this timescale as t L ≈ ω (2g 2 ν). Our perturbative expansion thus holds as long as Ω ≪ 1 t L . In the following we perform a detailed investigation of the effective spin dynamics. It turns out that the loss-gain parameter η is central in determining the qualitative dynamical behavior: we will identify an effective high-temperature regime in the asymptotic limits η → ∞ and η → 0 + . Furthermore, we find an effective low-temperature (cooling) dynamics when η < 1.
Large η: infinite temperature dynamics -As remarked above, the quantity ∆E i can be interpreted as the change in the energy function occurring when the i-th spin is flipped, i.e.
In passing, we remark that the energy levels defined by Eq. (4) are (at least) doubly degenerate, since E(⃗ σ) = E(−⃗ σ). For a large gain-loss difference, η ≫ 1, we find that in Eq. ( Therefore both functions are approximately zero and the parameter ν ≈ 4 1+θ 1−θ η −1 determines the leading behavior of the timescale t L . The validity of the perturbative requirement thus imposes an upper bound to loss-gain difference, which must satisfy 1 ≪ η ≪ 4g 2 (1+θ) ωΩ(1−θ) . With the above approximations the rate W ⃗ σ→⃗ σ ′ acquires 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 σ ′ → σ are equal. This gives rise to an infinite-temperature dynamics which populates all configurations uniformly. This behavior is highlighted in Fig. 2(a): we show that the average energy ⟨E⟩ (t) approaches (up to finite size corrections) zero, indicating a equal population of all spins states at stationarity.
Interestingly, for large η and up to second order in perturbation theory, the rate W ⃗ σ→⃗ σ ′ is formally equivalent to the dissipative dynamics of a fictitious transverse field Ising model. The corresponding Hamiltonian iŝ H eff = Ω eff ∑ iσ x i +ξE(σ z ) [with Ω eff = Ωλ, ξ = 8λ 2 (ωη 2 )] and the spins are subject to strong dephasing at a (site- ωη(1−θ) [30]. Here λ is an arbitrary factor that should be chosen consistently with the (perturbative) requirement Ω eff γ eff,i ≪ 1. Therefore, in this limit, the bosons can be interpreted as forming an infinite temperature bath causing dephasing of the spin degrees of freedom.
Small η: approximate low-temperature dynamics -For η < 1 there exists a regime in which the rate equation dynamics mimics a thermal process with finite temperature. To be precise, this limit is achieved by fixing the parameters κ (gain) and ω, while γ (loss) is varied. Accordingly, the limit η → 0 + has to be interpreted as γ → κ + , so that θ = (1 + ηω κ) −1 → 1 remains finite. Provided the parameters are chosen carefully, this leads to the cooling of the spins with respect to the Ising energy function (4) [see Fig. 2b]. To obtain approximate expressions for the transition rate 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 η, the integrand defining W ⃗ σ→⃗ σ ′ is rapidly suppressed for τ > 0 due to a fast initial increase of f (τ ) ≈ 2(1 − cos τ ) η. Thus, the integral is dominated by the contribution close to τ = 0. Hence, one can expand all arguments in powers of τ (see Appendix). Setting for simplicity ω = κ = 1 and keeping for brevity only the leading orders in η → 0, one obtains τ + f (τ ) ≈ τ 2 η − τ 3 6, ν ≈ 2 and s(τ ) ≈ −τ + τ 3 6. This implies that the suppression of the integrand occurs on a timescale τ ∼ η (4g 2 i ), whereas the cosine term oscillates with a frequency which is approximately Γ i = 4(∆E + g 2 i ). We thereby identify (i) a regime of "small energy jumps", where Γ 2 i ≪ 4g 2 i η and (ii) a "large energy jumps" one with Γ 2 i ≫ 4g 2 i η. In case (i), we obtain (see Appendix) where the index i reminds us of which spin is being flipped [31]. The rate reaches its maximum when ∆E = −g 2 i ≤ 0 and, in general, W i ( ∆E ) < W i (− ∆E ). This means 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 expansion of Fourier integrals [32] (see also Appendix). To leading order this yields a power-law decay W i (∆E) ≈ 8Ω 2 g 2 i Γ −4 i , which shares the same "cooling" properties. A numerical integration suggests that W i ( ∆E ) < W i (− ∆E ) holds also in between the asymptotic cases (i) and (ii).
To shed further light on the cooling dynamics we analyze this asymmetry of the rates through the ratio R i (∆E) = W i (∆E) W i (−∆E). This is depicted in Fig. 2(c,d) as a function of ∆E and η, respectively. In regime (i) we have R i (∆E) ≈ e −4η∆E , which implies a thermal dynamics with an effective inverse temperature β eff = 4η. Note that the r.h.s. has no dependence on the index i, implying the existence of a unique, well-defined temperature for all spin-flip processes. In Fig. 2(e) we display 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 of case (i). At η = 0 we have β eff approaches zero, leading to an infinite-temperature dynamics. This is reasonable, since in this limit the bosonic gain rate approaches the loss rate. This implies the population of arbitrarily high Fock states, effectively heating the bosons. The latter then act as a high-temperature bath on the spins. If, on the other hand, 1 β eff remains small or comparable with the energy gap from the ground states of (4) -which on average is of order g 2 , meaning 4ηg 2 ≥ 1) -an effective low-temperature dynamics is realized.
In case (ii), the ratio R i (∆E) ≈ (∆E −g 2 i ) 4 (∆E +g 2 i ) 4 tends to increase towards 1 as ∆E grows. Typically, the available ∆E i s populate both range (i) and (ii), implying the presence of type (ii) processes which do not follow the same 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 small perturbation for the following two reasons: first, the distribution of energy jumps is peaked around 0, implying that, if the parameters are adequately chosen, most jumps lie in regime (i). Second, since the rates are decreasing functions of ∆E + g 2 i , type (ii) processes occur at smaller rates than the type (i), thermal ones. A numerically exact analysis of the classical master equation for N = 10, displayed in Fig. 2(b), indeed shows that the effect of the non-thermal processes is sufficiently weak to avoid having a significant population of high-energy states in the long-time limit. The statistical energy distribution tends instead to become concentrated on lowenergy configurations, highlighting a clear bias of the dynamics towards cooling, as compared e.g., to the η ≫ 1 case in panel (a).
Breakdown of detailed balance -Outside the thermal regimes the dynamics is not an equilibrium one, i.e. it does not obey detailed balance. This can be proved via Kolmogorov's criterion [33] which we analyze for the loop formed in the configuration space of a two-spin system (see Fig. 3): (↑↑) → (↑↓) → (↓↓) → (↓↑) → (↑↑). To this end we investigate the ratio between the product of the rates for the clockwise (blue arrows) cycle and the corresponding product for the counter-clockwise (red arrows) one. This ratio goes to 1 when η → ∞ and also when η → 0, signalling the emergence of the infinitetemperature dynamics. For different values η the ratio is typically different from one, which indicates the persistence of probability currents in the stationary state and the absence of detailed balance.
Conclusions -We have studied a disordered dissipative spin-boson system in the limit of strong coupling and weak driving, which can for example be implemented on trapped ion quantum simulators. Many aspects of the emerging physics can be understood in terms of a disordered fully-connected Ising model whose state evolves according to a rate equation. In general the dynamics violates detailed balance, and the system is thus out of equilibrium. However, we could identify parameter regimes in which the evolution is effectively thermal. Among them is one where predominantly low-energy configurations are populated,which mimics the action of a low-temperature dynamics. In the future it would be interesting to see whether this effective cooling mechanism permits to access low-energy states or even ground states of complex spin networks. This might open an elegant way for encoding and solving computationally hard problems [17] or associative memories [18] through Ising energy functions and an appropriate (thermal) dynamics on quantum simulators.
Appendices for "Dynamics of strongly coupled disordered dissipative spin-boson systems"

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 of the state ρ(t) asρ = L 0 (ρ) + L 1 (ρ), with L 0 = L − L 1 and L 1 (⋅) = −iΩ[∑ iσ x i , ⋅]. Secondly, we assume the stationary state of L 0 of the form σ are a set of classical probabilities and ρ σ is the corresponding bosonic state, that we assume to be a gaussian state. Considering L 1 perturbatively with respect to L 0 , and projecting the dynamics onto the stationary manifold of L 0 , we exploit the Nakajima-Zwanzig formalism to write the evolution of the spin as where Pρ spin (t) = Tr B [ρ stat (t)] = ∑ ⃗ σṗ⃗ σ (t) ⃗ σ⟩ ⟨⃗ σ , Tr B is the partial trace over the boson, and we have defined the spin-configuration dependent superoperators with M i = ∑ l≠i g l σ l , and D γ , D κ the dissipative terms representing cooling and heating, respectively. By projecting Eq.(S1) on a state ⃗ σ ′ ⟩, the dynamics reduces to the evolution of the classical probabilities ruled by a master equation whose general form is the followingṗ where , is the transition rate for the switching ⃗ σ → ⃗ σ ′ and it allows only single 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 1)ρ ⃗ σ , with V ±, * ⃗ σ,i the adjoint superoperator of V ± ⃗ σ,i and I the identity operator. It is worth noticing that the identity operator can be expressed in terms of a generalised displacement operator of field coherent states as I =D(0), whereD(τ ) = e α(τ )â † −β(τ ) * â e −γ(τ ) with α(0) = β(0) = γ(0) = 0. We then verify that the displacement operatorD(0) is mapped into the generalised oneD(τ ) by applying the adjoint superoperator V ± ⃗ σ,i : indeed, by considering the differential equation where we have defined the dimensionless time τ = tω, and η = (γ − κ) ω, θ = κ γ ∈ [0, 1), and The previous steps allow us to write the partial trace over the boson as Tr . We recall that the bosonic state ρ ⃗ σ has been assumed to be a gaussian state. In this case, we recognise the quantity which can be integrated to give the closed expression It is worth remarking that the exponent can be rewritten as highlighting the "thermal" structure of the rates. Note that, in order to obtain the approximation (S13), we have assumed that we can resum the Taylor expansion of the original cosine to the function cos(Γ i τ ), whose series only coincides with the former up to O(τ 2 ). This is only valid as long as the cosine does not oscillate significantly before the other Gaussian term suppresses the integrand; in other words, Eq. (S14) should be valid up to values of Γ i of the order of ∼ 1 √ η. 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 this range. To do this, we extract the asymptotic behavior of the rate for Γ i → ∞. We start by rewriting the integrand in W as We now use the result that, if the function A admits a small τ expansion A(τ ) = ∞ n=0 a n τ n , (S18) then asymptotically in the limit Γ i → ∞ one finds i n n! a n Γ −n−1 i . (S19) The leading term in this expansion corresponds to the lowest n for which one finds a non-vanishing real part. In particular, we note that Re[a n i n+1 ] equals (−1) l+1 Re[a 2l+1 ] if n = 2l + 1 is odd, and (−1) l+1 Im[a 2l ] if n = 2l is even. For our function we find The leading behavior in the large Γ i limit is therefore determined by Re[a 3 ], implying 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 η sufficiently small): setting we find where A = 32g 2 i π(2κ + η)η −3 .