Diagnosing quantum chaos in many-body systems using entanglement as a resource

Classical chaotic systems exhibit exponentially diverging trajectories due to small differences in their initial state. The analogous diagnostic in quantum many-body systems is an exponential growth of out-of-time-ordered correlation functions (OTOCs). These quantities can be computed for various models, but their experimental study requires the ability to evolve quantum states backward in time, similar to the canonical Loschmidt echo measurement. In some simple systems, backward time evolution can be achieved by reversing the sign of the Hamiltonian; however in most interacting many-body systems, this is not a viable option. Here we propose a new family of protocols for OTOC measurement that do not require backward time evolution. Instead, they rely on ordinary time-ordered measurements performed in the thermofield double (TFD) state, an entangled state formed between two identical copies of the system. We show that, remarkably, in this situation the Lyapunov chaos exponent $\lambda_L$ can be extracted from the measurement of an ordinary two-point correlation function. As an unexpected bonus, we find that our proposed method yields the so-called"regularized"OTOC -- a quantity that is believed to most directly indicate quantum chaos. According to recent theoretical work, the TFD state can be prepared as the ground state of two weakly coupled identical systems and is therefore amenable to experimental study. We illustrate the utility of these protocols on the example of the maximally chaotic Sachdev-Ye-Kitaev model and support our findings by extensive numerical simulations.


I. INTRODUCTION
A key characteristic of chaotic quantum many-body systems is rapid dispersal of quantum information deposited among a small number of elementary degrees of freedom. After a short time the information is distributed among exponentially many degrees of freedom, whereby it becomes effectively lost to all local observables. This apparent loss of quantum information through unitary evolution, known also as "scrambling", lies at the heart of thermalization in closed systems and plays a key role in understanding quantum aspects of black holes as epitomized by Hawking's information loss paradox [1]. Black holes are believed to scramble at the fastest possible rate consistent with causality and unitarity [2]. Some strongly coupled quantum systems, such as the Sachdev-Ye-Kitaev (SYK) model [3][4][5][6] and its variants, are also known to be fast scramblers; this motivates their description as duals of gravitational theories containing a black hole. Scrambling in quantum theories can be quantified through the out-of-time order correlators defined below, which -for chaotic systems -show exponential growth at intermediate times with a characteristic Lyapunov exponent λ L . A universal upper bound on chaos, conjectured by Maldacena, Shenker and Stanford [7], posits that λ L ≤ 2πT and is saturated in the class of maximally chaotic systems which includes black holes and SYK models.
Diagnosing quantum chaos and scrambling in realistic * lantagne@phas.ubc.ca † plugge@phas.ubc.ca physical systems is a problem of fundamental importance that has been only partially addressed so far. As we review below, the great hurdle in conventional approaches to diagnosing chaos is the necessity to evolve the quantum system backward in time during the measurement [8][9][10]. So far this has been achieved in a very limited range of systems, mainly quantum simulators [11] and ion traps [12,13]. However there is no hope of applying this method to a broader class of naturally occurring quantum many-body systems, because in those one simply does not possess the level of control required to reverse the time evolution. In this paper we address this pressing challenge by introducing a new approach to diagnosing chaos in quantum many-body systems which does not require backward time evolution during the measurement. The approach is based on a procedure which creates an entangled resource state that permits chaos diagnosis through an ordinary measurement. Specifically, we show that in this resource state the chaos exponent λ L governs the exponential growth of an ordinary two-point correlator at intermediate times and can thus be experimentally accessed using routine spectroscopic techniques.
Rather than a Loschmidt echo, our scheme more closely resembles the approach of Refs. [14][15][16] to the measurement of Renyi entropy by readout of entangling operators between identical copies of a quantum state. Alternative approaches for the detection of scrambling in quantum systems include interferometry [17] or outof-equilibrium measurement protocols [18,19]. In our case, the entanglement between identical copies of the chaotic quantum system is generated as resource from a specifically engineered Hamiltonian. The ability to detect OTOCs hence emerges from measurements of simple arXiv:1907.01628v3 [cond-mat.str-el] 5 Mar 2020 operators in an a-priori complicated ground state of two coupled chaotic systems.
A quantitative measure of scrambling in quantum systems is given by the expectation value of a commutator squared [20], of two initially commuting Hermitian operators [W, V ] η = 0. In the following we allow for both bosonic (η = 1) and fermionic (η = −1) statistics, with [ · , · ] η denoting the commutator (anti-commutator) for η = 1 (η = −1). The operators evolve in time according to the system Hamiltonian H through W (t) = e iHt W e −iHt and ... denotes the thermal average at inverse temperature β = 1/T . The intuition behind the definition (1) is the following: As the operator W (t) evolves in time it becomes more and more complex until it eventually fails to commute with operator V . One thus expects C η (t) to grow as a function of time and eventually saturate at a value close to 2 V 2 W 2 for large t, regardless of the specific form of V and W . In chaotic many-body systems, at intermediate times, the growth of C η (t) follows an exponential dependence, as long as V and W are "simple" operators composed of products of a small number of elementary degrees of freedom. Expanding the commutator in Eq. (1) we obtain two types of thermal averages The averages on the first line represent naturally timeordered correlators (NTOC) that correspond to "sensible" experiments performed in an ordinary quantum system. For instance the second term V W (t)W (t)V describes a process in which we perturb the system at time t = 0 by applying operator V , evolve the perturbed system forward in time and then perform a measurement of a quantity represented by operator W 2 .
The averages on the second line of Eq.
(2) represent out-ot-time-ordered correlators. These correspond to less sensible experiments. The second term for instance, which is often denoted as describes the process in which we compare two states of the system: one obtained by first perturbing with V , then applying W (t) at a later time t; the other obtained by perturbing first with W (t), evolving backward in time, then applying V . We note that in the special case where operators V and W are also unitary, the quantity C η (t) can be expressed simply as The fundamental difference between NTOC and OTOC is best visualized by placing the operators on the Schwinger-Keldysh contour illustrated in imaginary time evolution is generated by the density matrix which becomes explicit if we rewrite Eq. (3) as F (t) = tr[W (t)V W (t)V e −βH ]. Notice that while NTOC can be represented by placing the operators on a conventional Schwinger-Keldysh contour with one forward and one backward evolving branch (Fig. 1a), OTOC require a doubled contour indicated in Fig. 1b.
As a rule, the physical implementation of an OTOC measurement requires (actual or effective) backward time evolution and is therefore difficult to achieve in most systems. Similar to Loschmidt echo [21][22][23][24], a measurement of F (t) is possible in situations where one controls the Hamiltonian at the microscopic level and can, in particular, reverse its sign to generate backward time evolution. As a practical matter this restriction greatly limits the types of systems in which the phenomenon of scrambling can be experimentally probed.
In the rest of the paper we introduce, discuss, and put to the test a family of protocols that probe OTOCs but do not require backward time evolution. They instead require two copies of the system prepared in a special entangled state called "thermofield double" (TFD). As we explain in detail below, TFD is a pure quantum state whose reduced density matrix coincides with the thermal density matrix of one copy of the system. The TFD state has been widely studied in quantum gravity theories as a description of traversable wormholes [25][26][27][28]. It has a remarkable property, which we review below, that time effectively flows in the opposite direction in two entangled systems. It is this property which underlies its usefulness in the proposed OTOC measurement proto-cols. Importantly, recent theoretical work has established a simple method that can be used to prepare the TFD state [29]. The method consists of weakly coupling two identical subsystems in a specific way, then cooling the combined system to its ground state. As an example, the ground state wavefunction of two identical SYK Hamiltonians coupled by simple bilinear tunneling terms has > 96% overlap with a TFD wavefunction [27] (see also Fig. 2a). We also note interesting recent works on how to prepare a TFD state using quantum circuits [30][31][32]; however, these approaches are limited to the moderate system sizes accessible on present-day quantum simulators, similar to the experiments of Refs. [11][12][13].
In the following we first review the concept of the TFD state, and then demonstrate how it can be used to diagnose quantum chaotic behavior via an OTOC measurement that does not require explicit backward time evolution of quantum states (Sec. II) We discuss in detail how the TFD state can be prepared and used to extract the chaos exponent λ L from an equilibrium measurement of a two-point correlation function. We then apply these general ideas to a pair of coupled SYK Hamiltonians, recently argued to be holographically dual to a traversable wormhole, and known to admit a TFD ground state (Sec. III). This simple model serves as a testbed for demonstrating the usefulness of our OTOC measurement protocols. Finally, in Sec. IV we discuss possible physical realizations of the coupled SYK models in the laboratory, and expand on the challenges of performing the necessary measurements. We conclude with an outlook onto interesting future work and outstanding challenges in Sec. V.

II. OTOC MEASUREMENT USING THE THERMOFIELD DOUBLE STATE
In this Section we review the concept of the thermofield double state, discuss its properties, and then show how it can be used to measure OTOCs in a way that does not require explicit backward time evolution.

A. TFD state: Definition and properties
Consider two copies of the same system, left and right, described by many-body Hamiltonians H L and H R , respectively. We assume that H α (α = L, R) are invariant under time reversal generated by an antiunitary operator Θ. The TFD state at inverse temperature β is then defined as where |n α is an eigenstate of H α with energy eigenvalue E n , and Z β = n e −βEn is the partition function. |n = Θ|n denotes the time-reversed partner of the eigenstate |n which shares the same energy eigenvalue E n . We note that time reversal is necessary here to define a unique TFD state, as each eigenstate |n is defined up to an overall phase e iφn . A direct product |n L ⊗ |n R is however well-defined because |n transforms with the opposite phase to |n under the corresponding U (1) transformation. In the limit of zero temperature |TFD β simply becomes a direct product of L and R ground states, whereas at infinite temperature it becomes a maximally entangled state between the two subsystems.
The TFD state has several important properties. The expectation value of any one-sided operator with respect to |TFD β is given by a thermal average, It is also important to note that |TFD β is not an eigenstate of the full system Hamiltonian H = H L + H R . It is, however, an eigenstate with eigenvalue zero of H − = H L − H R ; it can be easily checked that This has implications for the concept of time-translation invariance in the TFD state. Eq. (7) implies that |TFD β evolves trivially under H − , Hence the expectation value of a product of two operators acting in the L and R systems has the property valid for arbitrary t. The second line follows upon replacing |TFD β in the expectation value on the first line by e −it(H L −H R ) |TFD β using Eq. (8), and recalling that O R commutes with H L (and same for L ↔ R). Choosing t = t 2 we see that F(t 1 , t 2 ) is a function of t 1 + t 2 only. This should be compared to the statement of timetranslation invariance in a conventional (unentangled) state, where the expectation value only depends on t 1 − t 2 as long as the Hamiltonian is independent of time.
In the context of wormhole physics, Eq. (9) can be interpreted as time flowing in the opposite direction on the two sides of the wormhole, represented in the quantum theory by two entangled subsystems. It is this peculiar property that ultimately allows one to use a TFD as a resource for OTOC measurement without explicit backward time evolution.

B. Probing OTOC using TFD state
For the purposes of this subsection, we will assume that we have the ability to engineer two identical copies of an interesting quantum many-body system described by Hamiltonians H L and H R and prepare them in the TFD state Eq. (5). In the subsequent Sections, we will discuss how this can be achieved in practice, and study some concrete examples. Here we focus on elucidating how a two-sided measurement performed on the TFD state can be used to probe correlation functions that map onto thermal OTOCs with respect to one subsystem.
Consider a naturally time-ordered correlator evaluated with respect to the TFD state in Eq. (5). Here T denotes the time-ordering operator. Normal timeordered expectation values of this type correspond to physical quantities that are measurable, at least in principle. The specific average defined in Eq. (10) can be thought of as a component of the current-current correlator (where the current operator involves both sides of the composite system) that would arise in the calculation of the appropriate linear response conductance. We now assume t > t and write the average in Eq. (10) explicitly using the TFD state (5). We obtaiñ where we used the fact that L and R operators (anti-)commute at all times, and that they act only on L and R eigenstates, respectively. The right hand side has been written as a product of expectation values taken in L and R systems separately. Because each such expectation value is a complex number and the two subsystems are identical, we can now drop the L and R subscripts, recognizing that it does not matter in which subsystem they are evaluated [33] We thus havẽ where we used general properties of time-reversed states [34], namely ā|b = Θa|Θb = b|a , and As the final step we insert the Boltzmann factors into the expectation values, and replace them by powers of the density matrix, e.g. e −β(En+Em)/2 n|W (t)V (t )|m = Z β n|y 2 W (t)V (t )y 2 |m , with This allows us to perform the sum over n using the completness relation n |n n| = 1 and arrive at the result Note that the trace is now evaluated with respect to the eigenstates of a single-sided Hamiltonian.
For any t > 0 and t < 0 Eq. (15) has the structure of an OTOC. In the special case t = −t it becomes where we used the time-translation invariance to shift all temporal arguments by +t for clarity. We observe thatF (t, −t) coincides with the canonical OTOC function F (2t) defined in Eq. (3), except for the placement of the density matrix powers y (see also Fig. 1c). Expressions of this type are called "regularized" OTOCs and have been extensively studied in the literature [5,7]. Regularized OTOCs exhibit less singular behavior than unregularized OTOCs when evaluated analytically, and they have been argued to more reliably measure quantum chaos in many-body systems [35][36][37]. The universal upper bound on the Lyapunov exponent λ L has also been proven only for regularized OTOCs [7]. We see that, remarkably, naturally time-ordered correlators evaluated in the TFD state map onto regularized OTOCs with respect to the single system.
Following the sequence of steps between Eqs. (10) and (15) it is possible to derive other useful identities that relate expectation values of operators in the TFD state to single-sided OTOCs. For instance we find The first line can be interpreted as creating an excitation in the TFD state at time −t, evolving forward in time and performing a two-sided measurement of a quantity represented by the operator V L V R at time zero. This correlator also maps onto the canonical OTOC albeit with a different regularization, illustrated in Fig. 1d. Below we refer to this as "asymmetric" regularized OTOC. We note that similar relations have been anticipated in the high-energy community [7,20,38].

C. Initial state preparation
Our considerations above establish formal identities relating naturally time-ordered correlators evaluated in an entangled state of two identical systems to out-of-timeordered correlators in the single system, such as Eqs. (16) and (17). A sensible measurement performed in the TFD state can thus provide information on the OTOC and diagnose quantum chaotic behavior in a many-body system. We now discuss a method that allows one to prepare the TFD state. In addition, we see from the discussion in the previous subsection that in order to probe the OTOC one in fact needs |TFD β at a negative time, for only then is the correlator on the left hand side of Eq. (17) naturally time ordered. The same remark applies toF (t, t ) defined in Eq. (10). In the following we therefore discuss how a state closely approximating |TFD β (−t) can be prepared in a realistic setup.
The easiest way to prepare a TFD state in the laboratory would be to engineer a Hamiltonian H S which admits |TFD β as a ground state. A collection of such Hamiltonians were recently constructed in Ref. [29], with a unique |TFD β ground state separated from the rest of the spectrum by a gap of order 1/β. Thus, a TFD state can be prepared by engineering the system to obey Hamiltonian H S , and then cooling it down to a physical temperature T phys small compared to the gap. The form of H S required to obtain the TFD ground state exactly is complicated and therefore not practical from the standpoint of generating the state in a laboratory. However, the ground state |Ψ 0 of a simple Hamiltonian with and Θ representing the time-reversal operator, has been shown [29] to approximate |TFD β to good accuracy for appropriately chosen coefficients c j . Here O L/R j are arbitrary (but identical) operators acting on the L and R systems, respectively. This method is expected to apply to generic many-body Hamiltonians H L respecting the eigenstate thermalization hypothesis [29], and is thus of broad relevance in the study of quantum chaotic systems.
In another recent work, Maldacena and Qi [27] used an even simpler construction with a coupling of the form to describe an eternal traversable wormhole formed by two copies of the SYK model. We will review this construction in the next Section, verify numerically that it admits a ground state that closely approximates the TFD state, and discuss some of its intriguing properties. The above procedure allows one to prepare a state that closely approximates |TFD β by coupling two systems through H I defined in Eqs. (19) or (20), then cooling the combined system to reach its ground state. When the coupling H I is switched off the system begins to evolve forward in time according to the decoupled Hamiltonian This evolution is non-trivial because |TFD β is not an eigenstate of H 0 . In order to probe OTOC, however, we require a TFD state prepared at negative time −t, as discussed above. It would thus appear that our protocol requires backward time evolution after all. We show in Appendix A, that the required resource state |Ψ 0 (−t) |TFD β (−t) for short time durations t can be prepared by manipulating the strength of the coupling H I , without the need to reverse the sign of H 0 . Since the ability to introduce and control H I is necessary to prepare the TFD state in the first place, this method does not introduce any substantial additional complications.

D. OTOC from two-point functions
Here we discuss an approach that allows one to extract the OTOC from the measurement of a time-ordered twopoint function G LR (t, t ), in a generic chaotic system with coupling µ between the two sides given by Eq. (20). This method relies only on the fact that the ground state of the coupled system closely approximates the TFD state and, importantly, does not require varying µ before or during the measurement. In contrast to related previous work [39,40], our method comprises the measurement of a single (averaged) Green's function, rather than statistics on an ensemble of measurements. We outline the argument below and provide technical details in Appendix B.
We consider the two-point time-ordered LR correlation function in real time, The average is taken with respect to the (TFD) ground state of the coupled system and is therefore time- where V = V (0), y represents the fourth root of the thermal density matrix, Eq. (14), and the trace is performed with respect to the many-body eigenstates |n of H L . Operators O j enter through the coupling H I which is assumed to have the Maldacena-Qi form Eq. (20). The trace on the first line is a time-independent constant, and the trace on the second line is a naturally time-ordered four-point correlator. Crucially, the trace on the third line has the structure of a regularized OTOC for all s inside the integration bounds. Near the lower bound s → 0 it coincides with the regularized OTOC. For s = 0 the trace represents a more general form of the OTOC dependent on two time variables, F (t 1 , t 2 ) = tr[W (t 1 )V y 2 W (t 2 )V y 2 ]. In chaotic systems the latter is commonly assumed to behave according to with b(t) an even function of t and b(0) real positive [36,41,42]. The trace on the last line of Eq. (23) is then related to F (t + s, t − s) and, upon integration, becomes proportional to B(t)e λ L t with B(t) = µ t 0 ds b(2s). If we adopt another common assumption [5] that due to their exponential growth at intermediate times OTOCs tend to dominate over NTOCs, we may conclude that the intermediate-time behavior (J −1 t µ −1 ) of the LR two-point correlator of the coupled theory should be well approximated by where A and B are slowly varying functions of t and may be taken as constants.
Eq. (25) indicates that, remarkably, the Lyapunov exponent characteristic of a single chaotic system at inverse temperature β can be extracted by measuring the LR causal two-point correlator in the ground state of two identical such systems, coupled through static bilinear terms as in Eq. (20). We remark that an expression similar to Eq. (23) can be derived for the retarded twopoint correlator G ret LR (t, −t) which also contains a dominant OTOC contribution at intermediate times. Retarded correlators are often more directly related to measurable quantities, and we will employ them in Sec. III where we provide an explicit numerical calculation for the example of coupled SYK models. This shows approximate exponential growth of the retarded two-point correlator in the appropriate time interval, which lends support to the conclusions reached in this subsection.

E. Discussion and caveats
Our main results obtained in this Section can be summarized as follows. We showed that naturally time ordered correlators, such as the one defined in Eq. (10), evaluated in the TFD state map to out-of-time-ordered correlators with respect to a single system. This result is generic and requires only that the correlator in question involves operators from both L and R side of the coupled system. (Details of the correlators however change the regularization structure, that is, the placement of y 2 factors in the corresponding OTOCs.) The transmutation from a two-sided NTOC to a single-sided OTOC can be intuitively understood as a consequence of the fact that, effectively, time flows in the opposite direction in the two subsystems forming the TFD. Mathematically this unusual property follows from Eq. (9), which shows that a generic two-sided correlator with respect to the TFD state depends on the sum of the temporal arguments for each subsystem, and not on their difference as would normally be the case.
Further, we argued that ordinary Green's functions should also capture the behavior of OTOCs for small values of the coupling µ between the subsystems and short times |t − t | µ −1 . Crucially, such two-point correlation functions are in principle much easier to probe in the laboratory, because the measurement can be performed under equilibrium conditions and the protocol does not require varying any system parameters. We shall discuss a specific example of this in Sec. IV-B.

III. APPLICATION: TWO COUPLED SYK MODELS
In this Section, we apply the ideas presented above to a concrete model recently introduced by Maldacena and Qi [27] that realizes a quantum-mechanical dual to an eternal traversable wormhole in (1+1)-dimensional anti-de Sitter spacetime (AdS 2 ) by coupling two identical SYK models. This model is convenient for us for several reasons. First, the SYK model is known to be maximally chaotic: its OTOC exhibits exponential growth with an exponent that saturates the universal chaos bound [5,7]. Second, the ground state of two such coupled SYK models is well approximated by the TFD state [27], which allows us to apply the machinery developed in the previous Section in a relatively simple setting. Third, there are several proposals in the literature for experimental realizations of the SYK model and its variants [43][44][45][46][47], making it a potentially fruitful platform for laboratory explorations.

A. The model
The model introduced by Maldacena and Qi in Ref. [27] has the form where µ is a constant, and H SYK α with α = (L, R) describe two identical SYK models each involving N Majorana zero-mode operators. These obey the usual algebra The coupling constants J ijkl are random, independent Gaussian variables respecting and are independent of α -that is, the disorder in both SYK models is perfectly correlated. Referring to a potential experimental realization of the Maldacena-Qi model using quantum dots (see Sec. IV below), and for the sake of brevity, we henceforth dub the two SYK subsystems as L and R "dots". Without loss of generality, we define a complex fermion basis as to construct the many-body Hilbert space. This basis is helpful because the anti-unitary time-reserval symmetry of the model is manifest and simply represented by Θ = K, where K denotes complex conjugation. The Majorana operators then transform as Θχ j L Θ −1 = χ j L and Θχ j R Θ −1 = −χ j R . With such a choice H becomes purely real in the many-body basis defined by the operators c j . The coupled SYK system is also invariant under the total fermion parity where N f = j c † j c j is the total fermion number. What is less obvious is that the fermion number modulo 4, Q 4 ≡ (N f ) mod 4, is also a symmetry of H. This property relies on the perfectly correlated disorder between the two subsystems [48].
The model can be solved in the limit of large N by methods developed in the context of the original SYK model [3][4][5]. This involves formulating the theory as a Euclidean-space path integral, averaging over the disorder using the replica formalism, and finally writing the large-N saddle-point action for the averaged fermion . The resulting saddle-point action reads [27] where S 0 = −N ln Pf(δ αβ ∂ τ − Σ αβ ) and Σ αβ denotes the self energy associated with G αβ . The corresponding saddle-point equations are obtained by varying the action with respect to G αβ and Σ αβ . Using the time-translation invariance, so that G αβ (τ 1 , τ 2 ) = G αβ (τ 1 − τ 2 ), and the mirror symmetry between the L and R subsystems, one can write the saddle-point equations in terms of two independent correlators G LL (τ ) and G LR (τ ). Their frequency-space counterparts are given as 2 and ω n = πT (2n + 1) is the nth Matsubara frequency. The self energies are given by In the following, we support the ideas for OTOC measurement using the TFD state presented in Sec. II by analyzing numerical solutions of the Maldacena-Qi model defined by Eqs. (26) and (27). We perform exact diagonalizations of the Hamiltonian for systems sizes 2N as large as 32, and use the results to calculate various quantities of interest. We also numerically solve the large-N saddle-point equations (32) and (33) by analytically continuing to real time and frequency domain, and then employing the iterative procedure described in Refs. [5,49]. This yields retarded propagators G ret αβ (ω) and their time domain counterparts. Some useful analytical simplifications of the above Schwinger-Dyson (SD) equations and details of our numerical procedures are described in Appendix F.

B. TFD ground state
In Ref. [27] it was argued that the model in Eq. (26) admits an approximate TFD ground state for all values of the dimensionless parameter µ/J, and an exact TFD ground state in the limits of either small or large µ/J. This can be understood intuitively as follows. For µ/J → 0 the ground state of the system is simply given by |0 L ⊗ |0 R , which coincides trivially with the zerotemperature TFD state |TFD ∞ . For µ/J → ∞ the system is best understood as a collection of N decoupled two-level systems with the Hamiltonian given by the last term in Eq. (26), and energy levels ±µ corresponding to the presence or absence of a fermion in that state. The many-body ground state |Ψ 0 of this system is unique and such that all N single-particle states are empty, where the c j are defined in Eq. (29). This state is equivalent (see Ref. [48] for an explicit proof) to the infinitetemperature TFD state where |n are the eigenstates of H SYK . For intermediate values of µ/J one must resort to numerical exact diagonalization, which confirms that the TFD state is always a good approximation to the true ground state of the system, as summarized in Fig. 2a (see also Refs. [27,48]). The overlap is always greater than 0.96, with the minimum occurring around µ/J ∼ 0.1. This minimum was argued to indicate a phase transition of the Hawking-Page type [50] between a wormhole phase at small µ/J and low temperature, and a black hole phase at large temperature [27,48]. The parameter β characterizing |TFD β which best describes the ground state is monotonically decreasing as a function of µ/J (see Fig. 2b). The energy gap to the first excited state, displayed in Fig. 2c, scales as the temperature of the TFD state 1/β, as expected from the arguments of Ref. [29]. In Sec. III C we obtain the constant of proportionality as E gap ≈ 1.3T from the large-N  solution, which agrees well with the ED numerics. However, as shown in Fig. 2f, our ED calculation does not show the scaling E gap ∼ µ 2/3 J 1/3 at small µ/J, expected from the wormhole duality and confirmed by solving the imaginary-time SD equations (32) and (33) in Ref [27]. This is presumably due to finite-size effects which become important at energy scales smaller than ∼ J/N .
We can extract the gap amplitude more precisely from the numerical solution of the large-N saddle-point equations (32) and (33), but now solved in real time and frequency domain. This is most easily done by analyzing the spectral function defined using the retarded propagator G ret LL (ω), which is related to the Matsubara frequency propagator G LL (iω n ) by the standard analytical continuation iω n → ω + iδ [51]. The spectral function is shown for several values of µ in Fig. 2d. The spectral gap E gap , defined here as the position of the first peak in A(ω), is plotted in Fig. 2f. It shows E gap = µ 2/3 J 1/3 scaling (with numerical prefactor very close to 1) for small µ/J, with a crossover to a linear dependence occurring around µ/J ≈ 0.1. An extensive symmetry analysis and substantial simplifications of the SD equations (32)- (33), discussed in Appendix F, allows us to converge the numerical solution for smaller µ/J than was previously reported [27,48]. This procedure gives access to the conformal ∼ µ 2/3 scaling regime and is also crucial in providing accurate results for the dynamics of the left-right correlators shown in Fig. 2e.
Note that the spectral function A(ω) displayed in Fig. 2d shows intriguing additional structure, beyond what was reported in previous works. We find a sharp peak at ω = E gap followed by an sequence of peaks centered close to harmonics of the gap, with spacing ∆ω ∼ 3E gap . The peak at E gap appears to be infinitely sharp (i.e. resolution-limited in our numerics), while the harmonics get progressively broader as shown in the inset of Fig. 2d. This structure is reflected in the behavior of G ret (t) which shows non-decaying oscillations with a period 2π/E gap at long times, Fig. 2e. The presence of sharp quasiparticle peaks in A(ω) at low frequency suggests an emergent Fermi-liquid description at low energies and temperatures, which is yet to be developed and poses an interesting challenge for future work.  3), (17) and (16), respectively. The asymmetric and regularized forms of the OTOC correspond to the time-ordered correlators in a TFD state given by Eqs. (17) and (10) Taking the latter as a benchmark for finite-N numerical results, we conclude that the regularized OTOC allows to access lower temperatures more reliably than the asymmetric or standard forms.

C. Measuring OTOCs in coupled SYK models
It is known that, in the limit of N → ∞ and at strong coupling βJ 1, the SYK model is maximally chaotic with a Lyapunov exponent saturating the chaos bound λ L = 2πT . In numerical calculations at relatively small N , the maximally-chaotic nature of the SYK model, as seen through the Lyapunov exponents, was never reliably observed and the failure was attributed to finite-N effects [44,52]. Indeed, the exponential growth of the OTOC, parametrized by with A, B real constants, can be expected for times J −1 t < 1/λ L log(N/B) and βJ < N . However, previous numerical calculations were carried out using the standard OTOC in Eq. (3) which shows stronger finite-size effects [37,42]. We compare in Fig. 3 the OTOCs obtained numerically (in a single SYK model) for the three different regularizations discussed in Sec. II B above: standard, regularized and asymmetric. We then extract the Lyapunov exponent for each choice by fitting to the expected functional from, Eq. (37), for intermediate times. Inspired by Ref. [53], we define the lower bound of the fitting region by a time t − such that F (t − ) ∼ 0.98F (0) which marks the beginning of the exponential growth. Similarly, we define the upper bound t + as the time at which the second derivative F (t + ) < 0 and thus cannot describe an exponential. For each regularization we observe an exponential growth characteristic of quantum chaotic systems -however the Lyapunov exponent λ L extracted from our fitting procedure at low temperature differs drastically between the three regularizations. Specifically, the standard and asymmetric forms appear to violate the chaos bound (as also reported elsewhere [44,52]). This is of course not a physical effect, but rather reflects the breakdown of our fitting procedure which occurs because the separation of time scales is insufficient for the small system sizes N considered. The regularized form of the OTOC captures the expected trend for the SYK model (red line in Fig. 3d, cf. Refs. [5,49]) more accurately due to weaker finite-size effects (see also Ref [37]).
As discussed above, the TFD setup naturally leads to regularized OTOCs with a square-root of thermal density matrices inserted inside the trace as indicated in Eq. (16). This is an interesting feature, because such symmetric insertion of thermal factors does not naturally appear in most other measurement schemes such as the Lochsmidt echo or those described in Refs. [11][12][13].
As discussed in Sec. II-D, a possibly more convenient way to access the OTOC is through a measurement of the two-sided Green's function G LR (t−t ) in the ground state of the coupled system. This is clearly a more straightforward measurement, but is limited to weak couplings µ/J. To verify that this approach indeed works we adapt Eq. (23) to the Maldacena-Qi model by identifying O j = χ j . Following the steps outlined in Sec. II-D, we derive the short-time expansion of the retarded version of the averaged LR Majorana propagator, valid for 0 < t µ −1 . Similar to the time-ordered case, G ret LR (t) contains an OTOC contribution (last line of Eq. (38)), and we therefore expect an exponential growth at intermediate times.  39)). This allows to extract λL(T ) in (c) which is consistent with the chaos bound λL = 2πT at low temperature. The scaling of λL(T ) expected [5,49] for the SYK model (as in Fig. 3) is also shown.
In Fig. 4a we show the imaginary part of G ret LR (t) calculated numerically from the large-N saddle-point equations for several values of µ/J. For sufficiently weak couplings µ/J, we can fit an approximately exponential growth in the expected regime, from which we extract a putative Lyapunov exponent λ L (µ) as shown in the inset of Fig. 4a. The extracted exponents follow the ∼ µ 2/3 scaling of the energy gap (see Fig. 2f) at small µ/J. Given that E gap scales linearly with the effective temperature 1/β of the corresponding TFD state [29] (see Fig. 2c), our results imply that λ L ∼ T , consistent with the expectation for the SYK model at low temperatures.
In order to make quantitative statements, we need to establish the coefficient of proportionality of λ L (T ) which requires the knowledge of the function T (µ). This can be in principle obtained from our ED results shown in Fig. 2b. However, because ED does not accurately capture the E gap ∼ µ 2/3 scaling at small µ/J, we do not expect this approach to be quantitatively reliable. On the other hand, as we show below, it is possible to extract the T (µ) dependence directly from the large-N formalism which correctly captures the E gap ∼ µ 2/3 scaling. To do this we apply Eq.
The left-hand side is evaluated in the ground state of the Maldacena-Qi model and gives H L TFD as a function of µ (dropping the SYK superscript from here on). The right hand side is evaluated in the thermal ensemble of a single SYK model and gives H β as a function of temperature. Matching these two energies through Eq. (39) then yields the required function T (µ). The expectation value of the Hamiltonian operator can be extracted from the system Green's functions obtained from the large-N saddle point equations. A textbook procedure [51] applied to the Maldacena-Qi Hamiltonian yields where G αβ (τ ) is the imaginary-time Green's function. Fourier transforming into the Matsubara frequency space and using the spectral representation of G αβ (iω n ) this can be rewritten in the integral form (41) which is convenient for numerical evaluation. Here n(ω) = 1/(e βω + 1) denotes the Fermi-Dirac distribution and ρ αβ (ω) are the spectral functions defined in Appendix F. We use Eq. (41) to evaluate the left-hand side of Eq. (39). The right-hand side can be obtained in an analogous manner and is given by Eq. (41) with µ = 0 and ρ LL replaced by the spectral function ρ(ω) of the single SYK model. The results of this calculation are summarized in Fig. 4b. We find that T = µ for large µ/J while T /J ≈ 0.76(µ/J) 2/3 at small µ/J. Using this result we obtain Lyapunov exponents in good agreement with the prediction for maximal chaos at low temperatures, λ L = 2πT , as shown in Fig. 4c. However, we do not obtain a quantitative agreement with the full solution for the Lyapunov exponents of the SYK model (also shown in Fig. 3). Our method overestimates the value of λ L for intermediate T /J, which could be due to uncertainties in selecting the optimal time window for the fitting procedure, or contamination from the NTOC terms and higher-order terms in µ in the expansion leading to Eq. (38).

IV. PHYSICAL REALIZATIONS, MEASUREMENT SCHEMES
The general scheme to probe quantum chaos using the thermofield double state discussed in Sec. II is applicable to physical systems of essentially any type. The key requirement is to have two identical copies of the system which can be initialized into the TFD state and then measured. Although there are other known systems that exhibit quantum chaos, we continue focusing here on the SYK family of models which are exactly solvable in the large-N limit and have been widely studied in the literature. Below we discuss possible physical realizations of two coupled SYK models, as well as protocols that yield out-of-time ordered correlation functions by performing legitimate causal measurements.
A. Realizations of coupled SYK models

Quantum dots
Perhaps the most conceptually transparent realization of the Maldacena-Qi model [27] is depicted in Fig. 5. It consist of N semiconductor quantum wires proximitized to realize a topological superconductor phase with a pair of Majorana zero modes bound to their ends [54][55][56][57][58]. The wires are weakly coupled to a pair of identical quantum dots, such that the zero modes delocalize into them and form two identical SYK models when inter-actions between the underlying electrons are taken into account [45]. In a wire of finite length L, the two Majorana endmodes are weakly coupled due to the overlap of their exponentially decaying wavefunctions in the bulk of the wire. For the jth wire this coupling has the form iµχ j L χ j R with µ ∼ e −L/ξ cos (k F L), where ξ denotes the superconducting coherence length and k F is the Fermi wavevector of electrons in the wire. Both ξ and k F are sensitive to the gate voltage applied to the wire, which makes the coupling strength µ tunable, at least in principle. A similar term arises for long wires L ξ upon including capacitive effects in each quantum wire. Since the device in Fig. 5 serves as an instructive example below, we expand on some technical details of its realization, following Ref. [45], in Appendix D.
While the device depicted in Fig. 5 may look straightforward, its experimental realization presents a significant challenge for reasons that we now discuss. On the positive side there now exists compelling experimental evidence for Majorana zero modes in individual proximitized InAs and InSb wires. The initial pioneering study by the Delft group [59] has been confirmed and extended by several other groups [60][61][62][63][64][65]. Assembling and controlling large collections of such wires, as would be needed in the implementation of a single SYK model, represents a significant engineering challenge. Constructing an identical pair of SYK models entails another level of difficulty. In the proposal of Ref. [45], random structure of the SYK coupling constants J ijkl originates from microscopic disorder that is present in the quantum dot. It is clearly impossible to create two quantum dots that would have identical configurations of microscopic disorder. A possible solution to this problem would be to engineer quantum dots that are nearly disorder-free, and then introduce strong disorder by hand in a controlled and reproducible fashion. This could be achieved, e.g., by creating a rough boundary or implanting scattering centers in the dot's interior. In such a situation the electron scattering (and therefore the structure of J ijkl ) would be dominated by the artificially introduced defects, and two nearly identical quantum dots could conceivably be produced.

Fu-Kane superconductor
Another proposal to realize the SYK model starts from Majorana zero modes localized in vortices at the interface between a topological insulator (TI) and an ordinary superconductor (SC) [66]. Specifically, when N such vortices are trapped in a hole fabricated in the superconducting layer, and when the chemical potential of the TI surface state is tuned close to the Dirac point, the effective low-energy description of the system is given by the SYK Hamiltonian [44]. The random structure of J ijkl here comes from the randomly shaped hole boundary, and can be well-approximated by a Gaussian distribution in certain limits [44,67]. This setup can be turned into a realization of the Maldacena-Qi model, by using a thin film of a TI with a SC layer equipped with an identical hole on each surface, as illustrated in Fig. 6a. For a thick film this setup generates two decoupled identical copies of the SYK model. For a thin film (e.g. composed of several quintuple layers of Bi 2 Se 3 ), the tails of Majorana wavefunctions extending into the bulk from the two surfaces will begin to overlap. The leading term describing such an overlap will be of the form iµ j χ j L χ j R , as required for the Maldacena-Qi model. The coupling strength µ here will depend exponentially on the film thickness d, but cannot be easily tuned once the device is assembled.
The advantage of this proposal over the quantum dots in Fig. 5 is that randomness in J ijkl here comes from the shape of the hole and is, therefore, under experimental control. Two nearly identical SYK models can conceivably be fabricated in this setup. On the other hand the experimental status of Majorana zero modes in the Fu-Kane superconductor is not nearly as well developed as in quantum wires. Experimental signatures consistent with zero modes bound to individual vortices have been reported in Bi 2 Te 3 /NbSe 2 heterostructures [68,69], but this result remains unconfirmed by other groups. More recently signatures of Majorana zero modes have been observed in surfaces of the iron based superconductor FeTe 0.55 Se 0.45 [70][71][72][73][74]. It is thought that this material is a topological insulator in its normal state, and its surfaces realize the Fu-Kane model when the bulk enters the superconducting phase below the critical temperature T c 14K.
These experimental developments identify the Fu-Kane superconductor as a promising platform for Majorana device engineering. Future efforts might bring us closer to realizing the SYK and Maldacena-Qi models.

Graphene flake bilayers
The complex fermion version of the SYK model, sometimes abbreviated as cSYK, exhibits properties in many ways similar to the canonical SYK model with Majorana fermions [3,6]. It is defined by the Hamiltonian where c j annihilates a complex fermion andμ is the chemical potential. A realization of the cSYK model has been proposed using electrons in the lowest Landau level of a nanoscale graphene flake with an irregular boundary [46]. Once again randomness in J ij;kl originates from the irregular boundary of the flake. Two identical flakes forming a bilayer illustrated in Fig. 6b could realize a complex fermion version of the Maldacena-Qi model if electrons were permitted to tunnel, with weak tunneling amplitude, between the adjacent sites of the two flakes. The tunneling amplitude would depend sensitively on the distance d between the flakes (or on the number of layers in a multi-layer graphene sandwich), but again cannot be easily changed once the device is assembled. This proposed setup eliminates the need for Majorana zero modes, which is a significant potential advantage. On the other hand, the detailed theory of a TFD-like state and its relation to the ground state of the coupled system has not been worked out for the complex fermion version of the model, and we leave this as an interesting problem for future study.

B. OTOC measurement schemes
Because the available experimental probes will depend on the specific details of the physical realization, we offer here only general remarks on how OTOC may be measured using the protocols developed in Sec. II. For concreteness and simplicity, we focus again on the proposed coupled SYK dot realization of the Maldacena-Qi model, depicted in Fig. 5, but we expect our discussion to be valid more generally.
At the highest level, we may distinguish two types of situations when attempting to probe OTOC through a causal (time-ordered) measurement: we either have the ability to control the coupling strength µ on microscopic timescales (i.e. times of order /J), or we do not. In the first case we can manipulate µ to prepare the initial resource state |TFD β (−t) , and then perform a two-sided time-ordered measurement as discussed in Sec. II and Appendix A. This has the advantage of directly probing the regularized or asymmetric OTOCs. We give some concrete examples of this below. If µ cannot be controlled on microscopic timescales, it is still possible to extract the OTOC by measuring G ret LR (t) in a system with constant nonzero µ, as discussed in Sec. II-D. While this measurement is in principle easier, the quantitative interpretation is less clean, because it necessitates disentangling of the OTOC contributions from the NTOC terms in Eq. (23).

When µ can be controlled
Following existing theoretical work on eternal traversable wormholes [27,29], we discussed a method to reliably create a TFD state by cooling down to the ground state of a weakly coupled two-system Hamiltonian H(λ) = H 0 + λH I . Measuring the OTOC requires the TFD state evolved to negative time, |TFD β (−t) , which we demonstrate in Appendix A can be achieved, for short time durations at least, by tuning the dimensionless coupling λ. We emphasize that in a generic many-body system, this should be a much easier task than true backward time evolution of a many-body excited state that would normally be required to measure an OTOC. Such backward time evolution necessitates the reversal of the sign of the many-body Hamiltonian H 0 which is, in the vast majority of cases, not feasible by any known technique. On the other hand, manipulating the strength of couplings between two systems can often be achieved, e.g., by gating, as discussed in the previous Section and Appendix D for the setup of Fig. 5.
With the above caveats, the proposed protocol to measure OTOC using the TFD state could be defined follows. (i) Prepare two identical copies of the system that are weakly coupled and described by H(1) = H S . (ii) Cool the coupled system to a physical temperature T phys that is much smaller than the energy gap of the combined system, which puts it into its ground state well approximated by |TFD β (0) . (iii) Increase coupling λ to a value larger than 1 for a short period of time. This creates a good approximation of |TFD β (−t) where the time evolution is with respect to H 0 , cf. Appendix A. (iv) Decouple the system by setting λ = 0, and probe it by a conventional measurement. One possibility, mathematically expressed in Eq. (17), is to excite the system on one side at time −t and then preform a two-sided measurement at time zero. This procedure yields a direct measure of the asymmetrically-regularized OTOC shown in Fig. 3.
A way to measure two-sided two-body Majorana operators is via wire-charge measurements, cf. Refs. [75,76] and detailed in Appendix D, that only rely on a capacitive coupling between an external readout circuit and the nanowire charges. We consider the simplest case where a single, collective gate in Fig. 5 couples to all nanowires. Quantizing a fluctuating global gate voltage v g (t) → [a(t) + a † (t)], and assuming roughly isotropic capacitive coupling parameters ∼ g j g to all nanowires, one finds with a single photon species a and total nanowire-charge operatorQ(t) = jq j (t) = j iχ j L (t)χ j R (t). The term H res encodes the external resonator readout circuit, and generates the dynamics for resonator photons a(t). Either the transmission amplitude or -phase shifts in this external circuit, by means of the capacitive coupling in Eq. (43), then yield a probe of the total nanowire charge Q(t). The latter is directly related to the averaged equaltime left-right Green's function, as G LR (t) = −iQ(t)/N .

When µ is fixed
In Secs. II-D and III-D, we showed that the averaged LR Green's function of the coupled theory at fixed µ contains information on the OTOC for small couplings µ/J and short times. We now argue that the retarded version of the LR Green's function [Eq. (38)] can be probed by a straightforward tunneling measurement in the setup of Fig. 5. Consider a tunnel probe weakly coupled to one of the wires at a point distance x from its left end (represented as the central probe in Fig. 5). A standard tunneling experiment measures differential tunneling conductance g(V ) = dI/dV , which is proportional to the electron spectral function in the wire ρ x (ω) at point x and frequency ω = eV , where V is the applied bias voltage. In Appendix C we show that this quantity is related to the retarded LR two-point Majorana correlator by a simple relation The constant of proportionality K x depends on the tunneling matrix element and on the position x in the wire, but is time independent as the measurement is performed under equilibrium conditions. Therefore, time dependence of G ret LR (t) and the relevant Lyapunov exponent can be extracted from the measured spectral function, using Eq. (44).
The result given in Eq. (44) relies on two simple observations, discussed in more detail in Appendix C. First, a retarded time-domain correlator of Hermitian operators is purely imaginary. This fact follows directly from its definition and underlies the proportionality of iG ret LR (t) to a real quantity. Second, Majorana zero mode operators χ k L/R are simply related to the electron operators in the wire through the solution of the relevant Bogoliubovde Gennes equation [54][55][56][57][58], which is largely dictated by symmetries of the setup. This implies proportionality of iG ret LR (t) to the electron spectral function, Fouriertransformed into the time domain.
Based on these observations we expect Eq. (44) to be robust and independent of system details. Remarkably, in conjunction with Eq. (38) it connects the Lyapunov exponent of an OTOC with the electron spectral function in a proximitized semiconductor system, which is routinely measured in tunneling and other spectroscopic experiments.

V. CONCLUSIONS AND OUTLOOK
In this work, we introduced and extensively tested the concept of entanglement in the thermofield double state as a tool to measure out-of-time ordered correlators in quantum many-body systems. OTOCs have been of great interest recently because their exponential growth at intermediate times provides direct access to diagnosing quantum-chaotic behavior in many-body systems.
While previous work has implemented OTOC measurements in small-scale and highly controllable quantum systems [10][11][12][13], these approaches do not lend themselves to the analysis of large, complex many-body systems that realize quantum chaos in solid-state platforms. Based on the thermofield double state, one of the main workhorses for the theoretical description of black hole and wormhole quantum physics [7,25,26], we proposed and tested new protocols for OTOC measurement, where the preparation of a specific resource state -namely the TFDreplaces the need for complicated time-evolution or echo procedures at a later stage. The TFD entangled pair can be obtained as a unique ground state of two coupled copies of the interacting quantum system under investigation [27][28][29]. We showed that a conventional measurement with no or only minimal control of the system parameters can directly access the so-called regularized OTOCs. The latter have been introduced as mathematical objects in field-theoretical calculations, because in certain limits they are less singular than the canonical OTOCs. Regularized OTOCs have recently been argued to measure quantum chaos more reliably than canonical OTOCs [35][36][37], a result corroborated by our numerical analysis. However, regularized OTOCs are even more difficult to access than canonical OTOCs in a physical measurement, given that the insertion of square roots of the density matrix on their Schwinger-Keldysh contours (see Fig. 1) does not reflect a sensible thermal measurement, even if backward time evolution is considered possible. To our knowledge, only the interferometric approach of Ref. [17] potentially allows for their extraction. It is all the more exciting that they arise as naturally accessible objects in our TFD-based protocols.
Perhaps the most surprising outcome of our considerations is the realization, expressed mathematically in Eqs. (23) and (25), that the Lyapunov exponent λ L characterizing quantum chaos is in fact encoded in the intermediate-time behavior of the ordinary two-point correlator G LR (t). The latter is measured under equilibrium conditions, between operators drawn from the two subsystems forming the TFD. We confirmed this result through a numerical solution of the large-N saddle point equations associated with two coupled SYK models. These indeed show approximate exponential growth of G LR (t) at intermediate times with a Lyapunov exponent λ L ≈ 2πT consistent with the presence of maximal chaos, saturating the Maldacena-Shenker-Stanford bound [7] λ L ≤ 2πT in the weak coupling limit µ/J 1. This finding is significant because in many systems such two-point correlators can be probed without much difficulty by spectroscopic techniques. For example, in the proposed quantum dot realization of two coupled SYK models illustrated in Fig. 5, the retarded Majorana correlator G ret LR (t) is found to be proportional to the Fourier transform of the electron spectral function ρ x (ω) (see Eq. (44)) which is accessible through a routine tunneling measurement.
Last, in this work we have made substantial progress in understanding the structure of the large-N Schwinger-Dyson equations for the Maldacena-Qi model [27] comprised of two coupled SYK models, cf. Sec. III-A and Appendix F. We showed that it becomes possible to describe the full real-time dynamics in terms of a single (retarded) Green's function, the corresponding spectral function and a single self-energy. Finding an explicit analytical solution for the dynamics of such coupled quantum chaotic systems, at least in certain limiting cases, would clearly be very rewarding. Specifically, as we argued, it should be possible to extract the intermediate-time exponential growth of G LR (t) and the corresponding chaos exponent directly from the large-N saddle point equations. By contrast, in the single SYK model one has to go beyond the saddle point equations and sum an infinite series of ladder diagrams to evaluate the OTOC [4,5].
As an outlook, interesting future work includes the detailed investigation of physical platforms for coupled chaotic quantum systems, for example, based on the ideas presented in Sec. IV-A and Refs. [44][45][46]. The key challenge here will be to prepare two systems that are nearly identical in that the interaction coupling constants J ijkl are essentially the same in both. We discussed several possible approaches to this challenge in Sec. IV-A, but neither is fully satisfactory. Going forward, the most promising route appears to involve an exact microscopic symmetry that would relate two subsystems. For instance time-reversal Θ in a systems of spin-1 2 fermions would mandate identical J ijkl (up to a complex conjugation) for the two spin projections. A closely related attractive research direction is towards realizations of "wormholes" in coupled complex-fermion SYK models [3,6]. Use of complex fermions would alleviate the need for Majorana zero-modes as basic ingredients and would reintroduce electron spin as a potentially useful degree of freedom. While it is not obvious at present how to formulate the corresponding complex-fermion TFD state, guidance can be taken from the pedagogical discussion of Ref. [29].
Further, the generality of the construction in Ref. [29] suggests that many more interesting physical systems, including in higher spatial dimensions, might lend themselves to an investigation of their chaotic behavior using our proposed method. We hope that our findings will stimulate further developments on both theoretical and experimental fronts, which will eventually lead to practical tools for quantum chaos diagnosis in interacting many-body systems.
Using Eq. (A2), the last exponential evaluates to e −iλ 0t . For short time durations one can furthermore neglect the t 2 term in the first exponential, which leads to We see that increasing the coupling strength λ to a value larger than 1 has the same effect as evolving the state |Ψ 0 backward in time under the decoupled Hamiltonian H 0 = H L + H R . Manipulating the coupling strength λ can therefore be used to prepare the TFD state at negative times, and represents a simple alternative to engineering a sign inversion of the complicated, interacting Hamiltonian H 0 . Eq. (A6) remains true for sufficiently short times t such that one can neglect the commutator term in the exponential of Eq. (A5). This can be estimated from the condition Some intuitive understanding of the backward time evolution described above can be gained by analyzing an example of a simple system. Consider a spin-1/2 degree of freedom in a magnetic field described by the Hamiltonian H 0 = −Jσ z and H I = −µσ x . The ground state |Ψ 0 of the combined Hamiltonian H(λ) = H 0 +λH I with λ = 1 has the spin pointing along the direction parallel to the total magnetic field B = (µ, 0, J). If we switch off H I the spin will start precessing counterclockwise around the field direction B 0 = (0, 0, J) associated with H 0 . This is analogous to the TFD state evolving according to the decoupled Hamiltonian H 0 forward in time. On the other hand, if we instead increase the value of λ, the spin will start precessing counterclockwise around the new field direction B λ = (λµ, 0, J), as illustrated in Fig. 7a. At short times t we observe that the evolution for λ > 1 approximates backward time evolution under H 0 . This effect can be quantified by calculating the overlap between the ground state evolved backward in time according to H 0 , i.e. |Ψ 0 (−t) = e +iH0t |Ψ 0 and the ground state evolved forward in time according to H(λ), |Ψ λ (t) = e −iH(λ)t |Ψ 0 . Elementary but somewhat tedious calculation gives an explicit expression for P (t) that we plot in Fig. 7b, for several values of µ. We observe that for short times and µ J the overlap remains very close to 1, confirming that the method indeed yields an excellent approximation to the state |Ψ 0 evolved backward in time. Crucially, this backward time evolution does not require reversing the sign of H 0 , and is achieved solely by controlling the strength of the H I perturbation.
In Fig. 7c we present numerical evidence supporting this claim for the coupled Maldacena-Qi model. Similar behavior as described above is observed at short times in this interacting many-body system suggesting that it is generic and thus can be used to prepare the required initial state |Ψ 0 (−t) in a wide variety of settings.

Appendix B: Short-time expansion of LR two-point correlator
In this Appendix we derive Eq. (23), which we used in the main text to argue that the LR two-point correlator G LR (t, t ), defined in Eq. (22), contains at short times information on the OTOC. We proceed by evaluating the two-point correlator Here |Ψ 0 denotes the ground state of the combined system, described by the Hamiltonian H = H L + H R + H I , which we will approximate later on by |TFD β . We work in the Heisenberg picture where V α (t) is an arbitrary Hermitian operator evolving according to the full Hamiltonian H. G LR (t, t ) is a naturally time-ordered correlator that a physical probe would measure in the ground state of the combined system. We would like to know how this quantity is related to what a physical probe would measure in a thermal ensemble at inverse temperature β of a single, decoupled system. Mathematically, the goal is to express G(t, t ) as an average with respect to the TFD state of operators that evolve according to H 0 = H L + H R . To proceed, we pass to the interaction picture by writing [51] where superscript I denotes the interaction picture and U (t, t ) = e iH0t e −iH(t−t ) e −iH0t is the unitary operator that translates between the Heisenberg and interaction pictures. Using Eq. (B2), the correlator becomes where we used the property U (t, s)U (s, t ) = U (t, t ). For simplicity, we henceforth also assume that t > t .
We now employ a standard result of diagrammatic many-body theory [51] that express U (t, t ) as a series expansion in H I (t) of the form where . . . 0 denotes the expectation value with respect to the ground state |Ψ 0 . This expression is valid when one can neglect all higher order terms in the expansion (B4) of U (t, t ), U (0, t) and U (t , 0). This requires short time durations |t − t | as well as individually small |t| and |t |. In the following we focus on the symmetric case t = −t, which has a convenient property that short duration |t − t | automatically assures that |t| and |t | are small.
As the final step we substitute for H I the Maldacena-Qi form given in Eq. (20), and approximate |Ψ 0 by |TFD β . With these choices, the expectation values in Eq. (B5) can be expressed in terms of single-sided averages using the procedure explained in Sec. II-B of the main text, see especially the steps leading to Eq. (15). The correlator thus becomes where η = +/− for bosonic/fermionic operators. We dropped superscript I and subscripts R/L on all operators; it is understood that they now evolve according to the single-sided Hamiltonian, say H L , while the traces are taken with respect to the eigenstates |n of the same Hamiltonian. Each individual trace in Eq. (B6) is timetranslation invariant, and in the following we find it convenient to shift all temporal arguments of operators inside the traces by +t. In addition, using the cyclic property of the trace, we may combine the second and third lines and express the last line as an integral from 0 to t. This leads to the form quoted in the main text Eq. (23). Mathematically, Eq. (B6) can be viewed as an expansion of the propagator G in powers of a dimensionless time variablet = µt to first order. Higher order contributions that would result from the omitted terms in Eq. (B4) can be neglected when µt 1 which constrains the expected domain of validity of Eq. (B6) to short times or small values of coupling µ.

Appendix C: OTOC from electron spectral function
In this Appendix we derive Eq. (44), which provides a simple route to access OTOC and Lyapunov exponent through an equilibrium tunneling measurement of the electron spectral function in a wire that forms a part of the device shown in Fig. 5.
As a first step, we show that a time-domain retarded propagator of Hermitian operators is imaginary valued. Consider a retarded propagator defined as evaluated in the ground state |Ψ 0 of the system. Expanding the anticommutator and using the basic property of the inner product a|b = b|a * , we can rewrite the average as A(t)B(t ) 0 + [B(t )A(t)] † * 0 . For Hermitian operators A and B this equals 2Re A(t)B(t ) 0 , and G ret (t, t ) is therefore purely imaginary.
As we already mentioned, a tunneling measurement can be used to extract the electron spectral function ρ x (ω) which is related to the electron propagator G x (ω). We therefore start by considering the corresponding quantity defined in the time domain as where c x (t) annihilates electron at time t and at point x of the wire. We want to relate G(t) to the Majorana propagators G αβ (t) discussed in the main text. To this end recall the relations which follow from the solution of the Bogoliubov-de Gennes (BdG) equations for the Majorana wire [54][55][56][57][58].
Here Φ L/R (x) represent the Majorana wavefunctions. They are real-valued and peaked at the L or R end of the wire, respectively, with exponentially decaying tails extending into the wire. The form of Eqs. (C3) is constrained by the choice that χ L/R transform as even/odd under time reversal. We may invert Eqs. (C3) to express the electron operators in terms of the Majoranas, where the dots represent the remaining quasiparticle operators that form the complete set of solutions of the BdG equations. We will assume that these occur at non-zero energies separated from the zero-mode manifold by a gap, and will thus not affect the low-energy spectral function. Substituting into Eq. (C2) and neglecting these terms, we find the Majorana fermion contribution to the electron Green's function Making a further non-essential assumption of mirror symmetry between the L and R sides of the system, we have G LL (t) = G RR (t) and G LR (t) = −G RL (t). Also, we note that the same calculation can be repeated for the retarded function iG ret To complete this calculation, it remains to relate Re[G ret x (t)] to the electron spectral function which is the observable quantity. We use the spectral representation G ret x (ω) = ∞ −∞ dω ρ x (ω )/(ω − ω + iδ), which, upon Fourier transforming, gives Taking the real part and recalling that ρ x (ω) is strictly real, we obtain Finally combining with Eq. (C7), we arrive at Eq. (44) of the main text.
Appendix D: Majorana nanowire device One of our example realizations of two coupled chaotic systems is the two-sided SYK device shown in Fig. 5, inspired by the SYK setup of Chew, Essin and Alicea [45]. Since this is the most transparent and directly controllable realization we discuss, let us here expand on its technical underpinning in the framework of Ref. [45]. The system in Fig. 5 is described by N sets of Majorana modes χ j L , χ j R at the left and right ends of the nanowires, and by the N L/R N complex fermions c s,α=L/R hosted in disordered wave-functions of its quantum dots. We write the dot fermions in Majorana representation c s,α = (η s,α +iη s,α )/2, where η s,α is even under time-reversal (TR) whileη s,α is odd. Similarly we take nanowire Majoranas χ j L and χ j R to be even and odd under TR, respectively. Assuming that both quantum dots preserve the BDI symmetry class of the Majorana sector [45], the left and right ends of the device then are guaranteed to host N TR-even (TR-odd) Majorana zero modes. This statement holds unless the LR couplings ∼ iχ j L χ j R are introduced, where in the main text we discuss the full crossover from weak to strong bilinear couplings.
To express the toy-model Hamiltonian describing the device in Fig. 5, it is now convenient to introduce Majorana spinors χ L/R , η L/R and η L/R for the respective left/right groups of Majorana fermions. Following Ref. [45], we first describe the hybridization of Majoranas χ L/R into the left and right quantum dots as (D1) Here M λα are real rectangular matrices of couplings ∼ λ between the wire and dot Majoranas, cf. Fig. 5, and M α are real N α × N α matrices encoding the dots level structure ∼ . The coupling Hamiltonians H 0,α=L/R can be diagonalized by orthogonal rotations obtained from a singular-value decomposition of the coupling matrix encoded in M λα and M α [45], giving new Majorana spinors and η L = Oη ,L η L , and similar for right Majorana operators with L → R and η ↔ η. In these new operators, where η j,α andη j,α correspond to the rotated spinors η α and η α , and E j,α are associated hybridization energies obtained from rotating the coupling matrix. Note that as guaranteed by the BDI classification, each dots hosts N zero-modes encoded in spinors χ α that do not appear in H 0,α . Further, assuming strong wire-dot hybridizations λ N δ , where δ is the typical level spacing in the dot, all wire Majoranas χ α are absorbed into the respective left and right dots. Finite-energy modes in Eq. (D3) are split off to energies ∼ N δ by level repulsion. We now add the LR coupling between the two SYKdots. In the original basis, sensible couplings are pairwise between Majoranas χ j L and χ j R on each wire, yield- The specific form and possible tuning of couplings µ j is discussed below. Inverting the orthogonal transformations in Eq. (D2), one can represent the original wire Majoranas in the new variables as In the absence of left-right symmetry of the underlying disordered quantum dots in Fig. 5, there is no reason to assume that µ χχ is diagonal. Nevertheless, it is always possible to re-diagonalize the LR Hamiltonian by re-applying the orthogonal transformations to Majoranas χ L/R . Given that the original coupling matrix µ LR was diagonal, one finds H int,eff = i χ T L µ LR χ R with and similarly χ R = O T χχ,R χ R . Note O T χχ,s O χχ,s = 1, since only the full transformation in (D2) is orthogonal. We hence observe that under the approximations taken above, the interaction matrix µ LR = diag(µ j ) remains unchanged by the absorption of Majorana modes χ j,L/R into the SYK dots, but now acts on the new modes χ L/R .
Finally, for ease of notation, we relabel χ L/R → χ L/R and proceed to add on-site interactions for the remaining Majorana zero-modes, inherited from Coulomb interactions of the dot fermions c s,α=L/R [45]. Generically one then obtains a Hamiltonian H = H L + H R + H int,eff with site-dependent four-Majorana interactions as in Eq. (27) H SYK α=L/R = i<j<k<l J α ijkl χ i s χ j s χ k s χ l s .
From here, assuming left-right symmetry of the SYKdots in Fig. 5 such that J L ijkl = J R ijkl identically, one obtains the Maldacena-Qi Hamiltonian [27] in Sec. III.
We now discuss how one may realize (tunable) one-toone bilinear couplings across nanowires of the device in Fig. 5. While one may use residual Majorana hybridizations µ j ∼ µ 0 e −L/ξ that decay exponentially with wire length L, in practice it is desirable to tune, or at least switch on and off, the couplings in a more controllable fashion. To this end, consider the TS nanowires in Fig. 5 to be strongly coupled to a ground bulk superconductor, but not fully grounded. For a single wire, both its intrinsic single-electron charging energy E c and Josephson coupling E J to the ground bulk SC then become relevant. In the limit of large but finite E J /E c > 1, one finds an effective parity splitting between even and odd charge states on the wire, directly translating to a Majorana parity splitting for pairs χ j L , χ j R [77]. On the Hamiltonian level, this term can be incorporated as H int,j = µ(n j )iχ j L χ j R , µ(n j ) = µ 0 cos(πn j ) , (D8) where µ 0 depends on both E c and E J [77], and n j is a gate parameter set by a nearby electrostatic gate, thereby controlling the equilibrium charge (parity) on the wire. We thus see that a nearby collective gate as in Fig. 5, controlling charge on all nanowires in the device, can simultaneously switch on and off the bilinear Majoranacoupling across all pairs of modes χ j L/R . Further, at least in principle, one can also address (few) wires individually via additional gates not shown in Fig. 5.
Appendix E: Measurements in the nanowire device While one can measure some properties of the coupledwire SYK device with simple tunnel-probes as indicated in Fig. 5, here we mention another useful capability that comes with the inter-side coupling implementation in Appendix D. For each individually tunable gate voltage in the device, e.g. the collective gate or ones attached to the top/bottom-most nanowires, one can perform projective readouts of the nanowire paritiesq j = iχ j L χ j R or certain combinations thereof. The more individually addressable gates, the more completely one may map out the space of parity eigenvalues q j=1,...,N . For a detailed discussion of Majorana-parity readout via resonators attached to electrostatic gates, see e.g. Ref. [75,76].
To illustrate how this readout works, assume that a gate parameter n in Eq. (D8) is set such that the corresponding interaction is nearly switched off, n(t) = 1 2 + v(t) with a small fluctuating gate voltage |v(t)| 1. We then introduce a resonator circuit capacitively coupled to the gate, described by resonator photons a(t). Upon quantizing the fluctuating gate voltage v(t), assuming a capacitive interaction strength g between resonator gate and nanowire, one replaces v(t) → [a(t)+a † (t)]. The total nanowire-charge readout setup then is described by whereq(t) is the time-evolving nanowire charge parity. H res here encodes the resonator spectrum, and generates the dynamics for resonator photons a(t).
In the strong-coupling regime, leading to a net exchange of resonator photons with the system, one can directly access e.g. the transmission amplitudes or phase shifts of the resonator that depend on a(t) , and via Eq. (E1) also on q(t) [75]. This readout mode hence allows for a direct, time-resolved measurement of Majorana parities q(t) between the two coupled SYK dots in Fig. 5.
Here b 0 (t, t ) = b 0 (t − t ) is a bare photon Green's function of the uncoupled resonator Hamiltonian H res . Note that a time-dependent Green's functions as in Eq. (E2) is encoded by time-and frequency-resolved resonator occupations B a (τ, ω) = dτ e iωτ B a (τ + τ 2 , τ − τ 2 ), that can be measured in principle. Ignoring back-action of the resonator on the system (that generates an effective interaction as in Eq. (D8)), and to lowest order in chargeresonator interaction g, this gives information about d 0 (t, t ) = q(t)q(t ) 0 . The latter expression constitutes a four-Majorana two-sided correlator evaluated with respect to the bare SYK Hamiltonian H = H SYK L + H SYK R . In Sec. IV B of the main text, we discuss how charge or charge-correlator measurements can become useful tools to access OTOCs in coupled-wire SYK dots.
To find the fixed-point solution of the SD equations we hence perform numerical iteration starting from the ansatz for spectral functions given in Eqs. (F10) with small nonzero δ. The self-energy Σ ret + (ω) is calculated using Eqs. (F16) and (F17) with help of fast Fourier transform algorithms. G ret + is then computed from Eq. (F15) and used to reconstruct the full retarded propagator. New spectral functions are then extracted from Eqs. (F8) and the procedure is restarted from these. To improve the convergence properties we follow Ref. [5] and after each round of iteration we mix the initial propagator with the new solution obtained from the SD equation. We declare convergence to a physical solution once the propagators, spectral densities and, in particular, the energy gap stop changing within the specified accuracy. We also check that the solutions are stable with repect to increasing the number of iterations, frequency/time resolution and cutoffs, and other non-physical ingredients of the numerical solution, such as the initial broadening δ. Results of our numerics are discussed in Sec. III of the main text.