1 Introduction↩︎

In recent years, breakthroughs in quantum technology have ushered in a new era of experimental quantum physics. From ultracold atomic gases and trapped ions to the development of programmable quantum computers, these platforms provide unprecedented opportunities to probe the intricate dynamics of quantum many-body systems far from equilibrium [1], [2]. One of the most compelling frontiers in this field is understanding how information behaves in these complex systems. While traditional approaches often focus on entanglement and thermalization [3]–[7], a more refined concept known as quantum information scrambling has emerged as a crucial lens for characterizing non-equilibrium dynamics [8]–[17].

Quantum information scrambling describes the process by which local information becomes irreversibly encoded into non-local degrees of freedom throughout a system. It is a universal feature of chaotic, strongly interacting quantum systems. To grasp this concept, one can imagine dropping a single droplet of colored ink into a glass of water. Initially, the information about the ink droplet’s color and position is localized. As the ink disperses due to diffusion, the information becomes more difficult to recover from any local observation [18], [19]. In a quantum system, unitary time evolution causes a similar effect: an initial local perturbation becomes distributed across the entire system, rendering it unrecoverable by any local measurement. This is a more subtle process than thermalization, which simply describes how a local subsystem reaches a state of thermal equilibrium with its environment: scrambling however concerns where the lost information goes, propagating it to the system’s global, non-local degrees of freedom.

The out-of-time-order correlator (OTOC) [20] has been identified as a uniquely powerful probe for this phenomenon in diverse fields, ranging from high-energy physics to condensed matter physics [12], [15]–[17]. Analogous to the Lyapunov exponent in classical theory of chaos, which measures the sensitivity to initial conditions, the OTOC physically quantifies how a local perturbation spreads or scrambles, throughout a many-body quantum system and mathematically quantifies the growth of a commutator between two initially local operators \(W\) and \(V\). At time \(t=0\), the commutator \([W(0),V(0)]\) is zero (or small if the operators are spatially separated), as the operators have no shared degrees of freedom. As time evolves \(t=\tau\), however, the time-evolved operator \(W(\tau)\) spreads across the system due to interactions. When its support reaches and overlaps with \(V(0)\), the commutator \([W(\tau),V(0)]\) becomes non-zero and its magnitude grows. This operator spreading is a direct signature of quantum information scrambling. The rate of this growth provides a clear distinction between different dynamical regimes: chaotic systems typically exhibit an early-time exponential growth, while integrable systems show a much slower, polynomial growth. This observation makes the OTOC a powerful diagnostic tool for quantum chaos.

Despite its profound theoretical significance, the experimental measurement of the OTOC presents a formidable challenge [17]. The nature behind the OTOC—out-of-time-order of operators—implies a sequence of forward and backward time evolutions, which is notoriously difficult to implement with high fidelity, especially in large, complex, strongly interacting systems [21], [22]. Consequently, much of the research in this area has shifted from theoretical characterization to the development of ingenious experimental protocols that can indirectly access the OTOC, e.g., rewinding time [23], [24], interferometer [25], [26], weak-measurement [27]–[29], two-point measurement [30], Rényi entropy [31], multiple quantum coherence intensity [32], quantum teleportation [33], randomized measurement [34], thermofield double state [35], classical shadows [36], irreversibility [37] and so forth.

In this paper, we focus on the evaluation of the OTOC using the trapped-ion quantum computer reimei. We specifically investigate a method based on the principle of irreversibility proposed in the companion Letter [37], which re-frames the OTOC as the difficulty of recovering an initial state after a specific sequence of operations. Our main contribution is the first experimental demonstration of this irreversibility-susceptibility method for a non-trivial many-body system, the XXZ spin chain, prepared in a thermal state at finite temperature. To provide a comprehensive understanding of its performance, we also conduct a comparative analysis of this method against two other prominent techniques: the rewinding time method and the weak-measurement method. By implementing these distinct methods and analyzing their performance, our work offers crucial insights into their practical applicability on a near-term quantum computer. This research not only validates quantum computers as powerful platforms for exploring fundamental concepts in quantum chaos but also provides a practical framework for future investigations.

This paper is structured as follows. In Sec. 2, we provide a detailed mathematical definition of the OTOC and discuss its connection to information scrambling. In Sec. 3, we review three key methods for its evaluation: the rewinding time, weak-measurement, and irreversibility-susceptibility methods. In Sec. 4, we describe our experimental setup on a quantum computer, including the state preparation algorithm and the Hamiltonian, and present our experimental results. Finally, in Sec. 5, we conclude with a summary of our findings and a discussion of future prospects.

2 Theoretical framework of the OTOC↩︎

The OTOC provides a powerful theoretical framework for quantifying information scrambling and diagnosing quantum chaos in many-body systems [12], [20]. Its core concept revolves around the dynamics of two initially local operators.

2.1 Operator spreading and OTOC↩︎

Consider two initially local operators, \(W\) and \(V\). In the Heisenberg picture, an operator \(W\) evolves in time as \(W(\tau):= U^{\dagger}(\tau)W(0)U(\tau)\), where \(U(\tau) = \exp(-\ifmmode\mathrm{i}\else\i \fi H\tau)\) is the time-evolution operator generated by the Hamiltonian \(H\) of the system \(\bfS\) under consideration. While \(W(0)\) and \(V(0)\) might be simple operators (e.g., a Pauli operator on a single qubit), \(W(\tau)\) generally becomes a complex, non-local operator that acts on a larger portion of \(\bfS\) as time progresses. This phenomenon, known as operator spreading, is a key signature of information scrambling.

The OTOC is defined by the thermal expectation value of the squared commutator between \(W(\tau)\) and \(V(0)\): \[\begin{align} C_{\beta}(\tau):=-\langle [W(\tau),V(0)]^{2} \rangle_{\beta}.\label{eq:def95otoc} \end{align}\tag{1}\] Here, we abbreviate the Gibbs state \(\rho(\beta,H)=\exp(-\beta H)/Z_{\beta}\) and the partition function \(Z(\beta,H)=\mathop{\mathrm{Tr}}[\exp(-\beta H)]\) to \(\rho_{\beta}\) and \(Z_{\beta}\), respectively; \(\langle \bullet \rangle_{\beta}=\mathop{\mathrm{Tr}}(\rho_{\beta}\bullet)\) denotes the thermal average with respect to \(\rho_{\beta}\). This definition is sensitive to the commutator, which is small if the operators are spatially separated, but grows as the time-evolved operator \(W(\tau)\) spreads and overlaps with \(V(0)\).

A related quantity, the four-point correlator \(F_{\beta}(\tau)\), is often studied for its computational convenience: \[\begin{align} F_{\beta}(\tau):=\langle W^{\dagger}(\tau)V^{\dagger}(0)W(\tau)V(0) \rangle_{\beta}.\label{eq:def95correlator} \end{align}\tag{2}\] If \(W\) and \(V\) are unitary operators, the squared commutator \(C_{\beta}(\tau)\) can be expressed in terms of the real part of \(F_{\beta}(\tau)\) as \[\begin{align} C_{\beta}(\tau)=2\{1-\relax [F_{\beta}(\tau)]\}.\label{eq:otoc-correlator} \end{align}\tag{3}\] If, in addition, \(W\) and \(V\) are also Hermitian, e.g. Pauli operators, this relationship simplifies to \(C_{\beta}(\tau)=2[1-F_{\beta}(\tau)]\), where \(F_{\beta}(\tau)\) becomes a real quantity. The squared commutator can also be expressed in terms of the Frobenius norm of the commutator, \(C_{\beta}(\tau)=\|\sqrt{\rho_{\beta}}[W(\tau), V(0)]\|^{2}_{2}\).

2.2 OTOC as a probe of scrambling↩︎

The significance of the OTOC lies in its ability to diagnose operator spreading and chaos. In a quantum chaotic system, the operator \(W(\tau)\) grows in complexity and size. At \(t=0\), the commutator \([W(0),V(0)]\) is zero (or small if the operators are spatially separated). As time evolves \(t=\tau\), the support of \(W(\tau)\) expands. When this expanding operator front overlaps with the support of \(V(0)\), the commutator becomes non-zero, and its magnitude grows [38]. This growth signifies the scrambling of quantum information: an initial local perturbation by \(V\) becomes entangled with the entire system, and its effect can only be detected by a later, non-local operation \(W\) [39], [40].

In systems with a classical chaotic analog, the OTOC often exhibits an early-time exponential growth captured by \[\begin{align} C_{\beta}(\tau)\sim\mathrm{e}^{2\lambda_{L}\tau}, \end{align}\] where \(\lambda_{L}\) is the quantum Lyapunov exponent [12], [41]–[44]. This exponential growth marks the initial phase of scrambling. At long times, the OTOC saturates at a constant value due to the finite size of the system, indicating that the operator has spread across the entire system.

2.3 Operator expansion and finite temperature↩︎

A more microscopic understanding of the OTOC can be gained by expanding the Heisenberg operator \(W(\tau)\) in a complete basis of operators \(\{\Gamma_{i}\}\), such as the Pauli strings for a system of qubits: \[\begin{align} W(\tau)=\sum_{i}\gamma_{i}(\tau)\Gamma_{i}. \end{align}\] Here, \(|\gamma_{i}(\tau)|^{2}\) can be interpreted as the probability of finding the operator \(W(\tau)\) in the state \(\Gamma_{i}\). The thermal average of the OTOC is directly related to these probabilities. At infinite temperature, the expectation value \(\langle \bullet \rangle_{\beta}\) is replaced by a simple trace, \(\mathop{\mathrm{Tr}}(\bullet)/\mathop{\mathrm{Tr}}(\openone)\). In this case, it can be shown that the average OTOC over all local operators \(\Gamma_{i,r}\) at a site \(r\) is proportional to the probability that the operator \(W(\tau)\) has expanded to include that site as a non-trivial component. The average squared commutator is then a complementary measure: \[\begin{align} &\frac{1}{d^{2}-1}\sum_{i:\Gamma_{i,r}\neq\openone_{r}}\frac{1}{\mathop{\mathrm{Tr}}(\openone)}\mathop{\mathrm{Tr}}([W(\tau),\Gamma_{i,r}]^{2})\nonumber\\ &\quad=\frac{2d^{2}}{d^{2}-1}\sum_{i:\Gamma_{i,r}\neq\openone_{r}}|\gamma_{i}(\tau)|^{2}, \end{align}\] where \(d\) is the dimension of the local Hilbert space. This expression elegantly links the OTOC to the size of the Heisenberg operator, which is defined by the number of sites with a non-trivial Pauli operator in its expansion [17].

While infinite-temperature OTOCs are useful for studying the fundamental physics of scrambling, it is also crucial to consider the finite-temperature case, where the expectation value is taken with respect to a thermal Gibbs state \(\rho_{\beta}\). The physics in this regime is more complex due to correlations present in the thermal state. For instance, there are multiple regularized versions of the measure for the OTOC, e.g., \[\begin{align} \tilde{F}_{\beta}(\tau):=\mathop{\mathrm{Tr}}[\rho^{\kappa_{1}}W(\tau)\rho^{\kappa_{2}}V(0)\rho^{\kappa_{3}}W(\tau)\rho^{\kappa_{4}}V(0)], \end{align}\] where \(\sum_{i}\kappa_{i}=1\). These different regularizations correspond to displacing the operators in imaginary time and can lead to different physical results, such as a dependence of the butterfly velocity on the choice of the regularization [13], [39], [40], [45]–[52]. Nevertheless, the core principle of the OTOC as a probe of operator spreading remains intact, providing a powerful tool for diagnosing information scrambling across various physical regimes.

3 Evaluation methods↩︎

While the experimental evaluation of the OTOC is a non-trivial task, each method presents unique theoretically interesting advantage and experimental challenges. We now explain these experimental evaluation methods of the OTOC, specifically focussing on the rewinding time method (RTM) [23], [24], weak-measurement method (WMM) [27]–[29], and irreversibility-susceptibility method (ISM) [37] for the comparison in our numerical experiments.

3.1 Rewinding time method↩︎

The RTM [23], [24] provides a direct way to measure the OTOC by relating it to the outcome of a quantum interference experiment [21], [22]. The central idea of this method is to measure the overlap between two quantum states, which are prepared by applying local operators and time evolution in different orders.

Let the target system \(\bfS\) be initialized in a pure state \(\ket{\psi}\) for simplicity. We are interested in the overlap of the states \[\begin{align} \ket{\psi_{1}}&=V(0)W(\tau)\ket{\psi},\\ \ket{\psi_{2}}&=W(\tau)V(0)\ket{\psi}. \end{align}\] The OTOC is related to this overlap, and an interferometric setup allows us to precisely quantify the relationship.

To measure \(F_{\beta}(\tau)\), we employ a practical quantum circuit that uses an ancillary qubit \(\bfC\) as a control system.

1. Prepare \(\bfC\) in the eigenstate \(\ket{+}_{\bfC} = (\ket{0}_{\bfC}+\ket{1}_{\bfC})/\sqrt{2}\) of the Pauli-\(x\) operator \(\sigma^{x}\), which corresponds to the eigenvalue \(+1\), and \(\bfS\) in the desired state \(\ket{\psi}_{\bfS}\). The initial state of composite system \(\bfC+\bfS\) is \(\ket{\Psi_{\text{in}}}=\ket{+}_{\bfC}\otimes\ket{\psi}_{\bfS}\).

2. Apply the controlled-operator \(| {0} \rangle\langle {0} |_{\bfC}\otimes\openone_{\bfS}+| {1} \rangle\langle {1} |_{\bfC}\otimes V(0)_{\bfS}\). \(V(0)\) is a local operator, meaning it acts only on a specific part of \(\bfS\).

3. Apply the forward time evolution \(\openone_{\bfC}\otimes U(\tau)_{\bfS}\).

4. Apply the operator \(\openone_{\bfC}\otimes W_{\bfS}\).

5. Apply the backward time evolution \(\openone_{\bfC}\otimes U^{\dagger}(\tau)_{\bfS}\).

6. Apply the controlled-operator \(| {0} \rangle\langle {0} |_{\bfC}\otimes V(0)_{\bfS}+| {1} \rangle\langle {1} |_{\bfC}\otimes\openone_{\bfS}\).

After these operations, the state of \(\bfC+\bfS\) becomes an entangled state \[\ket{\Psi_{\text{out}}}=\frac{1}{\sqrt{2}}\left(\ket{0}_{\bfC}\ket{\psi_{1}}_{\bfS}+\ket{1}_{\bfC}\ket{\psi_{2}}_{\bfS} \right).\] The expression for this output state demonstrates the core principle of the interferometer: the two paths of the interferometer (the \(\ket{0}_{\bfC}\) and \(\ket{1}_{\bfC}\) branches) contain different evolutions of \(\bfS\). Its inclusion in the circuit is crucial for creating the out-of-time-ordered structure of operators.

In order to extract \(F_{\beta}(\tau)\) from output state \(\ket{\Psi_{\text{out}}}\), we perform a measurement on \(\bfC\). By measuring the expectation values of the Pauli operators \(\sigma^{x}\) and \(\sigma^{y}\) on \(\bfC\), we can directly determine the real and imaginary parts of \(F_{\beta}(\tau)\), i.e., \[\begin{align} F_{\beta}(\tau)=\langle \sigma^{x} \rangle_{\bfC}+\ifmmode\mathrm{i}\else\i \fi\langle \sigma^{y} \rangle_{\bfC}, \end{align}\] where \[\begin{align} \langle \sigma^{x} \rangle_{\bfC}&=\relax \left[\bra{\psi}V^{\dagger}(0)W(\tau)V(0)W^{\dagger}(\tau)\ket{\psi}\right]\nonumber\\ &=\relax F_{\beta}(\tau),\\ \langle \sigma^{y} \rangle_{\bfC}&=\relax \left[\bra{\psi}V^{\dagger}(0)W(\tau)V(0)W^{\dagger}(\tau)\ket{\psi}\right]\nonumber\\ &=\relax F_{\beta}(\tau). \end{align}\] Therefore, the OTOC can be calculated from Eq. 3 using \(F_{\beta}(\tau)\). The RTM was used in Ref. [53] for measuring the time-dependent evolution and fluctuation of the OTOC to experimentally investigate the dynamics of quantum scrambling. For the mixed state case as the input state, the specific generalization was proposed using the thermal field double state at finite temperature [54]; and it was employed in Ref. [55]

3.2 Weak-measurement method↩︎

The WMM offers an alternative approach for measuring the OTOC that elegantly circumvents the need for backward time evolution [27]–[29]. This approach relates the OTOC to a sequence of weak measurements, where the information is extracted from the target system with minimal disturbance. The method is particularly powerful when the operators of interest, \(W\) and \(V\), are unitary and satisfy the condition \(W^{2} = \openone\) and \(V^{2} = \openone\), which is called the dichotomic operator, such as Pauli operators, since it can be implemented with generalized measurements of any strength, not just weak ones. This property eliminates the issue of disturbance typically associated with weak measurements [56], [57].

The core idea is to measure the expectation values of nested commutators and anticommutators of the operators at different times. These values can be extracted by performing a sequence of weak measurements and averaging the outcomes with specific weights. To access the complex nature of the OTOC, we need a canonical type of measurement \[\begin{align} M^{A}_{\pm}(\phi)&:=\bra{z\pm}S_{A}(\phi)\ket{x+}\nonumber\\ &=\frac{1}{\sqrt{2}}\left[\cos\left(\frac{\phi}{2}\right)\openone\pm\sin\left(\frac{\phi}{2}\right)A\right], \end{align}\] where \[\begin{align} S_{A}(\phi)=\exp\left(-\ifmmode\mathrm{i}\else\i \fi\frac{\phi}{2}A\otimes Y\right) \end{align}\] for an arbitrary observable \(A\) and we use the notations \(\ket{\bullet\pm}\) for \(\bullet=x,y,z\) so that we can distinguish the bases of Pauli operators. \(M^{A}_{\pm}(\phi)\) causes a partial collapse of the state onto the eigenbasis of an operator \(A\). The outcome of the measurement \(M^{A}_{\pm}(\phi)\) is weighted by modified eigenvalues that correspond to the anticommutator \(\{A,B\}:=AB+BA\). The real part of the correlator is extracted by performing a sequence of the measurements and averaging the results with specific weights.

For clarity, we now replace the \(\pm\) notation with explicit labels, e.g., \(\pm1\to(-1)^{a}\) with \(a\in\{0,1\}\), which indicate the experimental outcome obtained when measuring the ancilla basis. The strength of the measurement \(M^{A}_{a}(\phi)\) is controlled by a coupling-strength \(\phi\in(0,\pi/2]\) that ranges from a near-identity transformation to a projective measurement, allowing the tuning of the backaction to the system from weak \(\phi=0\) to strong \(\phi=\pi/2\).

We define the rescaled value that the experimenter should assign each observed ancilla outcome \(a\in\{0,1\}\), \[\begin{align} \alpha_{a}(\phi)=\frac{(-1)^{a}}{\sin(\phi/2)}. \end{align}\] The values \(\alpha_{a}(\phi)\) act as modified eigenvalues of the self-adjoint operator \(A\), i.e., \(A\) can be decomposed into the POVM for the measurement implemented by \(M^{A}_{\pm}(\phi)\), \[\begin{align} \sum_{a=0,1}\alpha_{a}(\phi)[M^{A}_{a}(\phi)]^{\dagger}M^{A}_{a}(\phi)=A. \end{align}\] As a particularly important special case, when \(\phi=\pi/2\), the values \(\alpha_{a}(\phi)=(-1)^{a}\) reduce to the eigenvalues of \(A\) and the measurements are projective ones. These measurements can be implemented by coupling the system to an ancilla qubit with a controlled rotation gate and then measuring the ancilla. The strength of the measurement is controlled by a coupling angle \(\theta \in (0, \pi/2]\). Weak measurements correspond to a small angle \(\theta \to 0\).

Typically, determining complex quantities like operator correlators requires the use of weak measurements (\(\phi\approx0\)) to prevent disturbance. In special cases, however, relevant information may still be contained in the collected measurement statistics in spite of any disturbance. This is the case for qubits, where the following remarkable identities hold for any coupling-strength and thus enable the improved correlator measurement protocols; the anticommutator identities \[\begin{align} \sum_{a=0,1}\alpha_{a}(\phi)M^{A}_{a}(\phi)\rho[M^{A}_{a}(\phi)]^{\dagger}&=\frac{1}{2}\{A,\rho\},\\ \sum_{a=0,1}\alpha_{a}(\phi)[M^{A}_{a}(\phi)]^{\dagger}BM^{A}_{a}(\phi)&=\frac{1}{2}\{B,A\}. \end{align}\] To measure \(F_{\beta}(\tau)\), we consider the specific case where the operators \(W\) and \(V\) are dichotomic operators. The protocol to obtain the real part of \(F_{\beta}(\tau)\) is as follows.

1. Prepare each probe system \(\bfP_{\bullet}\) with \(\bullet=v,w,v',w'\) in the eigenstate \(\ket{+}\) of the Pauli-\(x\) operator \(\sigma^{x}\) and \(\bfS\) in the desired state \(\rho_{\bfS}\).

2. Perform a weak measurement \(M^{V}_{v}(\phi_{v})\) on the target system \(\bfS\).

3. Evolve \(\bfS\) forward in time by \(U(\tau)\) for a duration from \(t=0\) to \(t=\tau\).

4. Perform a weak measurement \(M^{W}_{w}(\phi_{w})\) on \(\bfS\).

5. Evolve \(\bfS\) backward in time by \(U(\tau)\) for a duration from \(t=0\) to \(t=\tau\).

6. Perform a weak measurement \(M^{V}_{v'}(\phi_{v'})\) on \(\bfS\).

7. Evolve \(\bfS\) forward in time by \(U(\tau)\) for a duration from \(t=0\) to \(t=\tau\).

8. Perform a weak measurement \(M^{W}_{w'}(\phi_{w'})\) on \(\bfS\).

Then, we average the products of the modified eigenvalues \(\alpha_{v}(\phi_{v}),\alpha_{w}(\phi_{w}),\alpha_{v'}(\phi_{v'}),\alpha_{w'}(\phi_{w'})\) from each measurement. Remarkably, the result of the averaging procedure directly yields a quantity related to the real part of \(F(t)\) being of the form \[\begin{align} &\sum_{\substack{\alpha_{v}(\phi_{v}),\alpha_{w}(\phi_{w}),\\\alpha_{v'}(\phi_{v'}),\alpha_{w'}(\phi_{w'})\in\{0,1\}}}\alpha_{v}(\phi_{v})\alpha_{w}(\phi_{w})\alpha_{v'}(\phi_{v'})\alpha_{w'}(\phi_{w'})\nonumber\\ &\quad\times\Pr\{\alpha_{v}(\phi_{v}),\alpha_{w}(\phi_{w}),\alpha_{v'}(\phi_{v'}),\alpha_{w'}(\phi_{w'})\|\rho\}\nonumber\\ &\quad=\frac{1}{2^{3}}\langle \{\{\{W(\tau),V(0)\},W(\tau)\},V(0)\} \rangle_{\rho}\nonumber\\ &\quad=\frac{1}{2}(1+\relax [F_{\beta}(\tau)])\nonumber\\ &\quad=1-\frac{1}{4}C_{\beta}(\tau). \end{align}\] Note that this identity holds for any strength of the measurements. Crucially, this method works even with strong, projective measurements, but the use of weak measurements allows for a reduced backaction on \(\bfS\), which can be useful in certain experimental contexts.

3.3 Irreversibility-susceptibility method↩︎

The ISM leverages a recently developed formulation that connects the OTOC to the fundamental notion of irreversibility. The core idea is that the OTOC, which quantifies quantum scrambling, can be reformulated as a specific measure of irreversibility of a quantum process [see the companion Letter [37] for more detail of the framework]. This formulation provides a thermodynamically meaningful perspective on the OTOC and offers an experimentally tractable way to evaluate it.

First, let us define the measure of irreversibility proposed in Ref. [58], which is a useful tool to investigate universal limitations [58]–[61] on the dynamics in various fields including black holes [8], [33], [58], [62]–[67], quantum thermodynamics [58], [59], [68], [69], quantum measurements [58], [60], [70]–[77], quantum computing [58], [78]–[81], and channel costs of general resource theories [61]. We consider a quantum process represented by a completely positive trace-preserving (CPTP) map \(\calL\) that takes a state from a target system \(\bfS\) to another system \(\bfS'\). For a given ensemble of input states \(\Omega=\{p_{k},\rho_{k}\}\), where \(\{\rho_{k}\}\) is a family of quantum states prepared with a probability distribution \(\{p_{k}\}\), the irreversibility of the process \(\calL\) is defined by \[\begin{align} \delta(\calL,\Omega)&:=\min_{\calR:\bfS'\rightarrow\bfS}\sqrt{\sum_{k}p_{k}\delta^{2}_{k}},\label{eq:irreversibility} \end{align}\tag{4}\] as the minimum recovery error over all possible recovery maps \(\calR\). Here, the minimization is performed over all CPTP maps \(\calR\), which attempt to reverse the process \(\calL\). \(\delta_k\) quantifies the distance between the original state \(\rho_k\) and the state after the loss and recovery processes \(\calR\circ\calL(\rho_{k})\): \[\begin{align} \delta_{k}&:=D_{F}(\rho_{k},\calR\circ\calL(\rho_{k})),\label{eq:distance} \end{align}\tag{5}\] where \(D_{F}(\rho,\sigma):=\sqrt{1-F(\rho,\sigma)^{2}}\) is the purified distance [82], which is based on the Uhlmann fidelity \(F(\rho,\sigma):=\mathop{\mathrm{Tr}}\left(\sqrt{\sqrt{\rho}\sigma\sqrt{\rho}}\right)\) [83]. This definition of irreversibility has the benefit, providing operational meanings of a variety of irreversibility measures: it gives a minimum threshold for the entropy production [84] in stochastic thermodynamics, a lower bound for the entanglement fidelity error [85] in quantum error correction, and a recovery error of the Petz recovery map [86]–[88].

The OTOC is defined by the following expression, which directly relates it to a specific instance of irreversibility: \[\begin{align} C_{\beta}(\tau)&=\lim_{\theta\rightarrow 0}\frac{\delta^{2}(\mathcal{L}_{\rho,V,\theta,\calD_{W}},\Omega_{1/2,\pm},\calR_{V,\bfS})}{\theta^{2}},\label{eq:specific95dist61otoc} \end{align}\tag{6}\] Next, we break down the components of this formula. The process \(\calL_{\rho,V,\theta,\calD_{W}}\) is constructed from a composite map that simulates the time evolution under the dynamics of the OTOC: \[\begin{align} \calL_{\rho,V,\theta,\calD_{W}}:=(\calD_{W}\otimes\openone_{\mathbf{Q}})\circ\mathcal{U}_{V,\theta}\circ\mathcal{A}_{\rho}, \end{align}\] which involves

1. Apply an appending process \(\calA_{\rho}(\bullet):=\bullet\otimes\rho\) that prepares the initial state \(\rho\) of \(\bfS\) on an ancilla qubit system \(\mathbf{Q}\), where the initial state of \(\bfQ\) is prepared from a specific test ensemble \(\Omega_{1/2,\pm}\).

2. Apply a weak interaction \(\mathcal{U}_{V,\theta}(\bullet):=U_{V}(\theta)(\bullet)U^{\dagger}_{V}(\theta)\) on \(\bfS+\bfQ\), where \(U_{V}(\theta)=\exp\left[-\ifmmode\mathrm{i}\else\i \fi\theta (Z\otimes V)\right]\).

3. Apply a scrambling process \(\calD_{W}(\bullet):=W(\tau)(\bullet)W^{\dagger}(\tau)\), where \(W(\tau)=U^{\dagger}(\tau)W(0)U(\tau)\).

4. Apply a recovery map \(\calR_{V,\bfS}:=\mathcal{J}_{\bfS}\circ\mathcal{U}^{\dagger}_{V,\theta}\), which is given by the inverse process of the weak interaction \(\mathcal{U}_{V,\theta}\) and \(\calJ_{\bfS}(\bullet):=\sum_{j=\pm}\bra{j}\mathop{\mathrm{Tr}}_{\bfS}(\bullet)\ket{j}| {j} \rangle\langle {j} |_{\bfQ}\).

Then, we calculate the irreversibility of the above process based on Eq. 4 and post-process its results using Eq. 6 to obtain \(C_{\beta}(\tau)\). Here, we assume that \(W\) is a self-adjoint and unitary operator (\(W=W^{\dagger}, WW^{\dagger}=\openone\)) and \(V\) is a self-adjoint operator. Note that we can relax the condition of \(W\) in general, while the condition is for the ease of the experimental implementation of the evaluation (see the companion Letter [37] and its Supplemental Material). The term \(\Omega_{1/2,\pm}\) prepares the initial states \(\ket{\pm}\) with the probabilities \(1/2\) in \(\bfQ\). In the numerical experiments, we only prepare the initial state \(\ket{+}\) of \(\bfQ\) and measure the Pauli-\(x\) operator in the output state \(\rho'\) of \(\bfQ\) so that we can calculate the irreversibility \[\begin{align} C_{\beta}(\tau)=\lim_{\theta\rightarrow 0}\frac{\langle \sigma^{x} \rangle_{\rho}-\langle \sigma^{x} \rangle_{\rho'}}{2\theta^{2}}, \end{align}\] as an equivalent and more experimentally convenient way [89], [90].

Essentially, we are asking: how difficult is it to recover the initial state after it has been subjected to a weak perturbation \(V(0)\) and a scrambling process \(W(\tau)\)? The OTOC is shown to be proportional to the second-order term of the generated irreversibility. The method offers a powerful framework for evaluating the OTOC on a quantum computer. Unlike many traditional approaches that require complex interferometry or multiple-point measurements, the irreversibility-susceptibility method simplifies the experimental setup by requiring only a single final measurement on \(\bfQ\). This method makes it a highly promising and resource-efficient technique for practical implementations and can be implemented by measuring the distinguishability of quantum states, which is a natural task for quantum information processing.

4 Experiments on the trapped-ion quantum computer, reimei↩︎

To overcome the challenges of measuring the OTOC, many experiments have been conducted with various ingenious techniques [55], [91]–[111]. However, these experiments have primarily focused on analyzing the properties of the OTOC within different models, with little to no consideration given to how the experimental evaluation methods themselves influence the results. Therefore, comparing the performance of the various methods used to measure the OTOC is crucial for a deeper understanding of its intrinsic properties. In this study, we compare the RTM, WMM, and ISM to investigate the method-dependent properties of the OTOC. We conducted all experiments on the trapped-ion quantum computer—Quantinuum reimei [112], which is accessible through the Quantinuum Nexus. We note that this work represents the first experimental application of both the WMM and ISM, using the actual hardware.

4.1 Setup↩︎

We performed numerical experiments on the evaluation of the OTOC using a cloud-based quantum computer; we conducted all experiments on the trapped-ion quantum computer—Quantinuum reimei [112]. The setup for our experiment, including the preparation of the initial state, is detailed below.

4.1.1 Initial state↩︎

For the initial state, we chose to use a Gibbs state. For a given Hamiltonian \(H\) and an inverse temperature \(\beta = 1/(k_{B}T)\), the Gibbs state is defined as \[\begin{align} \rho(\beta,H)=\frac{\mathrm{e}^{-\beta H}}{Z(\beta,H)}.\label{eq:gibbs} \end{align}\tag{7}\]

Preparing such a state on a quantum computer is a non-trivial task. This is because the state is generally a mixed state, and preparing a mixed state is not possible with only unitary operations on a closed system. To overcome this challenge, we utilized a variational quantum algorithm (VQA) proposed in Ref. [113]. This algorithm prepares a purified version of the Gibbs state on an enlarged system consisting of both a target system \(\bfS\) and an ancilla system \(\bfA\). The key technique of this algorithm is based on the principle that the Gibbs state is the unique state that minimizes the Helmholtz free energy for a given Hamiltonian. The free energy is of the form \[\begin{align} F(\rho)=\mathop{\mathrm{Tr}}(H\rho)-\beta^{-1}S(\rho), \end{align}\] where \[\begin{align} S(\rho):=-\mathop{\mathrm{Tr}}(\rho\log\rho) \end{align}\] is the von Neumann entropy. The Gibbs state is the solution to the minimization problem \[\begin{align} \rho_{\text{Gibbs}}=\arg\min_{\rho}F(\rho). \end{align}\] A significant challenge in realizing the minimization on a quantum computer is the direct measurement of the von Neumann entropy \(S(\rho)\), which is not a simple observable. The VQA we employed cleverly avoids this issue by preparing a special kind of entangled state, called a thermofield double (TFD) state, on \(\bfA+\bfS\). The algorithm uses a parameterized quantum circuit (PQC) composed of two parameterized unitaries \(U_{\bfA}(\boldsymbol{\theta})\) and \(U_{\bfS}(\boldsymbol{\phi})\). The unitary \(U_{\bfA}\) acts on the ancilla qubits \(\bfA\), and \(U_{\bfS}\) acts on the target qubits \(\bfS\). These two parts are connected by a series of CNOT gates between corresponding qubits of the two registers.

Starting from the initial state \(\ket{0}^{\otimes 2n}\), where \(n\) is the number of qubits in each register, the circuit performs the following operations.

1. Apply the unitary operator \(U_{\bfA}(\boldsymbol{\theta})\) to \(\bfA\): \[\begin{align} [U_{\bfA}(\boldsymbol{\theta})\otimes\openone_{\bfS}]\ket{0}^{\otimes 2n}=\ket{\Psi(\boldsymbol{\theta})}_{\bfA}\otimes\ket{0}^{\otimes n}_{\bfS}, \end{align}\] where \(\ket{\Psi(\boldsymbol{\theta})}_{\bfA}=\sum_{i=0}^{2^{n-1}}\sqrt{p_{i}(\boldsymbol{\theta})}\ket{\psi_{i}(\boldsymbol{\theta})}_{\bfA}\) and we denote by \(\{\ket{\psi_{i}(\boldsymbol{\theta})}\}\) by the computational basis on a quantum computer. The coefficients \(\sqrt{p_{i}(\boldsymbol{\theta})}\) are determined by the parameters \(\boldsymbol{\theta}\) of the unitary operator \(U_{\bfA}(\boldsymbol{\theta})\).

2. Apply CNOT gates between \(\bfA+\bfS\). The CNOT gates entangle the two registers, creating the state \[\begin{align} \ket{\text{TFD}(\boldsymbol{\theta})}=\sum_{i=0}^{2^{n-1}}\sqrt{p_{i}(\boldsymbol{\theta})}\ket{\psi_{i}(\boldsymbol{\theta})}_{\bfA}\ket{\psi_{i}(\boldsymbol{\theta})}_{\bfS}. \end{align}\] This state is the TFD state in the computational basis.

3. Apply the unitary operator \(U_{\bfS}(\boldsymbol{\phi})\) to \(\bfS\): \[\begin{align} \ket{\widetilde{\text{TFD}}(\boldsymbol{\theta},\boldsymbol{\phi})}&=[\openone_{\bfA}\otimes U_{\bfS}(\boldsymbol{\phi})]\ket{\text{TFD}(\boldsymbol{\theta})}\nonumber\\ &=\sum_{i=0}^{2^{n-1}}\sqrt{p_{i}(\boldsymbol{\theta})}\ket{\psi_{i}(\boldsymbol{\theta})}_{\bfA}\ket{E_{i}(\boldsymbol{\theta},\boldsymbol{\phi})}_{\bfS}. \end{align}\] Here, we assume that \(U_{\bfS}(\boldsymbol{\phi})\) transforms the computational basis states \(\ket{\psi_{i}(\boldsymbol{\theta})}\) to the eigenstates \(\ket{E_{i}(\boldsymbol{\theta},\boldsymbol{\phi})}\) of \(H\).

By tracing out the register \(\bfA\) from the final TFD state, we obtain a mixed state on \(\bfS\): \[\begin{align} \rho_{\bfS}(\boldsymbol{\theta},\boldsymbol{\phi})&=\mathop{\mathrm{Tr}}_{\bfA}[| {\widetilde{\text{TFD}}(\boldsymbol{\theta},\boldsymbol{\phi})} \rangle\langle {\widetilde{\text{TFD}}(\boldsymbol{\theta},\boldsymbol{\phi})} |]\nonumber\\ &=\sum_{i=0}^{2^{n-1}}p_{i}(\boldsymbol{\theta})| {E_{i}(\boldsymbol{\theta},\boldsymbol{\phi})} \rangle\langle {E_{i}(\boldsymbol{\theta},\boldsymbol{\phi})} |_{\bfS}, \end{align}\] where \(p_{i}(\boldsymbol{\theta})=\) The crucial insight of this algorithm is that we can get the probability distribution \(p_{i}(\boldsymbol{\theta})\) to calculate \(S(\rho)\) by measuring \(bfA\) in the computational basis without directly measuring \(\bfS\), i.e., applying the tracing out operation over \(\bfS\) to \(\ket{\text{TFD}(\boldsymbol{\theta})}\) yields a mixed state of \(\bfA\): \[\begin{align} \rho_{\bfA}(\boldsymbol{\theta})&=\mathop{\mathrm{Tr}}_{\bfS}[| {\text{TFD}(\boldsymbol{\theta})} \rangle\langle {\text{TFD}(\boldsymbol{\theta})} |]\nonumber\\ &=\sum_{i=0}^{2^{n-1}}p_{i}(\boldsymbol{\theta})| {\psi_{i}(\boldsymbol{\theta})} \rangle\langle {\psi_{i}(\boldsymbol{\theta})} |_{\bfA}. \end{align}\] This property means that the von Neumann entropy \(S(\rho_{\bfS})\) of \(\bfS\) is equal to the von Neumann entropy \(S(\rho_{\bfA})\) of \(\bfA\). Since \(\rho_{\bfA}\) is a diagonal state in the computational basis, its entropy can be calculated simply from \(\{p_{i}(\boldsymbol{\theta})\}\). Thus, we can avoid the need for complex task—direct entropy measurements on the quantum computer.

The cost function of the VQA is therefore defined as the free energy of \(\bfS\), which is minimized by optimizing the parameters \(\boldsymbol{\theta}\) and \(\boldsymbol{\phi}\): \[\begin{align} F(\boldsymbol{\theta},\boldsymbol{\phi})=\mathop{\mathrm{Tr}}[H\rho_{\bfS}(\boldsymbol{\theta},\boldsymbol{\phi})]-\beta^{-1}S(\rho_{\bfA}(\boldsymbol{\theta})). \end{align}\] The first term can be measured on \(\bfS\). The second term can be calculated classically from \(\{p_{i}(\boldsymbol{\theta})\}\). The VQA iteratively optimizes the parameters to minimize \(F(\boldsymbol{\theta},\boldsymbol{\phi})\) until convergence. Upon finding the optimal parameters \(\boldsymbol{\theta}^{*}\) and \(\boldsymbol{\phi}^{*}\), the circuit prepares the Gibbs state \(\rho_{\text{Gibbs}}\approx\rho_{\bfS}(\boldsymbol{\theta}^{*},\boldsymbol{\phi}^{*})\) on \(\bfS\).

4.1.2 Observables↩︎

In our experiment, we chose the one-dimensional Heisenberg XXZ model in a transverse magnetic field as the generator of dynamics for our system of interest. This is a canonical model in condensed matter physics known for its rich phase diagram, which includes both integrable and chaotic regimes. The Hamiltonian is given by \[\begin{align} H=-\frac{1}{4}\sum^{n-1}_{i=1}(\sigma^{x}_{i}\sigma^{x}_{i+1}+\sigma^{y}_{i}\sigma^{y}_{i+1}+\Delta\sigma^{z}_{i}\sigma^{z}_{i+1})-h\sum^{n}_{i=1}\sigma^{z}_{i},\label{eq:xxz} \end{align}\tag{8}\] where \(\sigma^{\bullet}_{i}\) for \(\bullet=x,y,z\) are the Pauli matrices acting on the \(i\)-th qubit, \(\Delta\) is the anisotropy parameter, and \(h\) is the strength of the transverse magnetic field. For the present paper, we set the magnetic field strength to \(h=(1-\Delta)/2\).

Our numerical experiments were conducted with the following specific parameters:

The size of \(\bfS\)

\(\bfS\) consisted of \(n=4\) qubits.

The parameters of \(H\)

We vary the anisotropy parameter \(\Delta\) from \(0.1\) to \(0.9\): \(\Delta=\{0,1,0.3,0.5,0.7,0.9\}\).

The time evolution

We measured the evolution of the OTOC from time \(t=0\) up to \(t=2.1\). This time interval was discretized into \(15\) time steps.

Observables

For the evaluation, we selected local Pauli operators as our observables. Specifically, we chose both \(W(0)\) and \(V(0)\) to be the Pauli-\(x\) operator on the first qubit: \(W(0)=V(0)=\sigma^{x}_{1}\). In terms of the full system, this corresponds to the operator \(X\otimes\openone\otimes\openone\otimes\openone\).

The value of the measured OTOC at each time step was averaged over \(10\) independent experimental iterations to mitigate statistical noise. For the statistical averaging, each quantum circuit was executed with a total of \(1000\) measurement shots. In the ISM, the interaction strength was set to \(\theta=0.4\) radians. The initial state for our experiment was prepared as an approximation of the Gibbs state defined in Eq. 7 , with the XXZ Hamiltonian specified in Eq. 8 .

4.2 Results↩︎

As shown in Fig. 4, the solid lines represent the ideal values obtained from numerical calculation. The circle dots show the values from noiseless simulations via aer-simulator, while the square dots represent the experimental data obtained from the reimei. The error bars indicate the standard deviation of the measured values.

Overall, a consistent observation across all the results obtained from three methods is the subtle deviation between the simulation and experimental results, which is primarily attributed to errors from extrinsic noise. In addition, the discrepancies between the theoretical and the simulation results are likely due to a combination of two factors: fidelity errors from potential imperfections in the preparation of the variational Gibbs state, and algorithmic errors from finite-value truncation in the Trotter decomposition. Despite these inaccuracies, the results clearly demonstrate the feasibility of experimentally exploring quantum dynamics on current quantum computers.

Next, we compare the results for each method in detail. A notable difference among the methods is the variation in the behavior of the measured OTOC values with respect to the parameter \(\Delta\) in the XXZ Hamiltonian. Figure 4 (a) shows the OTOC evaluated using the RTM. At early times, the results align well with the ideal values regardless of the values of \(\Delta\). However, as time progresses, the measured values for \(\Delta = 0.5\), \(0.7\), and \(0.9\) deviate downward from the theoretical predictions. Figure 4 (b) shows the OTOC evaluated using the WMM. Similar to the RTM, the results match the ideal values at early times. In the later time regime, however, the values for \(\Delta = 0.5\), \(0.7\), and \(0.9\) deviate upward from the theoretical predictions. Figure 4 (c) shows the OTOC evaluated using the ISM. For all values of \(\Delta\), the measured OTOC is in good agreement with the theoretical values on average. Because of the use of a weak interaction, the standard deviation for the ISM is larger compared to the other two methods. However, we must carefully interpret the apparent robustness of the ISM against variations in the anisotropy parameter \(\Delta\). As observed in Fig. 4 (c), the ISM exhibits larger standard deviations compared to the RTM and WMM. This variance is not a hardware artifact, but an intrinsic feature of protocols relying on weak interactions. Specifically, the useful signal extracted in the ISM—the shift in the expectation value of the ancilla’s Pauli-\(x\) operator—scales quadratically with the weak interaction strength: \(1 - \langle \sigma^{x} \rangle_{\rho'} \approx 2\theta^{2}C_{\beta}(\tau)\). Conversely, the shot noise for \(N\) measurement shots scales as \(1/\sqrt{N}\). Consequently, the signal-to-noise ratio (SNR) is inherently suppressed by the weak coupling parameter, yielding \(\text{SNR} \sim \mathcal{O}(\theta^{2}\sqrt{N})\). In our experiment, with \(\theta = 0.4\) radians and \(N = 1000\) shots, the statistical noise floor remains relatively high. Therefore, it is plausible that any subtle, method-dependent deviations induced by the physical parameter \(\Delta\)—which are clearly resolved in the \(\mathcal{O}(1)\) signals of the RTM and WMM—are currently being masked by the dominant statistical uncertainty in the ISM. While the mean values of the ISM faithfully track the ideal values across all tested regimes, this fundamental variance penalty precludes a definitive conclusion regarding the absence of \(\Delta\) dependence. Overcoming this limitation in future experiments will require significantly increasing the shot count \(N\) to compensate for the \(\theta^{2}\) suppression.

In summary, while the ISM shows little dependence on the parameter \(\Delta\), both the RTM and WMM exhibit reduced accuracy in the later time regime for larger values of \(\Delta\). This method-dependent behavior of the measured OTOC properties has not been previously reported. Further investigation and detailed analysis are needed to understand the underlying causes of this phenomenon.

5 Conclusion↩︎

In this study, we have presented the numerical evaluations of the OTOC on the reimei using three distinct protocols: the rewinding time method (RTM), the weak-measurement method (WMM), and the irreversibility-susceptibility method (ISM). By successfully implementing these methods on the reimei with a variationally prepared Gibbs state of the XXZ Heisenberg model, we have demonstrated the practical feasibility of exploring complex quantum dynamics on current-generation hardware. Our experimental results show a good overall agreement with theoretical predictions, validating these methods as practical tools for investigating quantum chaos.

A central aspect of our work was the comparative analysis of these three experimental schemes. We found that each method has unique strengths and limitations. The RTM is conceptually straightforward but requires precise measurement of two observables. The WMM offers physical insight into the OTOC based on quasiprobability, though it can be demanding in terms of the number of weak measurements required. The ISM provides a direct link between scrambling and irreversibility, but its performance is highly sensitive to the gate fidelity of the weak interactions.

The significance of our findings is twofold. First, we have provided the first experimental demonstration of the ISM for measuring the OTOC in a non-trivial many-body system, thereby expanding the toolkit of protocols for studying quantum information scrambling. Second, our work contributes to the growing evidence that current quantum processors are effective platforms for exploring complex, fundamental phenomena in quantum information science.

Our results also highlight a previously unreported, method-dependent behavior of the measured OTOC with respect to the XXZ Hamiltonian parameter \(\Delta\), particularly for the RTM and WMM. Future work should focus on clarifying the underlying cause of this dependence. Furthermore, to probe the onset of many-body chaos more directly, it is essential to scale these experiments to a larger number of qubits and investigate how the OTOC depends on system size. Applying these methods to different Hamiltonians would also provide deeper insights into the nature of quantum dynamics. Finally, improving the fidelity of variational Gibbs state preparation and incorporating advanced error mitigation techniques will be crucial steps toward achieving quantitatively precise OTOC measurements on near-term quantum hardware.

References↩︎