Stochastic Operator Variance: an observable to diagnose noise and scrambling

TL;DR

This paper introduces the stochastic operator variance (SOV), a new observable to diagnose noise and scrambling in quantum systems, linking it to out-of-time-order correlators and demonstrating its behavior in a stochastic Lipkin-Meshkov-Glick model.

Contribution

The paper presents the SOV as a novel observable that characterizes noise effects and quantum scrambling, with analytical and numerical analysis in a stochastic model.

Findings

SOV obeys an uncertainty relation.

SOV dynamics are linked to quantum Lyapunov exponent.

SOV behavior is demonstrated in a stochastic LMG model.

Abstract

Noise is ubiquitous in nature, so it is essential to characterize its effects. Considering a fluctuating Hamiltonian, we introduce an observable, the stochastic operator variance (SOV), which measures the spread of different stochastic trajectories in the space of operators. The SOV obeys an uncertainty relation and allows finding the initial state that minimizes the spread of these trajectories. We show that the dynamics of the SOV is intimately linked to that of out-of-time-order correlators (OTOCs), which define the quantum Lyapunov exponent $\lambda$. Our findings are illustrated analytically and numerically in a stochastic Lipkin-Meshkov-Glick (sLMG) Hamiltonian undergoing energy dephasing.

Equations

$${\hat{H}}_t={\hat{H}}_0+\sqrt{2\gamma }{\xi }_t\hat{L},$$

$${\hat{U}}_{\mathrm{d}t}=e^{-i{\hat{H}}_0\mathrm{d}t-i\sqrt{2\gamma }\mathrm{d}{W}_{t-\mathrm{d}t}\hat{L}}.$$

$$\mathrm{d}{\hat{A}}_t={\mathcal{L}}^{†}[{\hat{A}}_t]\mathrm{d}t+i\sqrt{2\gamma }[\hat{L},{\hat{A}}_t]\mathrm{d}{W}_{t-\mathrm{d}t},$$

$$\frac{\mathrm{d}(\Delta {\hat{A}}_t^2)}{\mathrm{d}t}={\mathcal{L}}^{†}[\Delta {\hat{A}}_t^2]-2\gamma [\hat{L},⟨{\hat{A}}_t⟩{]}^2.$$

$$\Delta {\hat{A}}_t^2\equiv ⟨{\hat{A}}_t^2⟩-⟨{\hat{A}}_t{⟩}^2=e^{{\mathcal{L}}^{†}t}[{\hat{A}}^2]-(e^{{\mathcal{L}}^{†}t}[\hat{A}]{)}^2.$$

$$\mathrm{Tr}(\Delta {\hat{A}}_t^2{\rho }_0)\mathrm{Tr}(\Delta {\hat{B}}_t^2{\rho }_0)$$

$$\ge \frac{1}{4}\left(D_{+}^2(\hat{A},\hat{B})-D_{-}^2(\hat{A},\hat{B})\right),$$

$$\frac{1}{N}\frac{\mathrm{d}\mathrm{Tr}(\Delta {\hat{A}}_t^2)}{\mathrm{d}t}=-\frac{2\gamma }{N}\mathrm{Tr}([\hat{L},⟨{\hat{A}}_t⟩{]}^2).$$

$$C_t=\frac{1}{2\gamma N}\frac{\mathrm{d}\mathrm{Tr}(\Delta {\hat{A}}_t^2)}{\mathrm{d}t}\sim ϵe^{{\lambda }_{\text{textsc}q}t},$$

$$C_t=\sum _{m,n}(l_m-l_n{)}^2{e}^{-2\gamma (l_m-l_n{)}^{2}t}|A_{nm}{|}^2,$$

$${\hat{H}}_{\text{textsc}lmg}=\Omega {\hat{S}}_z-\frac{2}{N}{\hat{S}}_x^2,$$

$${\hat{H}}_t={\hat{H}}_{\text{textsc}lmg}(1+\sqrt{2\gamma }{\xi }_t).$$

$$|\Psi ⟩=\underset{t\to \infty }{\lim }|v_0(t)⟩,$$

$$\frac{\mathrm{d}}{\mathrm{d}t}|{\psi }_t⟩=-i{\hat{H}}_t|{\psi }_t⟩,$$

$${\hat{A}}_t^{\text{textsc}(h)}={\hat{U}}_t^{†}\hat{A}{\hat{U}}_t.$$

$$\frac{\mathrm{d}}{\mathrm{d}t}{\hat{A}}_t^{\text{textsc}(h)}$$

$$⟨\psi ⟩{\mathcal{T}}_{\to }{e}^{+i{\int }_0^{t_1}{\hat{H}}_{\tau }\mathrm{d}\tau }{\mathcal{T}}_{\to }{e}^{+i{\int }_{t_1}^{t_2}{\hat{H}}_{\tau }\mathrm{d}\tau }\hat{A}{\mathcal{T}}_{\leftarrow }{e}^{-i{\int }_{t_1}^{t_2}{\hat{H}}_{\tau }\mathrm{d}\tau }{\mathcal{T}}_{\leftarrow }{e}^{-i{\int }_0^{t_1}{\hat{H}}_{\tau }\mathrm{d}\tau }\psi ,$$

$$\mathrm{d}{\hat{A}}_t=i[\hat{H},\hat{A}]\mathrm{d}t+\sqrt{\eta }({\hat{c}}^{†}{\hat{A}}_t+{\hat{A}}_t\hat{c})\mathrm{d}{Y}_{t-\mathrm{d}t}+\sum _m\left({\hat{L}}_m{\hat{A}}_t{\hat{L}}_m-\frac{1}{2}\{{\hat{L}}_m^{†}{\hat{L}}_m,{\hat{A}}_t\}\right)\mathrm{d}t,$$

$${\hat{U}}_{dt}=e^{-i\hat{H}\mathrm{d}t-i\sqrt{2\gamma }\hat{L}\mathrm{d}{W}_{t-\mathrm{d}t}},$$

$$\frac{d}{dt}⟨{\hat{A}}_t{⟩}^2={\mathcal{L}}^{†}[⟨{\hat{A}}_t⟩]⟨{\hat{A}}_t⟩+⟨{\hat{A}}_t⟩{\mathcal{L}}^{†}[⟨{\hat{A}}_t⟩].$$

$$\frac{\mathrm{d}}{\mathrm{d}t}⟨{\hat{A}}_t{⟩}^2-{\mathcal{L}}^{†}[⟨{\hat{A}}_t{⟩}^2]$$

$$\displaystyle=-2\gamma\bigg{(}\hat{L}\langle\hat{A}_{t}\rangle^{2}\hat{L}+\langle\hat{A}_{t}\rangle\hat{L}^{2}\langle\hat{A}_{t}\rangle-\hat{L}\langle\hat{A}_{t}\rangle\hat{L}\langle\hat{A}_{t}\rangle-\langle\hat{A}_{t}\rangle\hat{L}\langle\hat{A}_{t}\rangle\hat{L}\bigg{)}$$

$$=-2\gamma [\hat{L},⟨{\hat{A}}_t⟩{]}^2,$$

$$\frac{\mathrm{d}}{\mathrm{d}t}\Delta {\hat{A}}_t^2={\mathcal{L}}^{†}[\Delta {\hat{A}}_t^2]-2\gamma [\hat{L},⟨{\hat{A}}_t⟩{]}^2.$$

$$\left(\begin{array}[]{cc}\Delta\hat{A}_{t}^{2}&\Delta\widehat{AB}_{t}\\ \Delta\widehat{AB}_{t}^{\dagger}&\Delta\hat{B}^{2}_{t}\end{array}\right)\geq 0,$$

$$\left(\begin{array}[]{cc}\langle\Delta\hat{A}_{t}^{2}\rangle_{\rho}&\langle\Delta\widehat{AB}_{t}\rangle_{\rho}\\ \langle\Delta\widehat{AB}_{t}^{\dagger}\rangle_{\rho}&\langle\Delta\hat{B}^{2}_{t}\rangle_{\rho}\end{array}\right)\geq 0.$$

$$⟨\Delta {\hat{A}}_t^2{⟩}_{\rho }\ge 0,⟨\Delta {\hat{B}}_t^2{⟩}_{\rho }\ge 0,⟨\Delta {\hat{A}}_t^2{⟩}_{\rho }⟨\Delta {\hat{B}}_t^2{⟩}_{\rho }-|⟨\Delta {\widehat{AB}}_t{⟩}_{\rho }{|}^2\ge 0.$$

$$⟨\Delta {\hat{A}}_t^2{⟩}_{\rho }⟨\Delta {\hat{B}}_t^2{⟩}_{\rho }\ge |⟨\Delta {\widehat{AB}}_t{⟩}_{\rho }{|}^2.$$

$$|⟨\Delta {\widehat{AB}}_t{⟩}_{\rho }{|}^2$$

$$\displaystyle=\frac{1}{4}\bigg{(}\langle e^{\mathcal{L}^{\dagger}t}[\{\hat{A},\hat{B}\}]\rangle_{\rho}^{2}-\langle e^{\mathcal{L}^{\dagger}t}[[\hat{A},\hat{B}]]\rangle_{\rho}^{2}-2\langle e^{\mathcal{L}^{\dagger}t}[\{\hat{A},\hat{B}\}]\rangle_{\rho}\langle\{\hat{A}_{t},\hat{B}_{t}\}\rangle_{\rho}+2\langle e^{\mathcal{L}^{\dagger}t}[[\hat{A},\hat{B}]]\rangle_{\rho}\langle[\hat{A}_{t},\hat{B}_{t}]\rangle_{\rho}$$

$$\displaystyle\hskip 11.38109pt+\langle\{\hat{A}_{t},\hat{B}_{t}\}\rangle_{\rho}-\langle[\hat{A}_{t},\hat{B}_{t}]\rangle_{\rho}^{2}\bigg{)}$$

$$=\frac{1}{4}{\left(⟨e^{{\mathcal{L}}^{†}t}[\{\hat{A},\hat{B}\}]{⟩}_{\rho }-⟨\{{\hat{A}}_t,{\hat{B}}_t\}{⟩}_{\rho }\right)}^2-\frac{1}{4}{\left(⟨e^{{\mathcal{L}}^{†}t}[[\hat{A},\hat{B}]]{⟩}_{\rho }-⟨[{\hat{A}}_t,{\hat{B}}_t]{⟩}_{\rho }\right)}^2$$

$$=\frac{1}{4}\left(D_{+}^2(\hat{A},\hat{B})-D_{-}^2(\hat{A},\hat{B})\right),$$

$$\mathrm{Tr}(\Delta {\hat{A}}_t^2\rho )\mathrm{Tr}(\Delta {\hat{B}}_t^2\rho )\ge \frac{1}{4}\left(D_{+}^2(\hat{A},\hat{B})-D_{-}^2(\hat{A},\hat{B})\right).$$

Full text

Any realistic quantum system inevitably experiences some fluctuations induced by its surrounding environment. Current quantum technologies are limited by the action of this noise, motivating the pragmatic focus on Noisy-Intermediate-Scale Quantum (NISQ) devices [1, 2]. In any experimental setting, tunable parameters such as Hamiltonian coupling constants may exhibit fluctuations due to interactions with the environment [3, 4, 5]. In this context, the dynamics of an ensemble of noisy realizations can be described in terms of the noise-averaged density matrix, which evolves according to a master equation describing non-unitary evolution [3, 4, 6, 7]. Alternatively, noise can be utilized as a resource for the quantum simulation of open systems [4]. The study of fluctuations in noisy quantum systems is also connected to free probability [8, 9, 10].

The quest for understanding noise in chaotic systems has recently led to a flurry of activities exploring the signatures of quantum chaos when the dynamics is no longer unitary [11, 12, 6, 13, 14, 15, 16, 17, 18]. Out-of-time-order correlators (OTOCs) offer an important diagnostic tool, which was initially proposed in the theory of superconductivity [19]. Their use experienced renewed interest in defining a quantum analog of the Lyapunov exponent [20, 21], which measures the exponential sensitivity to the initial conditions in chaotic systems and is universally bounded by the system’s temperature [22]. The existence of a positive Lyapunov exponent classically is a necessary but not sufficient condition for the system to be chaotic—e.g. [23, 24]. Similarly, the exponential growth of the OTOC is not a sufficient signature for quantum chaos but rather indicates scrambling [25, 26]. OTOCs have been studied experimentally [27, 28, 29, 30, 31, 32] and in open systems where their evolution is changed by dissipation [33, 34, 35, 21, 36].

In this Letter, we consider the dynamics generated by a stochastic Hamiltonian and go beyond an average description of the state by introducing the Stochastic Operator Variance (SOV), an observable that characterizes the spread of stochastic trajectories of any operator. This notion is directly relevant to experiments, particularly in NISQ devices subject to various sources of noise. We compute the evolution of the SOV and show that it obeys a generalized uncertainty relation. It allows to identify the initial state that minimizes the deviation from Lindblad dynamics at long times. Surprisingly, we also find that the SOV evolution relates to that of an OTOC, which connects fluctuations of the system with scrambling. This is pictorially represented in Fig. 1. We illustrate our results in a stochastic generalization of the Lipkin-Meshkov-Glick (sLMG) model. The LMG model [37] describes an Ising spin chain with infinite range interactions and exhibits scrambling from an unstable fixed point [38, 39, 25]. It can be realized experimentally with trapped ions [40] and is amenable to dynamical control techniques such as shortcuts to adiabaticity [41]. Considering fluctuations in the energy scale, we compute the SOV in this model, show its connection to the OTOC, and use it to characterize the Lyapunov exponent in the classical limit.

The stochastic operator variance (SOV).— Let us consider a system evolving under a Hermitian Hamiltonian $\hat{H}_{0}$ and subject to classical noise $\xi_{t}$ modulating the coupling constant of a Hermitian operator $\hat{L}$, i.e.,

[TABLE]

where $\gamma$ measures the coupling strength between the system and noise—we set $\hbar=1$. The stochastic process is taken as real Gaussian white noise, that is, $\langle\xi_{t}\rangle=0$ and $\langle\xi_{t}\xi_{t^{\prime}}\rangle=\delta(t-t^{\prime})$. We will represent the stochastic averages by $\langle\bullet\rangle$; quantum expectation values will be written explicitly taking the trace over the state density matrix, $\mathrm{Tr}(\bullet\rho)$. We introduce the Wiener process $\mathrm{d}W_{t}\equiv\xi_{t}\mathrm{d}t$, which is convenient to deal with the formal treatment of stochastic differential equations (SDEs). It obeys Itō’s rules $\mathrm{d}W_{t}^{2}=\mathrm{d}t$ with vanishing higher-order terms, $\mathrm{d}t\,\mathrm{d}W_{t}=\mathrm{d}W_{t}^{k+1}=\mathrm{d}t^{k}=0\;\forall k>1$ [42]. Over a time increment $\mathrm{d}t$, the evolution of a Hermitian operator $\hat{A}$ in the Heisenberg picture is $\hat{A}_{t-\mathrm{d}t}=\hat{U}_{\mathrm{d}t}^{\dagger}\hat{A}_{t}\hat{U}_{\mathrm{d}t}$ (see App. A), since operators in the Heisenberg picture evolve backward in time [43, 44]. The associated propagator reads (see App. A)

[TABLE]

Expanding the propagator according to those rules and introducing the backward differential $\mathrm{d}\hat{A}_{t}=\hat{A}_{t-\mathrm{d}t}-\hat{A}_{t}$ [43, 44] yields the SDE for the evolution of $\hat{A}_{t}$ under the stochastic Hamiltonian (1),

[TABLE]

where $\mathcal{L}^{\dagger}[\hat{A}_{t}]=i[\hat{H}_{0},\hat{A}_{t}]-\gamma[\hat{L},[\hat{L},\hat{A}_{t}]]$ is the adjoint Lindbladian. Upon averaging (3), all the linear terms in $\mathrm{d}W_{t-\mathrm{d}t}$ vanish [42] and we find that the noise-averaged operator evolves with an adjoint Lindblad equation $\mathrm{d}_{t}\langle\hat{A}_{t}\rangle=\mathcal{L}^{\dagger}[\langle\hat{A}_{t}\rangle]$. This corresponds to the standard evolution of an observable in an open quantum system with a Hermitian jump operator $\hat{L}$ [45]. The formalism described so far has been introduced in [3] and used in [4] to engineer long-range and many-body interactions. Here, we focus on the stochastic variance of an observable.

In order to find the variance, the second stochastic moment $\langle\hat{A}_{t}^{2}\rangle$ is needed. Considering the reference operator to be $\hat{A}^{2}$ instead of $\hat{A}$, one finds that its evolution follows the Lindblad equation, $\mathrm{d}_{t}\langle\hat{A}_{t}^{2}\rangle=\mathcal{L}^{\dagger}[\langle\hat{A}_{t}^{2}\rangle]$. Recall that the average is over realizations of the noise and that $\langle\hat{A}_{t}^{2}\rangle$ is still an operator acting on the Hilbert space. Subtracting $\mathrm{d}_{t}\langle\hat{A}_{t}\rangle^{2}{=\langle\hat{A}_{t}\rangle\mathcal{L}^{\dagger}[\langle\hat{A}_{t}\rangle]+\mathcal{L}^{\dagger}[\langle\hat{A}_{t}\rangle]}\langle\hat{A}_{t}\rangle=\mathcal{L}^{\dagger}[\langle\hat{A}_{t}\rangle^{2}]+2\gamma[\hat{L},\langle\hat{A}_{t}\rangle]^{2}$ from both sides (see App. A), we find the evolution of the SOV, defined as $\Delta\hat{A}^{2}_{t}=\langle\hat{A}_{t}^{2}\rangle-\langle\hat{A}_{t}\rangle^{2}$, to be given by

[TABLE]

The SOV $\Delta\hat{A}_{t}^{2}$ is an observable that characterizes the deviation of any (stochastic) operator $\hat{A}_{t}$ from the noise-averaged operator in a stochastic evolution governed by the Hamiltonian (1)—see Fig. 1 for a scheme and App. D for a quantitative illustration. Although its equation of motion depends on out-of-time-order terms like $\hat{L}\langle\hat{A}_{t}\rangle\hat{L}\langle\hat{A}_{t}\rangle$, it can easily be computed from the evolution of $\langle\hat{A}_{t}\rangle$ and $\langle\hat{A}_{t}^{2}\rangle$. Indeed, the SOV evolves as

[TABLE]

The average $\langle\hat{A}_{t}\rangle$ follows the standard Lindblad dissipative dynamics from which the individual trajectories deviate as dictated by the variance $\Delta\hat{A}_{t}^{2}$. Since the SOV is Hermitian, it is an observable. It differs from the quantum variance of an operator $\hat{A}$ over a state $\rho_{t}$, commonly defined as $\mathrm{Var}(\hat{A},\rho_{t})=\mathrm{Tr}(\hat{A}^{2}\rho_{t})-\mathrm{Tr}(\hat{A}\rho_{t})^{2}$. For an initially pure state, this variance reads $\mathrm{Var}(\hat{A}_{t},\psi_{0})=\braket{\psi_{0}}{e^{\mathcal{L}^{\dagger}t}[\hat{A}^{2}]}{\psi_{0}}-\braket{\psi_{0}}{e^{\mathcal{L}^{\dagger}t}[\hat{A}]}{\psi_{0}}^{2}$. Such quantum variance cannot be obtained from the expectation value of an observable since it is nonlinear in the state $\ket{\psi_{0}}$. The difference between the SOV (5) evaluated on a pure state and the quantum variance reads $\braket{\psi_{0}}{\Delta\hat{A}_{t}^{2}}{\psi_{0}}-\mathrm{Var}(\hat{A}_{t},\psi_{0})=\braket{\psi_{0}}{e^{\mathcal{L}^{\dagger}t}[\hat{A}]\hat{\mathbb{Q}}e^{\mathcal{L}^{\dagger}t}[\hat{A}]}{\psi_{0}}$, where $\hat{\mathbb{Q}}=\mathbb{1}-\ket{\psi_{0}}\bra{\psi_{0}}$ is a projection operator on the complementary subspace of $\ket{\psi_{0}}\bra{\psi_{0}}$. Therefore, the SOV contains the information of the quantum variance plus the contribution from the projector over the complementary subspace, as the term governing the recombination of decay products in unstable systems [46].

A related notion of variance—without the stochastic interpretation—has been introduced in the theory of positive definite matrices [47], which provides us with tools to formally characterize the SOV. First, the SOV is positive semidefinite, $\Delta\hat{A}_{t}^{2}\geq 0$. This can be shown from Kadison’s inequality [48] which ensures that $e^{\mathcal{L}^{\dagger}t}[\hat{A}^{2}]\geq(e^{\mathcal{L}^{\dagger}t}[\hat{A}])^{2}$ for a positive and unital map, $e^{\mathcal{L}^{\dagger}t}[\mathbb{1}]=\mathbb{1}$. Second, our formalism leads to a generalization of the Robertson-Schrödinger uncertainty principle [49, 50]. Indeed, considering two operators $\hat{A},\;\hat{B}$ and their SOVs, we show in App. A that

[TABLE]

where $\Delta\widehat{AB}_{t}=e^{\mathcal{L}^{\dagger}t}[\hat{A}\hat{B}]-e^{\mathcal{L}^{\dagger}t}[\hat{A}]e^{\mathcal{L}^{\dagger}t}[\hat{B}]$ is the stochastic operator covariance of $\hat{A}$ and $\hat{B}$. The quantity D_{\eta}(\hat{A},\hat{B})\equiv\mathrm{Tr}\Big{(}e^{\mathcal{L}^{\dagger}t}\big{(}[\hat{A},\hat{B}]_{\eta}\big{)}\rho_{0}-[e^{\mathcal{L}^{\dagger}t}(\hat{A}),e^{\mathcal{L}^{\dagger}t}(\hat{B})]_{\eta})\rho_{0}\Big{)} measures the difference between evolving the commutator ($\eta=-$) and anti-commutator ($\eta=+$) as a whole or each separately (see App. A). Note that $D_{+}$ is purely real and $D_{-}$ is purely imaginary, so $D_{+}^{2}-D_{-}^{2}\geq 0$. For $\hat{A}=\hat{B}$, the uncertainty relation is saturated and Eq. (6) becomes an equality.

Bi-partite interpretation—Many quantum properties, such as entanglement or scrambling, are best understood in a bipartite system [51, 52, 53, 54]. The SOV can actually be interpreted analogously: consider a doubled Hilbert space $\mathscr{H}\otimes\mathscr{H}$, with an operator $\hat{A}$ living in each of the copies of the Hilbert space, denoted as $\mathscr{H}_{1}$ and $\mathscr{H}_{2}$. Using the swap operator $\mathbb{S}$ [53, 54], which introduces an interaction between the Hilbert spaces and is defined as $\mathbb{S}(\ket{i}_{1}\ket{j}_{2})=\ket{j}_{1}\ket{i}_{2}$, the product between operators can be interpreted as an operation over a doubled Hilbert space, namely \hat{X}\hat{Y}=\mathrm{Tr}_{\mathscr{H}_{2}}\big{(}(\hat{X}\otimes\hat{Y})\mathbb{S}\big{)}—see details in App. C. Thus, the first term in the SOV (5) can be interpreted as first letting the operators $\hat{A}$ interact to form $\hat{A}^{2}$ in the bipartite system and then letting that operator evolve under the dissipative evolution of a single bath, $e^{\mathcal{L}^{\dagger}t}[\hat{A}^{2}]$. By contrast, the second term corresponds to letting each of the uncoupled systems evolve with their own bath and then making them interact at time $t$. As expected, the difference between these terms is positive since $\Delta\hat{A}_{t}^{2}\geq 0$, and the first protocol always suffers less decoherence.

The SOV-OTOC connection—Remarkably, the expectation value of (4) over the completely-mixed state, $\hat{\rho}=\hat{\mathbb{1}}/N$, gives a dissipative version of the OTOC, namely

[TABLE]

OTOCs are typically defined from two operators as $C_{t}=-\mathrm{Tr}([\hat{B}_{0},\hat{A}_{t}]^{2})/N$, and measure the exponential sensitivity on initial conditions in quantum chaotic systems [20]. Indeed, in a quantum system with scrambling, one expects $C_{t}\sim\epsilon e^{\lambda_{\textsc{q}}t}$ in the time window $t_{s}\ll t\ll t_{\textsc{e}}$ between the saturation time of two-point functions, $t_{s}\sim 1/\lambda_{\textsc{q}}$, and that of the OTOC, known as the Ehrenfest time, $t_{\textsc{e}}\sim\ln(\hbar^{-1})/\lambda_{\textsc{q}}$ [22]. The main difference in our setting is that the evolved operator follows dissipative dynamics, $e^{\mathcal{L}^{\dagger}t}[\hat{A}]$, instead of unitary evolution, $e^{i\hat{H}t}\hat{A}e^{-i\hat{H}t}$. The connection between this OTOC and the SOV is pictorially shown in Fig. 1. It can be used to compute the Lyapunov exponent through

[TABLE]

where the exponential behavior holds only in systems with scrambling over the appropriate period $t_{s}\ll t\ll t_{\textsc{e}}$. The SOV-OTOC connection is complementary to the optimal-path approach to study chaos in continuously monitored systems [55]. Note that the SOV $\Delta\hat{A}_{t}^{2}$ is an observable constructed from the knowledge of the evolution under different noise realizations. The classical limit of the above equation is similar, the only difference being that $\Delta A_{t}^{2}$ becomes a function of time and the trace an average over a region of phase space [56].

The short-time decay of the OTOC, $C_{t}\sim C_{0}e^{-t/\tau_{D}}$, is characterized by the dissipation time \tau_{D}=\big{(}2\gamma\mathrm{Tr}([\hat{L},[\hat{L},\hat{A}]]^{2})/(C_{0}N)\big{)}^{-1} and the initial value $C_{0}=\mathrm{Tr}([\hat{L},\hat{A}]^{2})/N$. Interestingly, this dissipation time is related to the Hilbert-Schmidt norm of the dissipator acting on the initial operator (see App. B). It is analogous to the decoherence rate found in [6] but now obtained for operators in the Heisenberg picture. When $[\hat{H}_{0},\hat{L}]=0$, the operators share a common eigenbasis, and using $\hat{L}\ket{n}=l_{n}\ket{n}$, we can write the dissipative OTOC (7) as

[TABLE]

where $A_{nm}=\braket{n}{\hat{A}}{m}$. From this expression, the two exponentially decaying regimes reported in [57] appear from the largest and smallest eigenvalue differences governing the short- and long-time dynamics, respectively.

Note that, for simplicity, we have focused on the infinite temperature OTOC. However, the connection to OTOC (7) can be generalized to unregularized thermal OTOC [58] by tracing (4) over a thermal state $\rho_{\beta}=e^{-\beta\hat{H}}/\mathrm{Tr}(e^{-\beta\hat{H}})$. In this case, the SOV-OTOC relation also involves an additional expectation value of the Lindbladian $\mathrm{Tr}(\mathcal{L}[\Delta\hat{A}_{t}^{2}]\rho_{\beta})$, which is non-zero in general. Fidelity OTOCs [59] can be obtained by taking the jump operator to be a projector over a pure state. Regularized OTOCs [22] can also be obtained by modifying the jump operator as $\hat{L}\rightarrow\rho_{\beta}^{1/4}\hat{L}\rho_{\beta}^{1/4}$ and tracing over the completely mixed state.

We next illustrate our findings in the Lipkin-Meshkov-Glick (LMG) model subject to energy dephasing and characterize its Lyapunov exponent using the SOV-OTOC connection.

Stochastic Lipkin-Meshkov-Glick (sLMG) model.— The LMG model describes the collective motion of $N$ identical two-level systems fully connected to each other with the same coupling strength [37]. Its quantum Hamiltonian reads

[TABLE]

where $\Omega$ is the frequency in units of the coupling strength, and $\hat{S}_{j}$ are the general spin operators of dimension $2S+1$. We stay in the sector $S=N/2$. Since the total spin $\hat{\mathbf{S}}^{2}=\hat{S}_{x}^{2}+\hat{S}_{y}^{2}+\hat{S}_{z}^{2}$ commutes with the spin operators, $[\hat{S}_{j},\hat{\mathbf{S}}^{2}]=0$, the total angular momentum is conserved. Due to time-translational symmetry, energy is conserved, and since there is only one degree of freedom, this model is integrable. If this continuous symmetry is broken by periodic kicks in $\hat{S}_{x}^{2}$, the model turns into the known kicked top [11].

Here, we break time-translational symmetry by adding noise in the energy scale and consider

[TABLE]

This leads to dephasing in the energy eigenbasis at the ensemble level. The evolution of the SOV instantaneous (ordered) eigenvalues satisfying $\Delta\hat{A}_{t}^{2}\ket{v_{k}(t)}=\Lambda_{k}(t)\ket{v_{k}(t)}$, for this model is shown in Fig. 2(a). In analogy with the theory of quantum transport [60, 61], we find diffusive modes in which $\Lambda_{k}(t)\sim t$ and superdiffusive modes in which $\Lambda_{k}(t)\sim t^{3/2}$. In the more general case $[\hat{H}_{0},\hat{L}]\neq 0$, one also finds ballistic modes $\Lambda_{k}(t)\sim t^{2}$, as we detail in App. D. Based on the eigenvalues, one can find a state which minimizes the SOV at long times

[TABLE]

which corresponds to the steady state of the eigenvector with the smallest deviation, and is thus the state minimally affected by the noise. The expectation value over this state is shown in Fig. 2(a) (black dash-dotted) for the sLMG model, and we verify that it minimizes the spread at long times. Fig. 2(b) shows the evolution of the dissipative OTOC, computed from the SOV using Eq. (8). We observe several exponential decays, as apparent from (9) and in agreement with [57]. This illustrates how the SOV can be used to obtain the dissipative OTOC.

The classical limit of the Hamiltonian (10) is obtained by taking its expectation value over SU(2) coherent states $\ket{\zeta}$ [62] (see App. E) in the thermodynamic limit, $N\rightarrow\infty$. We introduce the canonical variables $Q$ and $P$ as $\zeta=\frac{Q-iP}{\sqrt{4-(Q^{2}+P^{2})}}$ [38, 39], which yields $H_{\textsc{lmg}}=\lim_{S\rightarrow\infty}\braket{\zeta}{\hat{H}_{\textsc{lmg}}}{\zeta}/S=\tfrac{\Omega}{2}P^{2}+\left(\tfrac{\Omega}{2}-1\right)Q^{2}+\tfrac{1}{4}(Q^{2}P^{2}+Q^{4}),$ where the terms of $\mathcal{O}(1/N)$ are neglected (see App. E). This model is integrable and exhibits an unstable fixed point at the origin, $Q^{*}=P^{*}=0$, for $0<\Omega<2$. Since scrambling originates from an unstable point, it is already present in the semiclassical limit [25, 38, 39].

Here, we consider the classical equivalent of (11), namely, $H_{t}=H_{\textsc{lmg}}(1+\sqrt{2\gamma}\xi_{t})$. The evolution of the noise-averaged observable displays the classical analog of energy dephasing, namely $\partial_{t}\langle A_{t}\rangle=-\{H_{\textsc{lmg}},\langle A_{t}\rangle\}_{\textsc{p}}+2\gamma\{H_{\textsc{lmg}},\{H_{\textsc{lmg}},\langle A_{t}\rangle\}_{\textsc{p}}\}_{\textsc{p}}$, where $\{f,g\}_{\textsc{p}}$ denotes the Poisson bracket of $f$ and $g$. We characterize the Lyapunov exponent using three complementary methods: Analytically, the approach proposed by van Kampen [63, 64] yields the Lyapunov for the sLMG as (see App. E) $\lambda^{(1)}=\sqrt{2\Omega-\Omega^{2}}-\gamma(2\Omega-\Omega^{2}).$ (ii) Second, numerically. The standard, classical definition of the average Lyapunov exponent gives $\lambda^{(2)}=\langle\lim_{t\rightarrow\infty}\ln(\ell_{t}/\ell_{0})/t\rangle,$ where $\ell_{t}^{2}=(Q_{t}{-}Q_{t}^{\prime})^{2}+(P_{t}{-}P_{t}^{\prime})^{2}$ is the distance between two initially close trajectories evolving with the same realization of the noise, which is found solving the stochastic Hamilton’s equations with a stochastic Runge-Kutta method [65] (see App. F). (iii) Finally, our formalism gives the Lyapunov from Eq. (7), by taking as observable the position $A_{t}=Q_{t}$, the jump operator being the Hamiltonian itself, $L=H_{0}$. Namely, $\lambda^{(3)}=\lim_{t\rightarrow\infty}\frac{1}{2t}\ln\left(\mathrm{d}_{t}\Delta Q_{t}^{2}/\epsilon\right)$.

Figure 3(a) shows the Lyapunov exponent $\lambda$ obtained from the above methods as a function of $\Omega$ for different noise strengths. We verify that the three methods are in good agreement up to numerical errors. For $\gamma=0$ (black line), we recover the behavior of LMG with $\lambda>0$ in the double well (DW) phase and $\lambda=0$ in the single well (SW) phase. Introducing a weak stochastic perturbation with $\gamma$, the Lyapunov exponent becomes smaller in the DW phase—trajectories diverge more slowly—while the SW phase acquires a positive $\lambda$. Increasing the level of noise causes the Lyapunov exponent in the DW phase to decrease and even reach negative values, where trajectories converge exponentially, while the value of $\lambda$ in the SW phase increases. This rich behavior is summarized in the phase diagram presented in Fig. 3(b). The latter further shows that, under the application of strong noise, the origin $(Q^{*},P^{*})=(0,0)$ is stable in the DW phase—trajectories converge exponentially to it—while it is unstable in the SW phase—trajectories diverge exponentially from it. Therefore the sLMG in the DW phase shows a noise-induced transition to stability [66], analogous to the stabilization seen for periodic driving [67].

In summary, we have introduced the stochastic operator variance and shown it is a valuable tool to study quantum systems driven by noise. We have shown that this observable obeys an uncertainty relation (6), and that it can be used to identify the state (12) which minimizes deviation from the Lindblad dissipative dynamics. We have provided a bipartite interpretation of the SOV. In addition, our results unveil a SOV-OTOC connection (8), which provides an operational protocol harnessing noise as a resource to probe OTOC and extract the Lyapunov exponent in noisy quantum chaotic systems. To illustrate our results, we introduced a stochastic generalization of LMG model and characterized its behavior in the quantum and classical realms. Our results provide the means to elucidate the fate of quantum chaos in noisy systems and benchmark NISQ devices.

Acknowledgements.— We thank Niklas Hörnedal, Howard Wiseman, Andrew Jordan, Federico Roccati, Federico Balducci, and Ruth Shir for insightful discussions and comments on the manuscript. This research was partly funded by the Luxembourg National Research Fund (FNR), Attract grant 15382998, and the John Templeton Foundation (Grant 62171). The opinions expressed in this publication are those of the authors and do not necessarily reflect the views of the John Templeton Foundation. For the purpose of open access, the authors have applied a Creative Commons Attribution 4.0 International (CC BY 4.0) license to any Author Accepted Manuscript version arising from this submission.

Appendix A Stochastic evolution: details on the main derivation

A.1 On the Heisenberg picture of an explicitly time-dependent Hamiltonian

We consider an explicitly time-dependent Hamiltonian $\hat{H}_{t}$, which in general does not commute with itself at different times, $[\hat{H}_{t},\hat{H}_{t^{\prime}}]\neq 0$. It governs the dynamics of a state through the Schrödinger equation

[TABLE]

whose solution $\ket{\psi_{t}}=\hat{U}_{t}\ket{\psi_{0}}=\mathcal{T}_{\leftarrow}e^{-i\int_{0}^{t}\hat{H}_{\tau}\mathrm{d}\tau}\ket{\psi_{0}},$ involves the propagator from time [math] to $t$, $\hat{U}_{t}$, where $\mathcal{T}_{\leftarrow}$ is the chronological time-ordering operator [45]. An (explicitly time-independent) operator $\hat{A}$ in the Heisenberg picture then evolves as

[TABLE]

Differentiating this evolution gives Heisenberg’s equation for the evolution of $\hat{A}_{t}$ as

[TABLE]

where $\hat{H}^{\textsc{(h)}}_{t}=\hat{U}_{t}^{\dagger}\hat{H}_{t}\hat{U}_{t}\neq\hat{H}_{t}$ is the Hamiltonian in the Heisenberg picture with respect to itself, which is not equal to $\hat{H}_{t}$ since it does not commute with itself at the different times $t^{\prime}<t$ included in $\hat{U}_{t}$. Note that the Heisenberg equation (15) is recovered when we introduce the Hamiltonian in the Heisenberg picture with respect to itself $\hat{H}^{\textsc{(h)}}_{t}$, but that its solution (14) only requires the bare Hamiltonian $\hat{H}_{t}$ in the propagator $\hat{U}_{t}$.

We bring attention to the fact that Heisenberg operators evolve “backwards in time”. This can be illustrated by splitting the evolution over two different times $0\leq t_{1}\leq t_{2}$. The expectation value of any general operator $\hat{A}$ then reads

[TABLE]

where $\mathcal{T}_{\rightarrow}$ denotes anti-chronological time-ordering. While states evolve forwards in time, $0\rightarrow t_{1}\rightarrow t_{2}$, the equivalent evolution at the level of operators is for them to evolve backwards in time, $t_{2}\rightarrow t_{1}\rightarrow 0$. Of course this is just a result of the Heisenberg representation and does not correspond to a backwards time evolution.

A.2 Heisenberg evolution with a stochastic Hamiltonian

Adjoint stochastic master equations (SMEs) have been studied in [43, 44] for continuous homodyne detection and read

[TABLE]

where we use the backwards differential $\mathrm{d}\hat{A}_{t}\equiv\hat{A}_{t-\mathrm{d}t}-\hat{A}_{t}$. $\eta$ is the efficiency of the detector, $\hat{c}$ is the measurement observable, $\hat{L}_{m}$ are the Lindblad operators describing coupling to a general bath, and $\mathrm{d}Y_{t}$ is the measurement record. $\mathrm{d}Y_{t}$ is a stochastic process. Note that this equation does not preserve trace. When $\mathrm{d}Y_{t}$ represents white noise without drift, i.e. $\langle\mathrm{d}Y_{t}\rangle=0$, the average over the noise recovers the standard adjoint master equation for the average operator $\langle\hat{A}_{t}\rangle$.

We draw attention to the fact that the measurement record appears as $\mathrm{d}Y_{t-\mathrm{d}t}$. To apply Itō’s rules, we need to evaluate the noise at the beginning of the time interval. In the Schrödinger picture, this corresponds to evaluating the noise at time $t$ in the interval $[t,t+\mathrm{d}t]$. But since we work in the Heisenberg picture we propagate the operator in the interval $[t-\mathrm{d}t,t]$ and we have to evaluate the noise at $t-\mathrm{d}t$. The propagator over $\mathrm{d}t$ thus reads

[TABLE]

which gives the SME (3). Note that this equation is trace preserving. So although in the case of driftless white noise $\langle\mathrm{d}Y_{t}\rangle=0$, the SME’s (3) and (17) give the same results at the average level, they describe different dynamics at the level of single trajectories.

A.3 Finding the equation of motion for the SOV

The equation of motion for the second stochastic moment simply reads $\frac{d}{dt}\langle\hat{A}^{2}_{t}\rangle=\mathcal{L}^{\dagger}[\langle\hat{A}^{2}_{t}\rangle]$, and that for the first moment squared reads

[TABLE]

In the latter, we want to introduce the term $\mathcal{L}^{\dagger}[\langle\hat{A}_{t}\rangle^{2}]$ in order to recover $\mathcal{L}^{\dagger}[\Delta\hat{A}^{2}_{t}]$ through linearity of the Lindbladian superoperator. The difference reads

[TABLE]

which yields the equation of motion (4) for the SOV given in the main text

[TABLE]

A.4 The SOV uncertainty principle

The theory of positive matrices (see pg. 75 [47]) tells us that given a completely positive unital map, like $e^{\mathcal{L}^{\dagger}t}[\bullet]$, for all operators $\hat{A},\;\hat{B}$, the following matrix is positive semidefinite

[TABLE]

i.e. the $2N\times 2N$ matrix has to be positive semidefinite. While this statement is powerful, it does not give any straightforward inequality. For this reason, we take the expectation value of each of the blocks over a certain state $\rho$, denoted $\langle\bullet\rangle_{\rho}=\mathrm{Tr}(\bullet\rho)$. This allows recasting the inequality for a 2x2 matrix,

[TABLE]

Let us now apply Sylvester’s criterion: a matrix is positive semi-definite if and only if all its principal minors, i.e. the determinants of the matrix in which we delete the columns and rows with the same index, are non-negative. Doing so, we can recast the complete positivity condition into

[TABLE]

The two first inequalities above are fulfilled because the SOV is positive semidefinite, the last condition yields an analog of the Schwarz inequality

[TABLE]

Let us now expand the r.h.s. of the inequality. Writing $\hat{A}\hat{B}=\frac{1}{2}(\{\hat{A},\hat{B}\}+[\hat{A},\hat{B}])$ and $\hat{B}\hat{A}=\frac{1}{2}(\{\hat{A},\hat{B}\}-[\hat{A},\hat{B}])$, we obtain

[TABLE]

In the last line, we introduced $D_{\eta}(\hat{A},\hat{B})=\langle e^{\mathcal{L}^{\dagger}t}([\hat{A},\hat{B}]_{\eta})\rangle_{\rho}-\langle[\hat{A}_{t},\hat{B}_{t}]_{\eta}\rangle_{\rho}$, where $[\hat{X},\hat{Y}]_{\eta}=\hat{X}\hat{Y}+\eta\hat{Y}\hat{X}$ with $\eta=\pm 1$ is the generalized commutator. Therefore the SOV uncertainty principle reads

[TABLE]

Appendix B Dissipative out-of-time-order correlators

B.1 Short-time expansion of OTOC

It is insightful to consider the short-time expansion of the OTOC

[TABLE]

The linear term in $t$ has an oscillatory contribution $-4i\gamma t\mathrm{Tr}([\hat{L},\hat{A}][\hat{L},[\hat{H}_{0},\hat{A}]])$ and a dissipative contribution

[TABLE]

where we used the cyclic property of the trace to write $\mathrm{Tr}(\hat{X}[\hat{Y},\hat{Z}])=\mathrm{Tr}([\hat{X},\hat{Y}]\hat{Z})$ with $\hat{X}=[\hat{L},\hat{A}],\;\hat{Y}=\hat{L},$ and $\hat{Z}=[\hat{L},[\hat{L},\hat{A}]]$. Then the inverse dissipation time, $\tau_{D}^{-1}\propto 2\gamma\mathrm{Tr}(\mathcal{D}[\hat{A}]^{2})=4(\mathcal{D}[\hat{A}],\mathcal{D}[\hat{A}])\geq 0$ is related to the Hilbert-Schmidt norm of the dissipator $\mathcal{D}[\bullet]=-[\hat{L},[\hat{L},\bullet]]$ acting on the initial operator $\hat{A}$. The short-time behavior is then determined by the exponential

[TABLE]

This approximation is compared to the full OTOC in Fig. 4, where the crossover between the different exponential regimes of the dissipative OTOC is apparent.

B.2 The dissipative OTOC for $[\hat{H}_{0},\hat{L}]=0$

In the simple case that $[\hat{H}_{0},\hat{L}]=0$, the operators $\hat{H}_{0}$ and $\hat{L}$ share a common eigenbasis, such that $\hat{H}_{0}=\sum_{n}E_{n}\ket{n}\bra{n}$ and $\hat{L}=\sum_{n}l_{n}\ket{n}\bra{n}$. The solution of the adjoint master equation then reads

[TABLE]

The behavior of the dissipative OTOC then is simply

[TABLE]

This shows that the short-time dynamics is governed by the two eigenvalues with largest difference $\max_{m,k}(l_{m}-l_{k})^{2}$ over which the operator has support, i.e. $A_{mk}\neq 0$, and the long-time dynamics is governed by the smallest possible non-zero difference $\min_{m,k}(l_{m}-l_{k})\neq 0$ over which the initial operator has support.

Appendix C Product of operators as a doubled Hilbert space operation

The standard product between operators in a Hilbert space $\mathscr{H}$ with closure relation $\hat{\mathbb{1}}_{\mathscr{H}}=\sum_{n}\ket{n}\bra{n}$ reads

[TABLE]

In turn, the tensor product of the operators over a doubled Hilbert space $\mathscr{H}=\mathscr{H}_{1}\otimes\mathscr{H}_{2}$ with closure relation $\mathbb{1}=\mathbb{1}_{\mathscr{H}_{1}}\otimes\mathbb{1}_{\mathscr{H}_{2}}=\sum_{n_{1},n_{2}}\ket{n_{1},n_{2}}\bra{n_{1},n_{2}}$ where $\ket{n_{1},n_{2}}=\ket{n_{1}}\otimes\ket{n_{2}}$ reads

[TABLE]

Thus, applying the swap operator, defined from $\mathbb{S}\ket{n_{1}}\otimes\ket{n_{2}}=\ket{n_{2}}\otimes\ket{n_{1}}$, yields

[TABLE]

The partial trace over the second Hilbert space $\mathrm{Tr}_{\mathscr{H}_{2}}(\bullet)=\sum_{n_{2}}\braket{n_{2}}{\bullet}{n_{2}}$ gives

[TABLE]

Appendix D Quantum stochastic Lipkin Meshkov Glick model

D.1 Transport in quantum sLMG

In the main text, we presented the eigenvalues of the SOV and their transport analog in the sLMG for $[\hat{H}_{0},\hat{L}]=0$, see Fig. 2 (a). The case in which $[\hat{H}_{0},\hat{L}]\neq 0$ shows a richer phenomenology, even for the simple sLMG model. The eigenvalues of the SOV are shown in Fig. 5. The smallest eigenvalue first evolves ballistically, $\Lambda_{0}(t)\sim t^{2}$, then turns into superdiffusive $\Lambda_{0}(t)\sim t^{3/2}$ before changing further to a slower growth $t^{\alpha}$ with $\alpha<1$, to ballistic again $\Lambda_{0}(t)\sim t^{2}$ and then saturates. Another interesting feature is that at long times all eigenvalues saturate to the same value, a feature not present in the case $[\hat{H}_{0},\hat{L}]=0$.

D.2 Visualizing the quantum stochastic operator variance

The SOV $\Delta\hat{A}_{t}^{2}$ that we introduced is an operator over the Hilbert space. We illustrate this quantity for the case of the quantum sLMG model. In general, this operator is $d$ dimensional, and it has non-zero components over all the elements of the considered operator basis. However, since the sLMG is a mean-field description, we are mainly interested in its projection over the subspace spanned by $\{\mathbb{1},\hat{S}_{x},\hat{S}_{y},\hat{S}_{z}\}$. We would like to have these elements be orthonormal with each other. For this reason, we choose the Hilbert-Schmidt inner product

[TABLE]

where the normalization has been chosen such that $(\hat{S}_{i},\hat{S}_{j})=\delta_{ij}$ and comes from $\mathrm{Tr}(\hat{S}_{z}^{2})=\sum_{j=0}^{2S}(S-j)^{2}=S(S+1)(2S+1)/3.$ Note that this normalization leads to $(\mathbb{1},\mathbb{1})=\frac{3}{S(S+1)}\neq 1$.

We further introduce a notion of Stochastic Operator Standard Deviation (SOSD) as the matrix square root of the SOV, $\Delta\hat{A}_{t}=\sqrt{\Delta\hat{A}_{t}^{2}}$. The SOV and its deviation are thus illustrated in Fig. 6.

Appendix E Classical stochastic LMG model

E.1 Classical limit of the LMG Hamiltonian

Let us introduce the SU(2) coherent states as [62]

[TABLE]

where $\hat{S}_{+}=\hat{S}_{x}+i\hat{S}_{y}$ is the raising operator and $\ket{S,-S}$ is the eigenstate of $\hat{S}^{2}$ and $\hat{S}_{z}$ with smallest $z$ component of the spin, namely $\hat{S}^{2}\ket{S,-S}=S(S+1)\ket{S,-S}$ and $\hat{S}_{z}\ket{S,-S}=-S\ket{S,-S}$. The coherent states $\ket{\zeta}$ represent a wavepacket with the minimum width allowed by the uncertainty principle localized around $\mathbf{n}=(\sin\theta\cos\phi,\sin\theta\sin\phi,\cos\theta)$ where $\zeta=-\tan\frac{\theta}{2}e^{-i\phi}$ and corresponds to a stereographic projection. Following [38], we define the classical Hamiltonian $H_{\textsc{lmg}}$ as

[TABLE]

where $S=N/2$. The expectation value of the relevant spin operators between coherent SU(2) states is [62]

[TABLE]

The classical LMG Hamiltonian then reads

[TABLE]

We then introduce the canonical variables $Q,\;P$ as

[TABLE]

that yield the relations

[TABLE]

Inverting theses equations gives

[TABLE]

that can be substituted in (39) and yield

[TABLE]

This corresponds to the classical Hamiltonian given in the main text, up to a constant shift in energy.

E.2 Heuristic reason for stabilization of sLMG

The stabilization we discussed in the main text can be understood in terms of the energy landscape associated with the stochastic Hamiltonian, illustrated in Fig. 7. In the limit that $\sqrt{2\gamma}\ll 1$, most realizations are as in the noiseless LMG. The behavior is then close to that of the deterministic model, with the double well (DW) phase showing an unstable point and the single well (SW) phase being stable. In the opposite limit, $\sqrt{2\gamma}\gg 1$, almost half of the realizations flip the Hamiltonian since $(1+\sqrt{2\gamma}\xi_{t})<0$, as shown in Fig. 7 (a). This brings a major difference between the DW and SW phases. In the SW phase, all points in phase space $(Q,P)$ either gain (d)$\rightarrow$(e) or lose (e)$\rightarrow$(d) energy. In the DW phase, however, when going from (b)$\rightarrow$(c), the points inside the wells lose energy while the points outside the wells gain it, and viceversa for (c)$\rightarrow$(b). This difference between the points of phase space provides a rationale for the stochastic stabilization seen for the DW phase—blue region in Fig. 3(b) of the main text.

E.3 Analytical determination of the Lyapunov exponent for the sLMG model

We consider a vector variable $\mathbf{u}_{t}$ subject to deterministic and fluctuating external perturbations, $\mathbb{A}_{d}$ and $\mathbb{A}_{s}(t)$, respectively. Its equation of motion is described by the stochastic differential equations

[TABLE]

Following van Kampen [63], we go in the interaction picture with respect to the deterministic evolution and consider $\mathbf{v}_{t}=e^{-\mathbb{A}_{d}t}\mathbf{u}_{t}$. In second-order of $\sqrt{\gamma}$, the average with fixed initial conditions evolves as

[TABLE]

valid for $\sqrt{2\gamma}t\ll 1$. We recognize the solution to order $\gamma$ of the linear differential equation, also known as Bourett’s integral equation, written back in the original representation as

[TABLE]

This equation is derived assuming the standard rules of calculus, and thus assumes Stratonovich formalism. The latter defines the stochastic integral as ${\int_{0}^{t}\delta(t-\tau)f(\tau)d\tau=\frac{1}{2}f(t)}$. So in the case that $\mathbb{A}_{d}$ and $\mathbb{A}_{s}$ commute, and for $\mathbb{A}_{s}(t)=\xi_{t}\mathbb{A}_{s}$ fluctuating with Gaussian white noise, Eq. (46) simplifies to

[TABLE]

where the change of integration variable $\tau^{\prime}=t-\tau$ brings the minus sign. In systems exhibiting chaos, the Lyapunov gives the exponential divergence of the trajectory. We interpret this as the maximum eigenvalue of $\mathbb{A}_{d}-\gamma\mathbb{A}_{d}^{2}$.

The LMG at the origin, $Q=P=0$, can be linearized into the harmonic oscillator $H=\frac{1}{2}[\Omega P^{2}+(\Omega-2)Q^{2}]$. Hamilton’s equation of motion gives $\dot{Q}_{t}$ and $\dot{P}_{t}$, from which we can compute the evolution of the quadratic terms as [64]

[TABLE]

The maximum eigenvalue of this matrix is $\lambda_{\textsc{lmg}}=\sqrt{2\Omega-\Omega^{2}},$ which recovers the Lyapunov exponent at the origin in the noiseless LMG model [38]. For the sLMG, we add noise in the energy scale and consider an evolution dictated by

[TABLE]

where $\xi_{t}$ is Gaussian white noise. The maximum eigenvalue of $\mathbb{A}_{d}-\gamma\mathbb{A}_{d}^{2}$ thus gives the average Lyapunov exponent as

[TABLE]

Appendix F Details on the numerical solutions

F.1 Vectorization

The formal solution to the master equation reads $\langle\hat{A}_{t}\rangle=e^{\mathcal{L}^{\dagger}t}[\hat{A}]$, where $\mathcal{L}^{\dagger}[\bullet]$ is the adjoint Liouvillian superoperator. To numerically solve it, we apply vectorization by which operators $\hat{A}_{t}$ become vectors in a larger Hilbert space $|\hat{A}_{t})$. Superoperators become operators over this larger Hilbert space through the mapping

[TABLE]

therefore the vectorized Lindbladian reads

[TABLE]

F.2 Solver of Stochastic Differential Equations

We solve the SDE using the explicit order 1.0 strong scheme [65], which briefly consists of the following. Consider a Itō SDE of the form

[TABLE]

where $X_{0}$ is the initial condition. Let $Y_{n}$ be the solution at time $t_{n}=n\delta,$ where $\delta$ is the time-step. We first set the initial condition $Y_{0}=X_{0}$ and compute recursively the solution as

[TABLE]

where $\Upsilon_{n}=Y_{n}+a_{n}\delta+b_{n}\sqrt{\delta}$, $a_{n}=a(Y_{n}),\;b_{n}=b(Y_{n})$, and $\Delta W_{n}=W_{n+1}-W_{n}$ are independent identically distributed normal random variables with zero average and variance $\delta$.