Dynamics of strongly coupled disordered dissipative spin-boson systems

# Dynamics 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 Centre for the Mathematics and Theoretical Physics of Quantum Non-equilibrium Systems, University of Nottingham, Nottingham NG7 2RD, UK
July 2, 2019
###### Abstract

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 Blais et al. (2004) or waveguide quantum electrodynamics Zheng et al. (2013) as well as trapped ions Cirac and Zoller (1995); Leibfried et al. (2003) 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 Kockum et al. (2019).

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 Farhi et al. (2001); Denchev et al. (2016).

Disordered spin-boson systems have only recently moved into the focus of theoretical investigations. References Strack and Sachdev (2011); Gopalakrishnan et al. (2011) explore the emergence of glassiness when many electromagnetic modes interact with an ensemble of qubits. In Refs. Rotondo et al. (2015a, b), 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 Torre et al. (2013); Fiorelli et al. (2019); Rotondo et al. (2018). This topic acquires further relevance in the light of recent experimental progress in multimodal cavity QED, which realize tunable range Vaidya et al. (2018) and sign-changing Guo et al. (2018) photon-mediated atomic interactions.

In this work we investigate a disordered and dissipative 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 dynamical 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 Barahona (1982) or associative memories Hopfield (1982) — together with a cooling protocol.

Model — We consider an ensemble of two-level systems interacting with a single bosonic mode described by the following Dicke Hamiltonian Garraway (2011); Kirton et al. (2019); Hepp and Lieb (1973a, b):

 ^H=ω^a†^a+N∑i=1gi^σzi(^a†+^a)+ΩN∑i=1^σxi. (1)

Here, 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 are assumed to be independent and randomly distributed with zero mean and variance .

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 equation

 ˙ρ=Lρ=−i[^H,ρ]+∑n=l,g^Lnρ^L†n−12{^L†n^Ln,ρ} (2)

with the jump operators , where () is the loss (gain) rate and .

A physical realization of this model can for instance be achieved on trapped-ion quantum simulators Porras and Cirac (2004); Aedo and Lamata (2018): Following the scheme represented in Fig. 1, such system would consist of ions coupling to the centre-of-mass phonon mode. As it has been shown for the quantum Rabi model Mezzacapo et al. (2014) and eventually generalized to the Dicke model Aedo and Lamata (2018), the application of multiple laser fields on the ions yields both the spin dependent coupling and the weak driving term, 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. Lin et al. (2009). 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 (). In the following, we sketch the perturbative technique we employ for this purpose. First, we split the Lindblad superoperator according to , where can be regarded as a small perturbation. Focusing now on , we notice that each commutes with all jump operators and Hamiltonian terms in it, implying that the -components of the spins constitute independent conserved quantities Albert and Jiang (2014). Hence, the dynamics can be separated in independent sectors labeled by the classical spin configurations (), where ; in other words, states belonging to different sectors never mix under the action of . In each sector, the bosonic mode evolves according to a Lindbladian 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 , no additional symmetries are present which could protect more complex subspaces. Hence, for any initial state of the spin-boson system the corresponding stationary state under is of the form , where the coefficients form a set of classical probabilities.

The perturbation couples sectors corresponding to different classical spin configurations . Its action can be incorporated perturbatively Degenfeld-Schonburg and Hartmann (2014); Marcuzzi et al. (2014) 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 via , where denotes the partial trace over the bosonic mode. This reduces the dynamics to the evolution of the classical probabilities according to a master equation . Here 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., only if and differ by a single spin). The rates read

 W→σ→→σ′=2Ω2ω∫∞0dτe−2g2iνω2(f(τ)+τ)\medsizecos[16ΔEiτ−g2is(τ)ω2(η2+4)], f(τ)=8−2η2η(η2+4)[1−e−η2τcos(τ)]−8e−η2τη2+4sin(τ), s(τ)=4η[e−η2τcos(τ)−1]+[η2−4]e−η2τsin(τ)η2+4, (3)

where the index 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 , the ratio and the parameter . Importantly, the sole dependence on the spin configuration is through the quantity , 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 . This corresponds to the typical timescale involved in the loss of coherence between sectors belonging to different classical spin configurations. Since the function is bounded, we can estimate this timescale as . Our perturbative expansion thus holds as long as . 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 . Furthermore, we find an effective low-temperature (cooling) dynamics when .

Large : infinite temperature dynamics — As remarked above, the quantity can be interpreted as the change in the energy function

 E(→σ)=−14∑i≠jgigjσiσj. (4)

occurring when the -th spin is flipped, i.e. . In passing, we remark that the energy levels defined by Eq. (4) are (at least) doubly degenerate, since . For a large gain-loss difference, , we find that in Eq. (3) . Therefore both functions are approximately zero and the parameter determines the leading behavior of the timescale . The validity of the perturbative requirement thus imposes an upper bound to loss-gain difference, which must satisfy .

With the above approximations the rate acquires a considerably simpler form: having neglected , it no longer depends on the sign of , 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 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 is formally equivalent to the dissipative dynamics of a fictitious transverse field Ising model. The corresponding Hamiltonian is [with , ] and the spins are subject to strong dephasing at a (site-dependent) rate Everest et al. (2016). Here is an arbitrary factor that should be chosen consistently with the (perturbative) requirement . 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 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 has to be interpreted as , so that 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 , which we now regard as a continuous parameter. Furthermore, we note that, for sufficiently small , the integrand defining is rapidly suppressed for due to a fast initial increase of . Thus, the integral is dominated by the contribution close to . Hence, one can expand all arguments in powers of (see Appendix). Setting for simplicity and keeping for brevity only the leading orders in , one obtains , and . This implies that the suppression of the integrand occurs on a timescale , whereas the cosine term oscillates with a frequency which is approximately . We thereby identify (i) a regime of “small energy jumps”, where and (ii) a “large energy jumps” one with . In case (i), we obtain (see Appendix)

 Wi(ΔE) ≈ Ω22√πηg2ie−ηg2i(ΔE+g2i)2, (5)

where the index reminds us of which spin is being flipped 111Note, that in this notation if then ).. The rate reaches its maximum when and, in general, . 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 Dai and Naylor (1992) (see also Appendix). To leading order this yields a power-law decay , which shares the same “cooling” properties. A numerical integration suggests that 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 . This is depicted in Fig. 2(c,d) as a function of and , respectively. In regime (i) we have , which implies a thermal dynamics with an effective inverse temperature . Note that the r.h.s. has no dependence on the index , implying the existence of a unique, well-defined temperature for all spin-flip processes. In Fig. 2(e) we display the ratio for different values of and show that different curves collapse to a single (negative) inverse temperature up to the edge of case (i). At we have 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, remains small or comparable with the energy gap from the ground states of (4) — which on average is of order , meaning ) — an effective low-temperature dynamics is realized.

In case (ii), the ratio tends to increase towards as grows. Typically, the available 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 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 , implying that, if the parameters are adequately chosen, most jumps lie in regime (i). Second, since the rates are decreasing functions of , type (ii) processes occur at smaller rates than the type (i), thermal ones. A numerically exact analysis of the classical master equation for , 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 low-energy configurations, highlighting a clear bias of the dynamics towards cooling, as compared e.g., to the 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 Zia and Schmittmann (2007) 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 when and also when , signalling the emergence of the infinite-temperature 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 Barahona (1982) or associative memories Hopfield (1982) through Ising energy functions and an appropriate (thermal) dynamics on quantum simulators.

Acknowledgments— The research leading to these results has received funding from the European Research Council under the European Unions Seventh Framework Programme (FP/2007-2013)/ERC Grant Agreement No. 335266 (ESCQUMA), the EPSRC Grant No.EP/R04340X/1 via the QuantERA project “ERyQSenS and from the University of Nottingham through a Nottingham Research Fellowship (M.M.). P.R. acknowledges funding by the European Union through the H2020 - MCIF No. . I.L. gratefully acknowledges funding through the Royal Society Wolfson Research Merit Award.

## References

Appendices for ”Dynamics of strongly coupled disordered dissipative spin-boson systems”

## Appendix A 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 as , with and . Secondly, we assume the stationary state of of the form , where , with and , are a set of classical probabilities and is the corresponding bosonic state, that we assume to be a gaussian state. Considering perturbatively with respect to , and projecting the dynamics onto the stationary manifold of , we exploit the Nakajima-Zwanzig formalism to write the evolution of the spin as

 P˙ρspin(t)=TrB∫+∞0dt′PL1eL0T′L1ρstat(t)==Ω2∑{→σ}p→σ(t)∑i∫+∞0dt′∑j=±TrB[eVj→σ,it′(ρ→σ)]××(^σxi|→σ⟩⟨→σ|^σxi−|→σ⟩⟨→σ|), (S1)

where , is the partial trace over the boson, and we have defined the spin-configuration dependent superoperators

 V±→σ,i(⋅)=−iω[^a†^a,⋅]+Dγ(⋅)+Dκ(⋅)−igiMi[(^a†+^a),⋅]±igiσi{(^a†+^a),⋅},

with , and , 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

 ˙p→σ=∑→σ′(W→σ′→→σp→σ′−W→σ→→σ′p→σ), (S2)

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 , with the adjoint superoperator of and 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 , where with . We then verify that the displacement operator is mapped into the generalised one by applying the adjoint superoperator : indeed, by considering the differential equation we obtain the solutions for the functions , , . For initial conditions , we get

 α+σi(τ)=[β+σi(t)]∗=i4giσiω(η−2i)[1−e(i−η2)τ],γ+σi(τ)=2g2iνω2η[f1(τ)+τ]+igiσiMiω2(η2+4)[s(τ)+τ],f1(τ)=1−eητη−4η[1−e−η2τcos(τ)]−8e−η2τsin(τ)η2+4,s(τ)=4η[e−η2τcos(τ)−1]+(η2−4)e−η2τsin(τ)η2+4, (S3)

where we have defined the dimensionless time , and , , and

 ν=4(1+θ)η(η2+4)(1−θ)=4(2κω+η)η2+4. (S4)

The previous steps allow us to write the partial trace over the boson as . We recall that the bosonic state has been assumed to be a gaussian state. In this case, we recognise the quantity as the characteristic function of the state . The expression of the characteristic function for a generic gaussian state reads

 (S5)

where , is the expectation value performed over the state , and represent the covariance matrix which reads

 (S6)

By applying the definition (S5), with expectation values of the operators obtained considering the lindblad operator , we get

 χ±→σ,i(τ)=exp{−1+θ2(1−θ)|α±σi(τ)|2+2igiMω(α±σi(τ)η−2i+α±∗σi(τ)η+2i)}, (S7)

where . Thus, the expression of the rate reads

 W→σ→→σ′=Ω2ω∫+∞0dτ∑j=±e−γjσi(τ)χj→σ,i(τ)=2Ω2ω∫∞0dτe−2g2iνω2[f(τ)+τ]\medsizecos[16ΔEiτ−g2is(τ)ω2(η2+4)],f(τ)=−2η2+8η(η2+4)[1−e−η2τcos(τ)]−8e−η2τη2+4sin(τ), (S8)

where retains the dependence on the spin configuration.

## Appendix B Approximate expression of the rate for η small

For sufficiently small , the exponent appearing in is dominated by

 f(τ)≈2η(1−cos(τ)), (S9)

implying that the integrand is quickly suppressed as grows. We therefore perform an expansion around which yields

 τ+f(τ)=(η2+4)(14ητ2−124τ3−4−η2192ητ4)+O(τ5) (S10)

and

 s(τ)=−τ+η2+424τ3−ηη2+496τ4+O(τ5). (S11)

From the first term in Eq. (S10) we see that the integrand is strongly suppressed on scales . Noticing that in the Taylor expansion of odd coefficients are finite, whereas even ones are for and introducing the rescaled integration variable we see that higher orders are perturbations of order and can be neglected. Similar considerations can be applied to , which can be therefore also approximated with its leading order . In the following, for simplicity we set , remembering that our “energy” is actually measured by construction in units of . Additionally, we introduce the shorthand

 Γi=16η2+4(ΔEi+g2i), (S12)

so that, by keeping only the lowest orders of the expansion in , we can approximate our rate as

 W→σ→→σ′≈2Ω2∫∞0dτe−2g2i(2κ+η)ητ2cos(Γiτ), (S13)

which can be integrated to give the closed expression

 W→σ→→σ′≈Ω2√πη2g2i(2κ+η)e−ηΓ2i8g2i(2κ+η). (S14)

It is worth remarking that the exponent can be rewritten as

 −ηΓ2i8g2i(2κ+η)=128η(2κ+η)(η2+4)2[−(ΔEi+g2i)24g2i]= =128η(2κ+η)(η2+4)2[E(→σ)−∑jg2j4]=βeffE(→σ)−const., (S15)

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 , whose series only coincides with the former up to . 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 of the order of . Since we wish to understand the behavior of the rates as functions of the energy difference 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 . We start by rewriting the integrand in as

 I(Γi,τ) ≡Re⎧⎪⎨⎪⎩e−2g2iνω2(f(τ)+τ)eiΓiτ−16ig2iη2+4(s(τ)+τ)⎫⎪⎬⎪⎭= (S16) =Re{A(τ)eiΓiτ}. (S17)

We now use the result that, if the function admits a small expansion

 A(τ)=∞∑n=0anτn, (S18)

then asymptotically in the limit one finds

 ∫∞0dτA(τ)eiΓiτ=∞∑n=0inn!anΓ−n−1i. (S19)

The leading term in this expansion corresponds to the lowest for which one finds a non-vanishing real part. In particular, we note that equals if is odd, and if is even. For our function we find

 a0 =1, (S20a) a1 =0, (S20b) a2 =−2g2i2κ+ηη, (S20c) a3 =g2i3[2κ+η−2i]. (S20d)

The leading behavior in the large limit is therefore determined by , implying

 W→σ→→σ′=2Ω2∫∞0dτI(Γi,τ)≈4Ω2g2i(2κ+η)[η2+416(ΔEi+g2i)]4. (S21)

To provide a very crude estimate of where the change from the two regimes characterized by Eqs. (S14) (“small ”) and (S21) (“large ”) occurs, we evaluate the point where the two asymptotic expressions cross (for sufficiently small): setting

 4g2i(2κ+η)Γ−4i=√πη2g2i(2κ+η)e−ηΓ2i8g2i(2κ+η) (S22)

we find

 Γ2i≈4g2i(2κ+η)η[logA+4log(12logA)], (S23)

where .

You are adding the first comment!
How to quickly get a good reply:
• Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
• Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
• Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
The feedback must be of minimum 40 characters and the title a minimum of 5 characters