Dynamics of strongly coupled disordered dissipative spin-boson systems

TL;DR

This paper explores the complex dynamics of strongly coupled, disordered, dissipative spin-boson systems, revealing thermalization, out-of-equilibrium behavior, and potential for cooling, with implications for quantum energy landscape engineering.

Contribution

It demonstrates that a spin Hamiltonian description remains valid in strong coupling and dissipation regimes, enabling new insights into their many-body dynamics.

Findings

Identification of approximately thermal regimes in dissipative strong coupling

Discovery of out-of-equilibrium dynamics breaking detailed balance

Potential for cooling spin degrees of freedom under specific conditions

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.

Figures

Equations

$$\hat{H}=\omega {\hat{a}}^{†}\hat{a}+\sum _{i=1}^N{g}_i{\hat{\sigma }}_i^z({\hat{a}}^{†}+\hat{a})+\Omega \sum _{i=1}^N{\hat{\sigma }}_i^x.$$

$$\dot{\rho }=\mathcal{L}\rho =-i[\hat{H},\rho ]+\sum _{n=l,g}{\hat{L}}_n^{}\rho {\hat{L}}_n^{†}-\frac{1}{2}\{{\hat{L}}_n^{†}{\hat{L}}_n^{},\rho \}$$

$$W_{\vec{\sigma }\to {\vec{\sigma }}^{\prime }}=\frac{2{\Omega }^2}{\omega }{\int }_0^{\infty }d\tau e^{-\frac{2g_i^2\nu }{{\omega }^2}\left(f(\tau )+\tau \right)}\text{medsize}\cos \left[16\frac{\Delta E_i\tau -g_i^{2}s(\tau )}{{\omega }^2({\eta }^2+4)}\right],$$

$$f(\tau )=\frac{8-2{\eta }^2}{\eta \left({\eta }^2+4\right)}\left[1-e^{-\frac{\eta }{2}\tau }\cos (\tau )\right]-\frac{8e^{-\frac{\eta }{2}\tau }}{{\eta }^2+4}\sin (\tau ),$$

$$s(\tau )=\frac{4\eta \left[e^{-\frac{\eta }{2}\tau }\cos (\tau )-1\right]+\left[{\eta }^2-4\right]e^{-\frac{\eta }{2}\tau }\sin (\tau )}{{\eta }^2+4},$$

$$E(\vec{\sigma })=-\frac{1}{4}\sum _{i\ne j}{g}_i{g}_j{\sigma }_i{\sigma }_j.$$

$$W_i(\Delta E)$$

$$\begin{array}{rl}P{\dot{\rho }}^{\mathrm{spin}}(t)=& {\mathrm{Tr}}_B{\int }_0^{+\infty }d{t}^{\prime }P{\mathcal{L}}_1{e}^{{\mathcal{L}}_0{T}^{\prime }}{\mathcal{L}}_1{\rho }_{\mathrm{stat}}(t)=\\ =& {\Omega }^2\sum _{\{\vec{\sigma }\}}{p}_{\vec{\sigma }}(t)\sum _i{\int }_0^{+\infty }d{t}^{\prime }\sum _{j=\pm }{\mathrm{Tr}}_B\left[e^{{\mathcal{V}}_{\vec{\sigma },i}^j{t}^{\prime }}({\rho }_{\vec{\sigma }})\right]\times \\ & \times \left({\hat{\sigma }}_i^x|\vec{\sigma }⟩⟨\vec{\sigma }|{\hat{\sigma }}_i^x-|\vec{\sigma }⟩⟨\vec{\sigma }|\right),\end{array}$$

$${\mathcal{V}}_{\vec{\sigma },i}^{\pm }(\cdot )=-i\omega [{\hat{a}}^{†}\hat{a},\cdot ]+{\mathcal{D}}_{\gamma }\left(\cdot \right)+{\mathcal{D}}_{\kappa }\left(\cdot \right)-i{g}_i{\mathcal{M}}_i\left[({\hat{a}}^{†}+\hat{a}),\cdot \right]\pm i{g}_i{\sigma }_i\left\{({\hat{a}}^{†}+\hat{a}),\cdot \right\},$$

$${\dot{p}}_{\vec{\sigma }}=\sum _{{\vec{\sigma }}^{\prime }}\left(W_{{\vec{\sigma }}^{\prime }\to \vec{\sigma }}{p}_{{\vec{\sigma }}^{\prime }}-W_{\vec{\sigma }\to {\vec{\sigma }}^{\prime }}{p}_{\vec{\sigma }}\right),$$

$$\begin{array}{rl}& {\alpha }_{{\sigma }_i}^{+}(\tau )=[{\beta }_{{\sigma }_i}^{+}(t){]}^{\ast }=\frac{i4g_i{\sigma }_i}{\omega (\eta -2i)}\left[1-e^{(i-\frac{\eta }{2})\tau }\right],\\ & {\gamma }_{{\sigma }_i}^{+}(\tau )=\frac{2g_i^2\nu }{{\omega }^2\eta }\left[f_1(\tau )+\tau \right]+\frac{i{g}_i{\sigma }_i{\mathcal{M}}_i}{{\omega }^2({\eta }^2+4)}[s(\tau )+\tau ],\\ & f_1(\tau )=\frac{1-e^{\eta \tau }}{\eta }-\frac{4\eta [1-e^{-\frac{\eta }{2}\tau }\cos (\tau )]-8e^{-\frac{\eta }{2}\tau }\sin (\tau )}{{\eta }^2+4},\\ & s(\tau )=\frac{4\eta \left[e^{-\frac{\eta }{2}\tau }\cos (\tau )-1\right]+\left({\eta }^2-4\right)e^{-\frac{\eta }{2}\tau }\sin (\tau )}{{\eta }^2+4},\end{array}$$

$$\nu =\frac{4(1+\theta )\eta }{({\eta }^2+4)(1-\theta )}=\frac{4\left(2\frac{\kappa }{\omega }+\eta \right)}{{\eta }^2+4}.$$

$${\chi }_{{\rho }_G}[\alpha (\tau )]=e^{-\frac{1}{4}{\vec{\alpha }}^T(\tau )\Sigma \vec{\alpha }(\tau )+\alpha (\tau ){⟨{\hat{a}}^{†}⟩}_G-{\alpha }^{\ast }(\tau ){⟨\hat{a}⟩}_G},$$

$$2\left(\begin{array}{ccc}-({⟨{\hat{a}}^2⟩}_G-{⟨\hat{a}⟩}_G^2)& & \frac{1}{2}({⟨{\hat{a}}^{†}\hat{a}⟩}_G+{⟨\hat{a}{\hat{a}}^{†}⟩}_G)-{⟨{\hat{a}}^{†}⟩}_G{⟨\hat{a}⟩}_G\\ \frac{1}{2}({⟨{\hat{a}}^{†}\hat{a}⟩}_c+{⟨\hat{a}{\hat{a}}^{†}⟩}_G)-{⟨{\hat{a}}^{†}⟩}_G{⟨\hat{a}⟩}_G& & -({⟨{\hat{a}}^{†2}⟩}_G-{⟨{\hat{a}}^{†}⟩}_G^2)\end{array}\right).$$

$${\chi }_{\vec{\sigma },i}^{\pm }(\tau )=\exp \left\{-\frac{1+\theta }{2(1-\theta )}|{\alpha }_{{\sigma }_i}^{\pm }(\tau ){|}^2+\frac{2i{g}_i\mathcal{M}}{\omega }\left(\frac{{\alpha }_{{\sigma }_i}^{\pm }(\tau )}{\eta -2i}+\frac{{\alpha }_{{\sigma }_i}^{\pm \ast }(\tau )}{\eta +2i}\right)\right\},$$

$$\begin{array}{rl}{W}_{\vec{\sigma }\to {\vec{\sigma }}^{\prime }}=& \frac{{\Omega }^2}{\omega }{\int }_0^{+\infty }d\tau \sum _{j=\pm }{e}^{-{\gamma }_{{\sigma }_i}^j(\tau )}{\chi }_{\vec{\sigma },i}^j(\tau )=\frac{2{\Omega }^2}{\omega }{\int }_0^{\infty }d\tau e^{-\frac{2g_i^2\nu }{{\omega }^2}\left[f(\tau )+\tau \right]}\text{medsize}\cos \left[16\frac{\Delta E_i\tau -g_i^{2}s(\tau )}{{\omega }^2({\eta }^2+4)}\right],\\ & f(\tau )=\frac{-2{\eta }^2+8}{\eta \left({\eta }^2+4\right)}\left[1-e^{-\frac{\eta }{2}\tau }\cos (\tau )\right]-\frac{8e^{-\frac{\eta }{2}\tau }}{{\eta }^2+4}\sin (\tau ),\end{array}$$

$$f(\tau )\approx \frac{2}{\eta }(1-\cos (\tau )),$$

$$\tau +f(\tau )=({\eta }^2+4)\left(\frac{1}{4\eta }{\tau }^2-\frac{1}{24}{\tau }^3-\frac{4-{\eta }^2}{192\eta }{\tau }^4\right)+O({\tau }^5)$$

$$s(\tau )=-\tau +\frac{{\eta }^2+4}{24}{\tau }^3-\eta \frac{{\eta }^2+4}{96}{\tau }^4+O({\tau }^5).$$

$${\Gamma }_i=\frac{16}{{\eta }^2+4}(\Delta E_i+g_i^2),$$

$$W_{\vec{\sigma }\to {\vec{\sigma }}^{\prime }}\approx 2{\Omega }^2{\int }_0^{\infty }d\tau e^{-2\frac{g_i^2\left(2\kappa +\eta \right)}{\eta }{\tau }^2}\cos \left({\Gamma }_i\tau \right),$$

$$W_{\vec{\sigma }\to {\vec{\sigma }}^{\prime }}\approx {\Omega }^2\sqrt{\frac{\pi \eta }{2g_i^2\left(2\kappa +\eta \right)}}{e}^{-\frac{\eta {\Gamma }_i^2}{8g_i^2\left(2\kappa +\eta \right)}}.$$

$$-\frac{\eta {\Gamma }_i^2}{8g_i^2\left(2\kappa +\eta \right)}=\frac{128\eta }{(2\kappa +\eta )({\eta }^2+4{)}^2}\left[-\frac{(\Delta E_i+g_i^2{)}^2}{4g_i^2}\right]=$$

$$=\frac{128\eta }{(2\kappa +\eta )({\eta }^2+4{)}^2}\left[E(\vec{\sigma })-\frac{{\sum }_j{g}_j^2}{4}\right]={\beta }_{\mathrm{eff}}E(\vec{\sigma })-const.,$$

$$I({\Gamma }_i,\tau )$$

$$=\mathrm{Re}\left\{A(\tau )e^{i{\Gamma }_i\tau }\right\}.$$

$$A(\tau )=\sum _{n=0}^{\infty }{a}_n{\tau }^n,$$

$${\int }_0^{\infty }d\tau A(\tau )e^{i{\Gamma }_i\tau }=\sum _{n=0}^{\infty }{i}^{n}n!a_n{\Gamma }_i^{-n-1}.$$

$$a_0$$

$$a_1$$

$$a_2$$

$$a_3$$

$$W_{\vec{\sigma }\to {\vec{\sigma }}^{\prime }}=2{\Omega }^2{\int }_0^{\infty }d\tau I({\Gamma }_i,\tau )\approx 4{\Omega }^2{g}_i^2(2\kappa +\eta ){\left[\frac{{\eta }^2+4}{16(\Delta E_i+g_i^2)}\right]}^4.$$

$$4g_i^2(2\kappa +\eta ){\Gamma }_i^{-4}=\sqrt{\frac{\pi \eta }{2g_i^2\left(2\kappa +\eta \right)}}{e}^{-\frac{\eta {\Gamma }_i^2}{8g_i^2\left(2\kappa +\eta \right)}}$$

Full text

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 $N$ 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):

[TABLE]

Here, $\hat{\sigma}_{i}^{x,y,z}$ are the Pauli operators and $\hat{a}$ and $\hat{a}^{\dagger}$ the bosonic annihilation and creation operators. The parameters $\omega$ and $\Omega$ 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 $\rho$ of this open quantum system therefore evolves under a Lindblad equation

[TABLE]

with the jump operators $\hat{L}_{l}=\sqrt{\gamma}\hat{a}$, $\hat{L}_{g}=\sqrt{\kappa}\hat{a}^{\dagger}$ where $\gamma$ ($\kappa$) is the loss (gain) rate and $\gamma>\kappa\geq 0$.

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 $N$ 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 $g_{k}\hat{\sigma}_{k}^{z}(\hat{a}+\hat{a}^{\dagger})$ and the weak driving term, $\Omega\hat{\sigma}_{k}^{x}$ 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 ($\Omega\ll g$). In the following, we sketch the perturbative technique we employ for this purpose. First, we split the Lindblad superoperator according to $\mathcal{L}=\mathcal{L}_{0}+\mathcal{L}_{1}$, where $\mathcal{L}_{1}(\cdot)=-i\Omega[\sum_{i}\hat{\sigma}_{i}^{x},\cdot]$ can be regarded as a small perturbation. Focusing now on $\mathcal{L}_{0}$, we notice that each $\hat{\sigma}_{i}^{z}$ commutes with all jump operators and Hamiltonian terms in it, implying that the $z$-components of the spins constitute $N$ independent conserved quantities Albert and Jiang (2014). Hence, the dynamics can be separated in $2^{N}$ independent sectors labeled by the classical spin configurations $\vec{\sigma}=(\sigma_{1},\ldots,\sigma_{N})$ ($\sigma_{i}\in\set{-1,1}$), where $\hat{\sigma}_{i}^{z}\ket{\vec{\sigma}}=\sigma_{i}\ket{\vec{\sigma}}$; in other words, states belonging to different sectors never mix under the action of $\mathcal{L}_{0}$. In each sector, the bosonic mode evolves according to a Lindbladian $\mathcal{L}_{0}(\hat{\sigma}^{z}_{i}\to\sigma_{i})$ which describes a damped quantum harmonic oscillator with a (spin-configuration-dependent) spatial displacement. This admits a single (bosonic) stationary state, denoted by $\rho_{\vec{\sigma}}$. 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 $\rho_{0}$ of the spin-boson system the corresponding stationary state under $\mathcal{L}_{0}$ is of the form $\rho_{\mathrm{stat}}=\sum_{\vec{\sigma}}p_{\vec{\sigma}}\rho_{\vec{\sigma}}\otimes\ket{\vec{\sigma}}\bra{\vec{\sigma}}$, where the coefficients $p_{\vec{\sigma}}$ form a set of classical probabilities.

The perturbation $\mathcal{L}_{1}$ couples sectors corresponding to different classical spin configurations $\vec{\sigma}$. Its action can be incorporated perturbatively Degenfeld-Schonburg and Hartmann (2014); Marcuzzi et al. (2014) as long as $\Omega$ 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 $\mathcal{L}_{0}$ via $P\rho(t)=\sum_{\vec{\sigma}}\mathrm{Tr}_{B}\left\{\bra{\vec{\sigma}}\rho(t)\ket{\vec{\sigma}}\right\}\rho_{\vec{\sigma}}\otimes\ket{\vec{\sigma}}\bra{\vec{\sigma}}$, where $\mathrm{Tr}_{B}$ denotes the partial trace over the bosonic mode. This reduces the dynamics to the evolution of the classical probabilities $p_{\vec{\sigma}}(t)=\mathrm{Tr}_{B}\left\{\bra{\vec{\sigma}}\rho(t)\ket{\vec{\sigma}}\right\}$ according to a master equation $\dot{p}_{\vec{\sigma}}=\sum_{\vec{\sigma}^{\prime}}W_{\vec{\sigma}^{\prime}\rightarrow\vec{\sigma}}p_{\vec{\sigma}^{\prime}}-W_{\vec{\sigma}\rightarrow\vec{\sigma}^{\prime}}p_{\vec{\sigma}}$. Here $W_{\vec{\sigma}^{\prime}\rightarrow\vec{\sigma}}$ is the rate for switching from configuration $\vec{\sigma}^{\prime}$ to $\vec{\sigma}$. Note, that up to second order in $\Omega$, the corresponding stochastic process includes only single spin flips (i.e., $W_{\vec{\sigma}^{\prime}\rightarrow\vec{\sigma}}\neq 0$ only if $\vec{\sigma}$ and $\vec{\sigma}^{\prime}$ differ by a single spin). The rates read

[TABLE]

where the index $i$ denotes which spin is being flipped and changes sign between configurations $\vec{\sigma}$ and $\vec{\sigma}^{\prime}$.

In Eq. (3) we have introduced the (scaled) difference between loss and gain rates $\eta=(\gamma-\kappa)/\omega\equiv\gamma/\omega(1-\theta)$, the ratio $\theta=\kappa/\gamma\in[0,1)$ and the parameter $\nu=4(1+\theta)\eta/[(\eta^{2}+4)(1-\theta)]$. Importantly, the sole dependence on the spin configuration is through the quantity $\Delta E_{i}~{}=g_{i}\sigma_{i}\sum_{l\neq i}g_{l}\sigma_{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_{\vec{\sigma}\rightarrow\vec{\sigma}^{\prime}}$. This corresponds to the typical timescale involved in the loss of coherence between sectors belonging to different classical spin configurations. Since the function $f(\tau)$ is bounded, we can estimate this timescale as $t_{\mathrm{L}}\approx\omega/(2g^{2}\nu)$. Our perturbative expansion thus holds as long as $\Omega\ll 1/t_{\mathrm{L}}$. In the following we perform a detailed investigation of the effective spin dynamics. It turns out that the loss-gain parameter $\eta$ is central in determining the qualitative dynamical behavior: we will identify an effective high-temperature regime in the asymptotic limits $\eta\to\infty$ and $\eta\to 0^{+}$. Furthermore, we find an effective low-temperature (cooling) dynamics when $\eta<1$.

Large $\eta$: infinite temperature dynamics — As remarked above, the quantity $\Delta E_{i}$ can be interpreted as the change in the energy function

[TABLE]

occurring when the $i$-th spin is flipped, i.e. $\Delta E_{i}=E(-\sigma_{i})-E(\sigma_{i})$. In passing, we remark that the energy levels defined by Eq. (4) are (at least) doubly degenerate, since $E(\vec{\sigma})=E(-\vec{\sigma})$. For a large gain-loss difference, $\eta\gg 1$, we find that in Eq. (3) $f(\tau)\sim s(\tau)\sim O(1/\eta)$. Therefore both functions are approximately zero and the parameter $\nu\approx 4\tfrac{1+\theta}{1-\theta}\eta^{-1}$ determines the leading behavior of the timescale $t_{\mathrm{L}}$. The validity of the perturbative requirement thus imposes an upper bound to loss-gain difference, which must satisfy $1\ll\eta\ll\tfrac{4g^{2}(1+\theta)}{\omega\Omega(1-\theta)}$.

With the above approximations the rate $W_{\vec{\sigma}\rightarrow\vec{\sigma}^{\prime}}$ acquires a considerably simpler form: having neglected $s(\tau)$, it no longer depends on the sign of $\Delta E_{i}$, implying that the rates for inverse processes $\sigma\to\sigma^{\prime}$ and $\sigma^{\prime}\to\sigma$ 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 $\braket{E}(t)$ approaches (up to finite size corrections) zero, indicating a equal population of all spins states at stationarity.

Interestingly, for large $\eta$ and up to second order in perturbation theory, the rate $W_{\vec{\sigma}\rightarrow\vec{\sigma}^{\prime}}$ is formally equivalent to the dissipative dynamics of a fictitious transverse field Ising model. The corresponding Hamiltonian is $\hat{H}_{\rm eff}=\Omega_{\rm eff}\sum_{i}\hat{\sigma}^{x}_{i}+\xi E(\hat{\sigma}^{z})$ [with $\Omega_{\rm eff}=\Omega\lambda$, $\xi=8\lambda^{2}/(\omega\eta^{2})$] and the spins are subject to strong dephasing at a (site-dependent) rate $\gamma_{\rm eff,i}=\tfrac{8g_{i}^{2}\lambda^{2}(1+\theta)}{\omega\eta(1-\theta)}$ Everest et al. (2016). Here $\lambda$ is an arbitrary factor that should be chosen consistently with the (perturbative) requirement $\Omega_{\rm eff}/\gamma_{\rm eff,i}\ll 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 $\eta$: approximate low-temperature dynamics — For $\eta<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 $\kappa$ (gain) and $\omega$, while $\gamma$ (loss) is varied. Accordingly, the limit $\eta\to 0^{+}$ has to be interpreted as $\gamma\to\kappa^{+}$, so that $\theta=(1+\eta\omega/\kappa)^{-1}\to 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 $\Delta E$, which we now regard as a continuous parameter. Furthermore, we note that, for sufficiently small $\eta$, the integrand defining $W_{\vec{\sigma}\rightarrow\vec{\sigma}^{\prime}}$ is rapidly suppressed for $\tau>0$ due to a fast initial increase of $f(\tau)\approx 2(1-\cos\tau)/\eta$. Thus, the integral is dominated by the contribution close to $\tau=0$. Hence, one can expand all arguments in powers of $\tau$ (see Appendix). Setting for simplicity $\omega=\kappa=1$ and keeping for brevity only the leading orders in $\eta\to 0$, one obtains $\tau+f(\tau)\approx\tau^{2}/\eta-\tau^{3}/6$, $\nu\approx 2$ and $s(\tau)\approx-\tau+\tau^{3}/6$. This implies that the suppression of the integrand occurs on a timescale $\tau\sim\sqrt{\eta/(4g_{i}^{2})}$, whereas the cosine term oscillates with a frequency which is approximately $\Gamma_{i}=4(\Delta E+g_{i}^{2})$. We thereby identify (i) a regime of “small energy jumps”, where $\Gamma_{i}^{2}\ll 4g_{i}^{2}/\eta$ and (ii) a “large energy jumps” one with $\Gamma_{i}^{2}\gg 4g_{i}^{2}/\eta$. In case (i), we obtain (see Appendix)

[TABLE]

where the index $i$ reminds us of which spin is being flipped 111Note, that in this notation if $W_{\vec{\sigma}\rightarrow\vec{\sigma}^{\prime}}\equiv W_{i}(\Delta E)$ then $W_{\vec{\sigma}^{\prime}\rightarrow\vec{\sigma}}\equiv W_{i}(-\Delta E)$).. The rate reaches its maximum when $\Delta E=-g_{i}^{2}\leq 0$ and, in general, $W_{i}(\left|\Delta E\right|)<W_{i}(-\left|\Delta E\right|)$. 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 $W_{i}(\Delta E)\approx 8\Omega^{2}g_{i}^{2}\Gamma_{i}^{-4}$, which shares the same “cooling” properties. A numerical integration suggests that $W_{i}(\left|\Delta E\right|)<W_{i}(-\left|\Delta E\right|)$ 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}(\Delta E)=W_{i}(\Delta E)/W_{i}(-\Delta E)$. This is depicted in Fig. 2(c,d) as a function of $\Delta E$ and $\eta$, respectively. In regime (i) we have $R_{i}(\Delta E)\approx e^{-4\eta\Delta E}$, which implies a thermal dynamics with an effective inverse temperature $\beta_{\rm eff}=4\eta$. 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}(\Delta E)]/\Delta E$ for different values of $\Delta E$ and show that different curves collapse to a single (negative) inverse temperature $-\beta_{\mathrm{eff}}$ up to the edge of case (i). At $\eta=0$ we have $\beta_{\rm 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/\beta_{\rm 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\eta g^{2}\geq 1$) — an effective low-temperature dynamics is realized.

In case (ii), the ratio $R_{i}(\Delta E)\approx(\Delta E-g_{i}^{2})^{4}/(\Delta E+g_{i}^{2})^{4}$ tends to increase towards $1$ as $\Delta E$ grows. Typically, the available $\Delta 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 [math], implying that, if the parameters are adequately chosen, most jumps lie in regime (i). Second, since the rates are decreasing functions of $\left|\Delta E+g_{i}^{2}\right|$, 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 low-energy configurations, highlighting a clear bias of the dynamics towards cooling, as compared e.g., to the $\eta\gg 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 Zia and Schmittmann (2007) which we analyze for the loop formed in the configuration space of a two-spin system (see Fig. 3): $(\uparrow\uparrow)\rightarrow(\uparrow\downarrow)\rightarrow(\downarrow\downarrow)\rightarrow(\downarrow\uparrow)\rightarrow(\uparrow\uparrow)$. 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 $\eta\to\infty$ and also when $\eta\to 0$, signalling the emergence of the infinite-temperature dynamics. For different values $\eta$ 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. $766442$. I.L. gratefully acknowledges funding through the Royal Society Wolfson Research Merit Award.

I 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 $\rho(t)$ as $\dot{\rho}=\mathcal{L}_{0}(\rho)+\mathcal{L}_{1}(\rho)$, with $\mathcal{L}_{0}=\mathcal{L}-\mathcal{L}_{1}$ and $\mathcal{L}_{1}(\cdot)=-i\Omega[\sum_{i}\hat{\sigma}^{x}_{i},\cdot]$. Secondly, we assume the stationary state of $\mathcal{L}_{0}$ of the form $\rho_{\mathrm{stat}}=\sum_{\vec{\sigma}}p_{\vec{\sigma}}\rho_{\vec{\sigma}}\ket{\vec{\sigma}}\bra{\vec{\sigma}}$, where $\ket{\vec{\sigma}}=\{\sigma_{1},...,\sigma_{N}\}$, with $\sigma_{i}=\pm 1$ and $\hat{\sigma}_{i}^{z}\ket{\vec{\sigma}}=\sigma_{i}\ket{\vec{\sigma}}$, $p_{\vec{\sigma}}$ are a set of classical probabilities and $\rho_{\sigma}$ is the corresponding bosonic state, that we assume to be a gaussian state. Considering $\mathcal{L}_{1}$ perturbatively with respect to $\mathcal{L}_{0}$, and projecting the dynamics onto the stationary manifold of $\mathcal{L}_{0}$, we exploit the Nakajima-Zwanzig formalism to write the evolution of the spin as

[TABLE]

where $P\dot{\rho}^{\mathrm{spin}}(t)=\mathrm{Tr}_{B}[\dot{\rho}_{\mathrm{stat}}(t)]=\sum_{\vec{\sigma}}\dot{p}_{\vec{\sigma}}(t)\ket{\vec{\sigma}}\bra{\vec{\sigma}}$, $\mathrm{Tr}_{B}$ is the partial trace over the boson, and we have defined the spin-configuration dependent superoperators

[TABLE]

with $\mathcal{M}_{i}=\sum_{l\neq i}g_{l}\sigma_{l}$, and $\mathcal{D}_{\gamma}$, $\mathcal{D}_{\kappa}$ the dissipative terms representing cooling and heating, respectively. By projecting Eq.(S1) on a state $\ket{\vec{\sigma}^{\prime}}$, the dynamics reduces to the evolution of the classical probabilities ruled by a master equation whose general form is the following

[TABLE]

where $W_{\vec{\sigma}\rightarrow\vec{\sigma}^{\prime}}=\Omega^{2}\int_{0}^{+\infty}dt\sum_{j=\pm}\mathrm{Tr}_{B}\left[e^{\mathcal{V}_{\vec{\sigma},i}^{j}\tau}(\rho_{\vec{\sigma}})\right]$, is the transition rate for the switching $\vec{\sigma}\rightarrow\vec{\sigma}^{\prime}$ 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^{\mathcal{V}_{\vec{\sigma},i}^{\pm}}(\rho_{\vec{\sigma}})=(e^{\mathcal{V}_{\vec{\sigma},i}^{\pm,*}}\mathbb{1})\rho_{\vec{\sigma}}$, with $\mathcal{V}_{\vec{\sigma},i}^{\pm,*}$ the adjoint superoperator of $\mathcal{V}_{\vec{\sigma},i}^{\pm}$ and $\mathbb{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 $\mathbb{I}=\hat{D}(0)$, where $\hat{D}(\tau)=e^{\alpha(\tau)\hat{a}^{\dagger}-\beta(\tau)^{*}\hat{a}}e^{-\gamma(\tau)}$ with $\alpha(0)=\beta(0)=\gamma(0)=0$. We then verify that the displacement operator $\hat{D}(0)$ is mapped into the generalised one $\hat{D}(\tau)$ by applying the adjoint superoperator $\mathcal{V}_{\vec{\sigma},i}^{\pm}$: indeed, by considering the differential equation $\frac{d}{d\tau}\left[\hat{D}_{\vec{\sigma},i}^{\pm}(\tau)\right]=\mathcal{V}_{\vec{\sigma},i}^{\pm*}\left[\hat{D}_{\vec{\sigma},i}^{\pm}(\tau)\right]$ we obtain the solutions for the functions $\alpha^{+}_{\sigma_{i}}(\tau)=\alpha^{-}_{-\sigma_{i}}(\tau)$, $\beta^{+}_{\sigma_{i}}(\tau)=\beta^{-}_{-\sigma_{i}}(\tau)$, $\gamma^{+}_{\sigma_{i}}(\tau)=\gamma^{-}_{-\sigma_{i}}(\tau)$. For initial conditions $\alpha_{\sigma_{i}}^{\pm}(0)=\beta_{\sigma_{i}}^{\pm}(0)=\gamma_{\sigma_{i}}^{\pm}(0)=0$, we get

[TABLE]

where we have defined the dimensionless time $\tau=t\omega$, and $\eta=(\gamma-\kappa)/\omega$, $\theta=\kappa/\gamma\in[0,1)$, and

[TABLE]

The previous steps allow us to write the partial trace over the boson as $\mathrm{Tr}_{B}[e^{\mathcal{V}_{\vec{\sigma},i}^{\pm}}(\rho_{\vec{\sigma}})]=e^{-\gamma_{\sigma_{i}}^{\pm}(\tau)}\mathrm{Tr}_{B}[\hat{D}^{\pm}_{\vec{\sigma},i}(\tau)\rho_{\vec{\sigma}}]$. We recall that the bosonic state $\rho_{\vec{\sigma}}$ has been assumed to be a gaussian state. In this case, we recognise the quantity $\mathrm{Tr}_{B}[\hat{D}^{\pm}_{\vec{\sigma},i}(\tau)\rho_{\vec{\sigma}}]$ as the characteristic function $\chi^{\pm}_{\vec{\sigma},i}(\tau)$ of the state $\rho_{\vec{\sigma}}$. The expression of the characteristic function for a generic gaussian state $\rho_{G}$ reads

[TABLE]

where $\vec{\alpha}^{T}=(\alpha(\tau),\alpha^{*}(\tau))$, $\braket{\cdot}_{G}$ is the expectation value performed over the state $\rho_{G}$, and $\Sigma$ represent the covariance matrix which reads

[TABLE]

By applying the definition (S5), with expectation values of the operators obtained considering the lindblad operator $\mathcal{L}_{0}$, we get

[TABLE]

where $\mathcal{M}=\sum_{l}g_{l}\sigma_{l}$. Thus, the expression of the rate reads

[TABLE]

where $\Delta E_{i}~{}=g_{i}\sigma_{i}\sum_{l\neq i}g_{l}\sigma_{l}$ retains the dependence on the spin configuration.

II Approximate expression of the rate for $\eta$ small

For sufficiently small $\eta\lesssim 1$, the exponent appearing in $W_{\vec{\sigma}\rightarrow\vec{\sigma}^{\prime}}$ is dominated by

[TABLE]

implying that the integrand is quickly suppressed as $\tau$ grows. We therefore perform an expansion around $\tau=0$ which yields

[TABLE]

and

[TABLE]

From the first term in Eq. (S10) we see that the integrand is strongly suppressed on scales $\tau\gtrsim\sqrt{\eta}$. Noticing that in the Taylor expansion of $\tau+f(\tau)$ odd coefficients are finite, whereas even ones are $O(1/\eta)$ for $\eta\to 0$ and introducing the rescaled integration variable $z=\tau\sqrt{\eta}$ we see that higher orders are perturbations of order $O(\eta^{n+1/2}z^{2n+1},\eta^{n-1}z^{2n})$ and can be neglected. Similar considerations can be applied to $s(\tau)$, which can be therefore also approximated with its leading order $-\tau$. In the following, for simplicity we set $\omega=1$, remembering that our “energy” $\Delta E_{i}$ is actually measured by construction in units of $\omega^{2}$. Additionally, we introduce the shorthand

[TABLE]

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

[TABLE]

which can be integrated to give the closed expression

[TABLE]

It is worth remarking that the exponent can be rewritten as

[TABLE]

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(\Gamma_{i}\tau)$, whose series only coincides with the former up to $O(\tau^{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 $\Gamma_{i}$ of the order of $\sim 1/\sqrt{\eta}$. Since we wish to understand the behavior of the rates as functions of the energy difference $\Delta E_{i}$ without restrictions imposed by the other parameters (like $\eta$), we need to account for energies which exceed this range. To do this, we extract the asymptotic behavior of the rate for $\Gamma_{i}\to\infty$. We start by rewriting the integrand in $W$ as

[TABLE]

We now use the result that, if the function $A$ admits a small $\tau$ expansion

[TABLE]

then asymptotically in the limit $\Gamma_{i}\to\infty$ one finds

[TABLE]

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 ${\mathrm{Re}}[a_{n}i^{n+1}]$ equals $(-1)^{l+1}{\mathrm{Re}}[a_{2l+1}]$ if $n=2l+1$ is odd, and $(-1)^{l+1}{\mathrm{Im}}[a_{2l}]$ if $n=2l$ is even. For our function we find

[TABLE]

The leading behavior in the large $\Gamma_{i}$ limit is therefore determined by ${\mathrm{Re}}[a_{3}]$, implying

[TABLE]

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

[TABLE]

we find

[TABLE]

where $A=32g_{i}^{2}\pi(2\kappa+\eta)\eta^{-3}$.