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

TL;DR

This paper introduces a novel protocol for diagnosing quantum chaos in many-body systems using entanglement and thermofield double states, avoiding the need for backward time evolution and enabling experimental measurement of chaos indicators.

Contribution

It proposes a new method to measure out-of-time-ordered correlators via ordinary measurements on thermofield double states, facilitating experimental detection of quantum chaos.

Findings

Lyapunov exponent can be extracted from two-point correlations in TFD states.

The method yields the regularized OTOC, a key indicator of quantum chaos.

Numerical simulations on the Sachdev-Ye-Kitaev model support the protocol's effectiveness.

Abstract

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.

Equations

$$C_{\eta }(t)=-⟨[W(t),V(0){]}_{\eta }^2⟩$$

$$C_{\eta }(t)$$

$$F(t)=⟨W(t)VW(t)V⟩,$$

$$C_{\eta }(t)=2\eta -2\mathrm{Re}F(t).$$

$$|{\mathrm{TFD}}_{\beta }⟩=\frac{1}{\sqrt{Z_{\beta }}}\sum _n{e}^{-\beta E_n/2}|\bar{n}{⟩}_L\otimes |n{⟩}_R$$

$$⟨{\mathcal{O}}_L{⟩}_{\mathrm{TFD}}=Z_{\beta }^{-1}\sum _n{e}^{-\beta E_n}{}_L⟨n|{\mathcal{O}}_L|n{⟩}_L.$$

$$(H_L-H_R)|{\mathrm{TFD}}_{\beta }⟩=0.$$

$$e^{-it(H_L-H_R)}|{\mathrm{TFD}}_{\beta }⟩=|{\mathrm{TFD}}_{\beta }⟩.$$

$$\mathcal{F}(t_1,t_2)$$

$$\tilde{F}(t,t^{\prime })=⟨\mathcal{T}[V_L(t)W_R(t)V_R(t^{\prime })W_L(t^{\prime })]{⟩}_{\mathrm{TFD}}$$

$$\begin{gathered} \tilde{F}(t,t^{\prime })=\frac{1}{Z_{\beta }}\sum _{n,m}{e}^{-\beta (E_n+E_m)/2}{}_L⟨\bar{n}|V_L(t)W_L(t^{\prime })|\bar{m}{⟩}_L \\ \times {}_R⟨n|W_R(t)V_R(t^{\prime })|m{⟩}_R, \end{gathered}$$

$$\begin{gathered} \tilde{F}(t,t^{\prime })=\frac{1}{Z_{\beta }}\sum _{n,m}{e}^{-\beta (E_n+E_m)/2}⟨m|W(-t^{\prime })V(-t)|n⟩ \\ \times ⟨n|W(t)V(t^{\prime })|m⟩, \end{gathered}$$

$$⟨\bar{n}|\mathcal{O}(t)|\bar{m}⟩=⟨m|\mathcal{O}(-t{)}^{†}|n⟩.$$

$$y^4=e^{-\beta H}/Z_{\beta }.$$

$$\tilde{F}(t,t^{\prime })=\mathrm{tr}[W(-t^{\prime })V(-t)y^{2}W(t)V(t^{\prime })y^2].$$

$$\tilde{F}(t,-t)=\mathrm{tr}[W(2t)V(0)y^{2}W(2t)V(0)y^2],$$

$$\begin{gathered} ⟨W_L(-t)V_L(0)V_R(0)W_L(-t){⟩}_{\mathrm{TFD}}= \\ =\mathrm{tr}[W(t)V(0)W(t)y^{2}V(0)y^2]. \end{gathered}$$

$$H_S=H_L+H_R+H_I,$$

$$H_I=\sum _j{c}_j{d}_j^{†}{d}_j,d_j={\mathcal{O}}_L^j-\Theta {\left({\mathcal{O}}_R^j\right)}^{†}{\Theta }^{-1}$$

$$H_I=i\mu \sum _j{\mathcal{O}}_L^j{\mathcal{O}}_R^j$$

$$H_0=H_L+H_R.$$

$$i{G}_{LR}(t,t^{\prime })=⟨\mathcal{T}{V}_L(t)V_R(t^{\prime })⟩.$$

$$i{G}_{LR}(t,-t)$$

$$F(t_1,t_2)\simeq b(t_1-t_2)e^{{\lambda }_L(t_1+t_2)/2}$$

$$i{G}_{LR}(t,-t)\simeq A+B{e}^{{\lambda }_{L}t},$$

$$H=H_L^{\mathrm{SYK}}+H_R^{\mathrm{SYK}}+i\mu \sum _j{\chi }_L^j{\chi }_R^j$$

$$H_{\alpha }^{\mathrm{SYK}}=\sum _{i<j<k<l}{J}_{ijkl}{\chi }_{\alpha }^i{\chi }_{\alpha }^j{\chi }_{\alpha }^k{\chi }_{\alpha }^l,$$

$$\overline{J_{ijkl}}=0,\overline{J_{ijkl}^2}=\frac{3!}{N^3}{J}^2$$

$$c_j=\frac{1}{\sqrt{2}}({\chi }_L^j-i{\chi }_R^j)$$

$$P=(-i{)}^N\prod _{j=1}^N{\chi }_L^j{\chi }_R^j\cong \left(N_f\right)\text{mod}2,$$

Full text

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 Hawking (1976). Black holes are believed to scramble at the fastest possible rate consistent with causality and unitarity Hayden and Preskill (2007). Some strongly coupled quantum systems, such as the Sachdev-Ye-Kitaev (SYK) model Sachdev and Ye (1993); Kitaev (2015); Maldacena and Stanford (2016); Sachdev (2015) 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 $\lambda_{L}$. A universal upper bound on chaos, conjectured by Maldacena, Shenker and Stanford Maldacena et al. (2016), posits that $\lambda_{L}\leq 2\pi 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 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 Swingle et al. (2016); Zhu et al. (2016); Swingle (2018). So far this has been achieved in a very limited range of systems, mainly quantum simulators Li et al. (2017) and ion traps Gärttner et al. (2017); Landsman et al. (2019). 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 $\lambda_{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. Daley et al. (2012); Abanin and Demler (2012); Islam et al. (2015) 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 Yao et al. (2016) or out-of-equilibrium measurement protocols Campisi and Goold (2017); Yunger Halpern (2017). 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 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 Shenker and Stanford (2014),

[TABLE]

of two initially commuting Hermitian operators $[W,V]_{\eta}=0$. In the following we allow for both bosonic ($\eta=1$) and fermionic ($\eta=-1$) statistics, with $[\>\cdot\>,\>\cdot\>]_{\eta}$ denoting the commutator (anti-commutator) for $\eta=1$ ($\eta=-1$). The operators evolve in time according to the system Hamiltonian $H$ through $W(t)=e^{iHt}We^{-iHt}$ and $\langle...\rangle$ denotes the thermal average at inverse temperature $\beta=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_{\eta}(t)$ to grow as a function of time and eventually saturate at a value close to $2\langle V^{2}\rangle\langle W^{2}\rangle$ for large $t$, regardless of the specific form of $V$ and $W$. In chaotic many-body systems, at intermediate times, the growth of $C_{\eta}(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

[TABLE]

The averages on the first line represent naturally time-ordered correlators (NTOC) that correspond to “sensible” experiments performed in an ordinary quantum system. For instance the second term $\langle VW(t)W(t)V\rangle$ 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. (I) represent out-ot-time-ordered correlators. These correspond to less sensible experiments. The second term for instance, which is often denoted as

[TABLE]

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_{\eta}(t)$ can be expressed simply as

[TABLE]

The fundamental difference between NTOC and OTOC is best visualized by placing the operators on the Schwinger-Keldysh contour illustrated in Fig. 1. The imaginary time evolution is generated by the density matrix which becomes explicit if we rewrite Eq. (3) as $F(t)={\rm tr}[W(t)VW(t)Ve^{-\beta 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 Hahn (1950); Peres (1984); Jalabert and Pastawski (2001); Gorin et al. (2006), 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 Gao et al. (2017); Maldacena et al. (2017); Maldacena and Qi (2018); Gao and Liu (2019). 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 protocols. Importantly, recent theoretical work has established a simple method that can be used to prepare the TFD state Cottrell et al. (2019). 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 Maldacena and Qi (2018) (see also Fig. 2a). We also note interesting recent works on how to prepare a TFD state using quantum circuits Wu and Hsieh (2018); Martyn and Swingle (2019); Zhu et al. (2019); however, these approaches are limited to the moderate system sizes accessible on present-day quantum simulators, similar to the experiments of Refs. Li et al. (2017); Gärttner et al. (2017); Landsman et al. (2019).

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 $\lambda_{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.

II.1 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_{\alpha}$ ($\alpha=L,R$) are invariant under time reversal generated by an antiunitary operator $\Theta$. The TFD state at inverse temperature $\beta$ is then defined as

[TABLE]

where $|n\rangle_{\alpha}$ is an eigenstate of $H_{\alpha}$ with energy eigenvalue $E_{n}$, and $Z_{\beta}=\sum_{n}e^{-\beta E_{n}}$ is the partition function. $|\bar{n}\rangle=\Theta|n\rangle$ denotes the time-reversed partner of the eigenstate $|n\rangle$ 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\rangle$ is defined up to an overall phase $e^{i\phi_{n}}$. A direct product $|\bar{n}\rangle_{L}\otimes|n\rangle_{R}$ is however well-defined because $|\bar{n}\rangle$ transforms with the opposite phase to $|n\rangle$ under the corresponding $U(1)$ transformation. In the limit of zero temperature $|{\rm TFD}_{\beta}\rangle$ 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 $|{\rm TFD}_{\beta}\rangle$ is given by a thermal average,

[TABLE]

It is also important to note that $|{\rm TFD}_{\beta}\rangle$ 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

[TABLE]

This has implications for the concept of time-translation invariance in the TFD state. Eq. (7) implies that $|{\rm TFD}_{\beta}\rangle$ evolves trivially under $H_{-}$,

[TABLE]

Hence the expectation value of a product of two operators acting in the $L$ and $R$ systems has the property

[TABLE]

valid for arbitrary $t$. The second line follows upon replacing $|{\rm TFD}_{\beta}\rangle$ in the expectation value on the first line by $e^{-it(H_{L}-H_{R})}|{\rm TFD}_{\beta}\rangle$ using Eq. (8), and recalling that ${\cal O}_{R}$ commutes with $H_{L}$ (and same for $L\leftrightarrow R$). Choosing $t=t_{2}$ we see that ${\cal F}(t_{1},t_{2})$ is a function of $t_{1}+t_{2}$ only. This should be compared to the statement of time-translation invariance in a conventional (unentangled) state, $\langle{\cal O}(t_{1}){\cal O}(t_{2})\rangle=\langle{\cal O}(t_{1}-t){\cal O}(t_{2}-t)\rangle$, 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.

II.2 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

[TABLE]

evaluated with respect to the TFD state in Eq. (5). Here ${\cal T}$ denotes the time-ordering operator. Normal time-ordered 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^{\prime}$ and write the average in Eq. (10) explicitly using the TFD state (5). We obtain

[TABLE]

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 111For fermionic systems, there is an additional subtlety in identifying the L and R operators. In order to preserve the fermionic commutation relations, one needs to define (say) $V_{L}\equiv V\otimes\mathcal{P}_{R}$ and $V_{R}\equiv I_{L}\otimes V$ (similarly for $W_{L}$ and $W_{R}$), where $\mathcal{P}_{R}$ is the fermionic parity operator acting on the $R$ system. Whenever a two-sided correlator contains an even number of fermionic operators on each side (such as in Eq. (10)), the $\mathcal{P}_{R}$ factors cancel out because $\mathcal{P}_{R}^{2}=I$. In more general cases one has to keep track of this factor, e.g. in deriving the analog of Eqs. (17) and (23) for fermionic operators. We thus have

[TABLE]

where we used general properties of time-reversed states 222Here we assume that time-reversal takes the form $\Theta=\mathcal{K}$ where $\mathcal{K}$ is complex conjugation. In general time-reversal can also include a unitary part $\mathcal{U}$, $\Theta=\mathcal{U}\mathcal{K}$, in which case our discussion still applies with the appropriate insertions of the unitary in the time-ordered correlators., namely $\langle\bar{a}|\bar{b}\rangle=\langle\Theta a|\Theta b\rangle=\langle b|a\rangle$, and

[TABLE]

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^{-\beta(E_{n}+E_{m})/2}\langle n|W(t)V(t^{\prime})|m\rangle=Z_{\beta}\langle n|y^{2}W(t)V(t^{\prime})y^{2}|m\rangle$, with

[TABLE]

This allows us to perform the sum over $n$ using the completness relation $\sum_{n}|n\rangle\langle n|=1$ and arrive at the result

[TABLE]

Note that the trace is now evaluated with respect to the eigenstates of a single-sided Hamiltonian.

For any $t>0$ and $t^{\prime}<0$ Eq. (15) has the structure of an OTOC. In the special case $t^{\prime}=-t$ it becomes

[TABLE]

where we used the time-translation invariance to shift all temporal arguments by $+t$ for clarity.

We observe that $\tilde{F}(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 Maldacena and Stanford (2016); Maldacena et al. (2016). 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 Liao and Galitski (2018); Romero-Bermúdez et al. (2019); Kobrin et al. (2020). The universal upper bound on the Lyapunov exponent $\lambda_{L}$ has also been proven only for regularized OTOCs Maldacena et al. (2016). 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

[TABLE]

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 Shenker and Stanford (2014, 2015); Maldacena et al. (2016).

II.3 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-time-ordered 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 $|{\rm TFD}_{\beta}\rangle$ 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 to $\tilde{F}(t,t^{\prime})$ defined in Eq. (10). In the following we therefore discuss how a state closely approximating $|{\rm TFD}_{\beta}(-t)\rangle$ 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 $|{\rm TFD}_{\beta}\rangle$ as a ground state. A collection of such Hamiltonians were recently constructed in Ref. Cottrell et al. (2019), with a unique $|{\rm TFD}_{\beta}\rangle$ ground state separated from the rest of the spectrum by a gap of order $1/\beta$. 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_{\text{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 $|\Psi_{0}\rangle$ of a simple Hamiltonian

[TABLE]

with

[TABLE]

and $\Theta$ representing the time-reversal operator, has been shown Cottrell et al. (2019) to approximate $|{\rm TFD}_{\beta}\rangle$ to good accuracy for appropriately chosen coefficients $c_{j}$. Here ${\cal 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 Cottrell et al. (2019), and is thus of broad relevance in the study of quantum chaotic systems.

In another recent work, Maldacena and Qi Maldacena and Qi (2018) used an even simpler construction with a coupling of the form

[TABLE]

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 $|{\rm TFD}_{\beta}\rangle$ 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

[TABLE]

This evolution is non-trivial because $|{\rm TFD}_{\beta}\rangle$ 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 $|\Psi_{0}(-t)\rangle\simeq|{\rm TFD}_{\beta}(-t)\rangle$ 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.

II.4 OTOC from two-point functions

Here we discuss an approach that allows one to extract the OTOC from the measurement of a time-ordered two-point function $G_{LR}(t,t^{\prime})$, in a generic chaotic system with coupling $\mu$ 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 $\mu$ before or during the measurement. In contrast to related previous work Gharibyan et al. (2019); Vermersch et al. (2019), 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,

[TABLE]

The average is taken with respect to the (TFD) ground state of the coupled system and is therefore time-translation invariant, $G_{LR}(t,t^{\prime})=G_{LR}(t-t^{\prime})$. The operators inside the average evolve according to the full coupled Hamiltonian $H_{S}=H_{L}+H_{R}+H_{I}$.

At weak coupling $\mu$ (compared to the energy scale $J$ of the chaotic Hamiltonian $H_{0}$), and for short time durations $\mu|t-t^{\prime}|\ll 1$, it is possible to rewrite the two-point correlator in Eq. (22) as a single-sided thermal average of operators evolving according to $H_{L}$ (or, equivalently, $H_{R}$). Formally, this is done by passing from the Heisenberg picture to the interaction picture and expanding the corresponding time-evolution operator $U(t,t^{\prime})$ to leading order in the small parameter $\mu|t-t^{\prime}|$. Details of this calculation are given in Appendix B, and the result is

[TABLE]

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\rangle$ of $H_{L}$. Operators ${\cal 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\to 0$ it coincides with the regularized OTOC. For $s\neq 0$ the trace represents a more general form of the OTOC dependent on two time variables, $F(t_{1},t_{2})={\rm tr}[W(t_{1})Vy^{2}W(t_{2})Vy^{2}]$. In chaotic systems the latter is commonly assumed to behave according to

[TABLE]

with $b(t)$ an even function of $t$ and $b(0)$ real positive Kitaev and Suh (2018); Romero-Bermúdez et al. (2019); Gu and Kitaev (2019). The trace on the last line of Eq. (II.4) is then related to $F(t+s,t-s)$ and, upon integration, becomes proportional to $B(t)e^{\lambda_{L}t}$ with $B(t)=\mu\int_{0}^{t}ds\ b(2s)$.

If we adopt another common assumption Maldacena and Stanford (2016) 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}\ll t\ll\mu^{-1}$) of the $LR$ two-point correlator of the coupled theory should be well approximated by

[TABLE]

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 $\beta$ 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. (II.4) can be derived for the retarded two-point correlator $G_{LR}^{\rm ret}(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.

II.5 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 $\mu$ between the subsystems and short times $|t-t^{\prime}|\ll\mu^{-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 Maldacena and Qi (2018) that realizes a quantum-mechanical dual to an eternal traversable wormhole in (1+1)-dimensional anti-de Sitter spacetime (AdS2) 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 Maldacena and Stanford (2016); Maldacena et al. (2016). Second, the ground state of two such coupled SYK models is well approximated by the TFD state Maldacena and Qi (2018), 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 Danshita et al. (2017); Pikulin and Franz (2017); Chew et al. (2017); Chen et al. (2018a); Franz and Rozali (2018), making it a potentially fruitful platform for laboratory explorations.

III.1 The model

The model introduced by Maldacena and Qi in Ref. Maldacena and Qi (2018) has the form

[TABLE]

where $\mu$ is a constant, and $H_{\alpha}^{\rm SYK}$ with $\alpha=(L,R)$ describe two identical SYK models

[TABLE]

each involving $N$ Majorana zero-mode operators. These obey the usual algebra $\{\chi_{\alpha}^{i},\chi_{\beta}^{j}\}=\delta^{ij}\delta_{\alpha\beta}~{},~{}(\chi_{\alpha}^{j})^{\dagger}=\chi_{\alpha}^{j}$. The coupling constants $J_{ijkl}$ are random, independent Gaussian variables respecting

[TABLE]

and are independent of $\alpha$ – 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

[TABLE]

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 $\Theta={\cal K}$, where ${\cal K}$ denotes complex conjugation. The Majorana operators then transform as $\Theta\chi_{L}^{j}\Theta^{-1}=\chi_{L}^{j}$ and $\Theta\chi_{R}^{j}\Theta^{-1}=-\chi_{R}^{j}$. 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

[TABLE]

where $N_{f}=\sum_{j}c_{j}^{\dagger}c_{j}$ is the total fermion number. What is less obvious is that the fermion number modulo 4, $Q_{4}\equiv\left(N_{f}\right)\;\text{mod}\;4$, is also a symmetry of $H$. This property relies on the perfectly correlated disorder between the two subsystems García-García et al. (2019).

The model can be solved in the limit of large $N$ by methods developed in the context of the original SYK model Sachdev and Ye (1993); Kitaev (2015); Maldacena and Stanford (2016). 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 propagator $G_{\alpha\beta}(\tau_{1},\tau_{2})={1\over N}\sum_{j}\langle{\cal T}\chi_{\alpha}^{j}(\tau_{1})\chi_{\beta}^{j}(\tau_{2})\rangle$. The resulting saddle-point action reads Maldacena and Qi (2018)

[TABLE]

where $S_{0}=-N\ln{\rm Pf}(\delta_{\alpha\beta}\partial_{\tau}-\Sigma_{\alpha\beta})$ and $\Sigma_{\alpha\beta}$ denotes the self energy associated with $G_{\alpha\beta}$. The corresponding saddle-point equations are obtained by varying the action with respect to $G_{\alpha\beta}$ and $\Sigma_{\alpha\beta}$. Using the time-translation invariance, so that $G_{\alpha\beta}(\tau_{1},\tau_{2})=G_{\alpha\beta}(\tau_{1}-\tau_{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}(\tau)$ and $G_{LR}(\tau)$. Their frequency-space counterparts are given as

[TABLE]

where $D(i\omega_{n})=[i\omega_{n}-\Sigma_{LL}(i\omega_{n})]^{2}+[i\mu-\Sigma_{LR}(i\omega_{n})]^{2}$ and $\omega_{n}=\pi T(2n+1)$ is the $n$th Matsubara frequency. The self energies are given by

[TABLE]

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. Maldacena and Stanford (2016); Banerjee and Altman (2017). This yields retarded propagators $G_{\alpha\beta}^{\rm ret}(\omega)$ 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.

III.2 TFD ground state

In Ref. Maldacena and Qi (2018) it was argued that the model in Eq. (26) admits an approximate TFD ground state for all values of the dimensionless parameter $\mu/J$, and an exact TFD ground state in the limits of either small or large $\mu/J$. This can be understood intuitively as follows. For $\mu/J\rightarrow 0$ the ground state of the system is simply given by $|0\rangle_{L}\otimes|0\rangle_{R}$, which coincides trivially with the zero-temperature TFD state $|\rm TFD_{\infty}\rangle$. For $\mu/J\rightarrow\infty$ 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 $\pm\mu$ corresponding to the presence or absence of a fermion in that state. The many-body ground state $|\Psi_{0}\rangle$ of this system is unique and such that all $N$ single-particle states are empty,

[TABLE]

where the $c_{j}$ are defined in Eq. (29). This state is equivalent (see Ref. García-García et al. (2019) for an explicit proof) to the infinite-temperature TFD state

[TABLE]

where $|n\rangle$ are the eigenstates of $H^{\text{SYK}}$. For intermediate values of $\mu/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. Maldacena and Qi (2018); García-García et al. (2019)). The overlap is always greater than $0.96$, with the minimum occurring around $\mu/J\sim 0.1$. This minimum was argued to indicate a phase transition of the Hawking-Page type Hawking and Page (1983) between a wormhole phase at small $\mu/J$ and low temperature, and a black hole phase at large temperature Maldacena and Qi (2018); García-García et al. (2019).

The parameter $\beta$ characterizing $|\rm TFD_{\beta}\rangle$ which best describes the ground state is monotonically decreasing as a function of $\mu/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/\beta$, as expected from the arguments of Ref. Cottrell et al. (2019). In Sec. III.3 we obtain the constant of proportionality as $E_{\rm gap}\approx 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_{\text{gap}}\sim\mu^{2/3}J^{1/3}$ at small $\mu/J$, expected from the wormhole duality and confirmed by solving the imaginary-time SD equations (32) and (33) in Ref Maldacena and Qi (2018). This is presumably due to finite-size effects which become important at energy scales smaller than $\sim 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

[TABLE]

defined using the retarded propagator $G_{LL}^{\mathrm{ret}}(\omega)$, which is related to the Matsubara frequency propagator $G_{LL}(i\omega_{n})$ by the standard analytical continuation $i\omega_{n}\to\omega+i\delta$ Fetter (1971). The spectral function is shown for several values of $\mu$ in Fig. 2d. The spectral gap $E_{\text{gap}}$, defined here as the position of the first peak in $A(\omega)$, is plotted in Fig. 2f. It shows $E_{\mathrm{gap}}=\mu^{2/3}J^{1/3}$ scaling (with numerical prefactor very close to 1) for small $\mu/J$, with a crossover to a linear dependence occurring around $\mu/J\approx 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 $\mu/J$ than was previously reported Maldacena and Qi (2018); García-García et al. (2019). This procedure gives access to the conformal $\sim\mu^{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(\omega)$ displayed in Fig. 2d shows intriguing additional structure, beyond what was reported in previous works. We find a sharp peak at $\omega=E_{\rm gap}$ followed by an sequence of peaks centered close to harmonics of the gap, with spacing $\Delta\omega\sim 3E_{\rm gap}$. The peak at $E_{\rm 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^{\rm ret}(t)$ which shows non-decaying oscillations with a period $2\pi/E_{\rm gap}$ at long times, Fig. 2e. The presence of sharp quasiparticle peaks in $A(\omega)$ 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.

III.3 Measuring OTOCs in coupled SYK models

It is known that, in the limit of $N\rightarrow\infty$ and at strong coupling $\beta J\gg 1$, the SYK model is maximally chaotic with a Lyapunov exponent saturating the chaos bound $\lambda_{L}=2\pi 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 Fu and Sachdev (2016); Pikulin and Franz (2017). Indeed, the exponential growth of the OTOC, parametrized by

[TABLE]

with $A$, $B$ real constants, can be expected for times $J^{-1}\lesssim t<1/\lambda_{L}\log(N/B)$ and $\beta J<N$.

However, previous numerical calculations were carried out using the standard OTOC in Eq. (3) which shows stronger finite-size effects Gu and Kitaev (2019); Kobrin et al. (2020). 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. Shen et al. (2017), we define the lower bound of the fitting region by a time $t_{-}$ such that $F(t_{-})\sim 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^{\prime\prime}(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 $\lambda_{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 Fu and Sachdev (2016); Pikulin and Franz (2017)). 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. Maldacena and Stanford (2016); Banerjee and Altman (2017)) more accurately due to weaker finite-size effects (see also Ref Kobrin et al. (2020)).

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. Li et al. (2017); Gärttner et al. (2017); Landsman et al. (2019).

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^{\prime})$ in the ground state of the coupled system. This is clearly a more straightforward measurement, but is limited to weak couplings $\mu/J$. To verify that this approach indeed works we adapt Eq. (II.4) to the Maldacena-Qi model by identifying ${\cal O}^{j}=\chi^{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,

[TABLE]

valid for $0<t\ll\mu^{-1}$. Similar to the time-ordered case, $G_{LR}^{\rm ret}(t)$ contains an OTOC contribution (last line of Eq. (38)), and we therefore expect an exponential growth at intermediate times.

In Fig. 4a we show the imaginary part of $G^{\rm ret}_{LR}(t)$ calculated numerically from the large-$N$ saddle-point equations for several values of $\mu/J$. For sufficiently weak couplings $\mu/J$, we can fit an approximately exponential growth in the expected regime, from which we extract a putative Lyapunov exponent $\lambda_{L}(\mu)$ as shown in the inset of Fig. 4a. The extracted exponents follow the $\sim\mu^{2/3}$ scaling of the energy gap (see Fig. 2f) at small $\mu/J$. Given that $E_{\rm gap}$ scales linearly with the effective temperature $1/\beta$ of the corresponding TFD state Cottrell et al. (2019) (see Fig. 2c), our results imply that $\lambda_{L}\sim 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 $\lambda_{L}(T)$ which requires the knowledge of the function $T(\mu)$. This can be in principle obtained from our ED results shown in Fig. 2b. However, because ED does not accurately capture the $E_{\rm gap}\sim\mu^{2/3}$ scaling at small $\mu/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(\mu)$ dependence directly from the large-$N$ formalism which correctly captures the $E_{\rm gap}\sim\mu^{2/3}$ scaling. To do this we apply Eq. (6) with ${\cal O}_{L}=H_{L}^{\rm SYK}$ to the Maldacena-Qi model, obtaining

[TABLE]

The left-hand side is evaluated in the ground state of the Maldacena-Qi model and gives $\langle H_{L}\rangle_{\rm TFD}$ as a function of $\mu$ (dropping the SYK superscript from here on). The right hand side is evaluated in the thermal ensemble of a single SYK model and gives $\langle H\rangle_{\beta}$ as a function of temperature. Matching these two energies through Eq. (39) then yields the required function $T(\mu)$.

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 Fetter (1971) applied to the Maldacena-Qi Hamiltonian yields

[TABLE]

where $G_{\alpha\beta}(\tau)$ is the imaginary-time Green’s function. Fourier transforming into the Matsubara frequency space and using the spectral representation of $G_{\alpha\beta}(i\omega_{n})$ this can be rewritten in the integral form

[TABLE]

which is convenient for numerical evaluation. Here $n(\omega)=1/(e^{\beta\omega}+1)$ denotes the Fermi-Dirac distribution and $\rho_{\alpha\beta}(\omega)$ 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 $\mu=0$ and $\rho_{LL}$ replaced by the spectral function $\rho(\omega)$ of the single SYK model. The results of this calculation are summarized in Fig. 4b. We find that $T=\mu$ for large $\mu/J$ while $T/J\approx 0.76(\mu/J)^{2/3}$ at small $\mu/J$. Using this result we obtain Lyapunov exponents in good agreement with the prediction for maximal chaos at low temperatures, $\lambda_{L}=2\pi 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 $\lambda_{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 $\mu$ 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.

IV.1 Realizations of coupled SYK models

IV.1.1 Quantum dots

Perhaps the most conceptually transparent realization of the Maldacena-Qi model Maldacena and Qi (2018) 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 Alicea (2012); Beenakker (2013); Leijnse and Flensberg (2012); Stanescu and Tewari (2013); Elliott and Franz (2015). 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 interactions between the underlying electrons are taken into account Chew et al. (2017). 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 $j$th wire this coupling has the form $i\mu\chi_{L}^{j}\chi_{R}^{j}$ with $\mu\sim e^{-L/\xi}\cos{(k_{F}L)}$, where $\xi$ denotes the superconducting coherence length and $k_{F}$ is the Fermi wavevector of electrons in the wire. Both $\xi$ and $k_{F}$ are sensitive to the gate voltage applied to the wire, which makes the coupling strength $\mu$ tunable, at least in principle. A similar term arises for long wires $L\gg\xi$ 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. Chew et al. (2017), 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 Mourik et al. (2012) has been confirmed and extended by several other groups Das et al. (2012); Deng et al. (2012); Rokhinson et al. (2012); Finck et al. (2013); Deng et al. (2016); Zhang et al. (2018). 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. Chew et al. (2017), 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.

IV.1.2 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) Fu and Kane (2008). 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 Pikulin and Franz (2017). 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 Pikulin and Franz (2017); Lantagne-Hurtubise et al. (2018). 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 Bi2Se3), 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\mu\sum_{j}\chi_{L}^{j}\chi_{R}^{j}$, as required for the Maldacena-Qi model. The coupling strength $\mu$ 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 ${\mathrm{Bi}}_{2}{\mathrm{Te}}_{3}/{\mathrm{NbSe}}_{2}$ heterostructures Xu et al. (2015); Sun et al. (2016), 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 FeTe0.55Se0.45 Wang et al. (2018); Chen et al. (2018b); Liu et al. (2018); Machida et al. (2019); Chiu et al. (2019). 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}\simeq 14$K.

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.

IV.1.3 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 Sachdev and Ye (1993); Sachdev (2015). It is defined by the Hamiltonian

[TABLE]

where $c_{j}$ annihilates a complex fermion and $\tilde{\mu}$ 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 Chen et al. (2018a). 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.

IV.2 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 $\mu$ on microscopic timescales (i.e. times of order $\hbar/J$), or we do not. In the first case we can manipulate $\mu$ to prepare the initial resource state ${|\rm TFD}_{\beta}(-t)\rangle$, 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 $\mu$ cannot be controlled on microscopic timescales, it is still possible to extract the OTOC by measuring $G_{LR}^{\mathrm{ret}}(t)$ in a system with constant nonzero $\mu$, 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. (II.4).

IV.2.1 When $\mu$ can be controlled

Following existing theoretical work on eternal traversable wormholes Maldacena and Qi (2018); Cottrell et al. (2019), 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(\lambda)=H_{0}+\lambda H_{I}$. Measuring the OTOC requires the TFD state evolved to negative time, $|{\rm TFD}_{\beta}(-t)\rangle$, which we demonstrate in Appendix A can be achieved, for short time durations at least, by tuning the dimensionless coupling $\lambda$. 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_{\text{phys}}$ that is much smaller than the energy gap of the combined system, which puts it into its ground state well approximated by $|{\rm TFD}_{\beta}(0)\rangle$. (iii) Increase coupling $\lambda$ to a value larger than 1 for a short period of time. This creates a good approximation of $|{\rm TFD}_{\beta}(-t)\rangle$ where the time evolution is with respect to $H_{0}$, cf. Appendix A. (iv) Decouple the system by setting $\lambda=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. Plugge et al. (2017); Karzig et al. (2017) 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)\to[a(t)+a^{\dagger}(t)]$, and assuming roughly isotropic capacitive coupling parameters $\sim g_{j}\simeq g$ to all nanowires, one finds

[TABLE]

with a single photon species $a$ and total nanowire-charge operator $\hat{Q}(t)=\sum_{j}\hat{q}_{j}(t)=\sum_{j}i\chi_{L}^{j}(t)\chi_{R}^{j}(t)$. The term $H_{\mathrm{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 equal-time left-right Green’s function, as $G_{LR}(t)=-iQ(t)/N$.

IV.2.2 When $\mu$ is fixed

In Secs. II-D and III-D, we showed that the averaged $LR$ Green’s function of the coupled theory at fixed $\mu$ contains information on the OTOC for small couplings $\mu/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 $\rho_{x}(\omega)$ at point $x$ and frequency $\omega=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

[TABLE]

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_{LR}^{\rm ret}(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_{LR}^{\rm ret}(t)$ to a real quantity. Second, Majorana zero mode operators $\chi^{k}_{L/R}$ are simply related to the electron operators in the wire through the solution of the relevant Bogoliubov-de Gennes equation Alicea (2012); Beenakker (2013); Leijnse and Flensberg (2012); Stanescu and Tewari (2013); Elliott and Franz (2015), which is largely dictated by symmetries of the setup. This implies proportionality of $iG_{LR}^{\rm ret}(t)$ to the electron spectral function, Fourier-transformed 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 Gärttner et al. (2017); Li et al. (2017); Landsman et al. (2019); Swingle (2018), 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 Maldacena et al. (2016); Gao et al. (2017); Maldacena et al. (2017), we proposed and tested new protocols for OTOC measurement, where the preparation of a specific resource state – namely the TFD – replaces 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 Maldacena and Qi (2018); Gao and Liu (2019); Cottrell et al. (2019). 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 Liao and Galitski (2018); Romero-Bermúdez et al. (2019); Kobrin et al. (2020), 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. Yao et al. (2016) 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. (II.4) and (25), that the Lyapunov exponent $\lambda_{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 $\lambda_{L}\approx 2\pi T$ consistent with the presence of maximal chaos, saturating the Maldacena-Shenker-Stanford bound Maldacena et al. (2016) $\lambda_{L}\leq 2\pi T$ in the weak coupling limit $\mu/J\ll 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^{\rm ret}_{LR}(t)$ is found to be proportional to the Fourier transform of the electron spectral function $\rho_{x}(\omega)$ (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 Maldacena and Qi (2018) 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 Kitaev (2015); Maldacena and Stanford (2016).

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. Pikulin and Franz (2017); Chew et al. (2017); Chen et al. (2018a). 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 $\Theta$ in a systems of spin-${1\over 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 Sachdev and Ye (1993); Sachdev (2015). 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. Cottrell et al. (2019).

Further, the generality of the construction in Ref. Cottrell et al. (2019) 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.

Acknowledgements

The authors are indebted to J. Alicea, V. Galitski, F. Haehl, A. Kitaev, B. Kobrin, X.L. Qi, C. Li, J. Maldacena, S. Sahoo and B. Swingle for useful and stimulating discussions. We would like to thank B. Kobrin and coworkers for pointing us to subtleties in the numerical analysis of (regularized) OTOCs in the SYK model Kobrin et al. (2020). We thank NSERC and CIFAR for financial support. SP and MF are grateful to KITP for hospitality during the conference and program “Order from Chaos”, where part of the research was conducted with support of the National Science Foundation under Grant No. NSF PHY-1748958. Further SP is grateful to the Aspen Center for Physics, supported by National Science Foundation grant PHY-1607611, for hospitality during the conference “Many-Body Quantum Chaos”.

Appendix A TFD preparation at negative times

To demonstrate the preparation of an initial TFD state $|\Psi_{0}(-t)\rangle$ at negative times, we consider a generalization of the coupled systems’ Hamiltonian $H_{S}$ to

[TABLE]

where $\lambda$ is a dimensionless parameter used to control the strength of coupling between the $L$ and $R$ systems. We have $H(1)=H_{S}$ and $H(0)=H_{0}$. First assume that the system has been cooled down and is in the ground state $|\Psi_{0}\rangle$ of $H(1)=H_{S}$,

[TABLE]

Now imagine we increase $\lambda$ to a value larger than one. The system will start evolving according to $|\Psi_{\lambda}(t)\rangle=e^{-iH(\lambda)t}|\Psi_{0}\rangle$. Using $H(\lambda)=\lambda H_{S}-(\lambda-1)H_{0}$, which follows from Eq. (45), this evolution can be rewritten as

[TABLE]

We next employ the Baker-Campbell-Hausdorff formula $e^{A}e^{B}=e^{A+B+{1\over 2}[A,B]+\dots}$ to separate the two terms in the exponential, and multiply from the right by $e^{-B}$ to obtain

[TABLE]

where dots represent higher-order commutators. Taking $A=-i\lambda H_{S}t+i(\lambda-1)H_{0}t$ and $B=i\lambda H_{S}t$, Eq. (47) becomes

[TABLE]

Using Eq. (46), the last exponential evaluates to $e^{-i\lambda\epsilon_{0}t}$. For short time durations one can furthermore neglect the $t^{2}$ term in the first exponential, which leads to

[TABLE]

We see that increasing the coupling strength $\lambda$ to a value larger than 1 has the same effect as evolving the state $|\Psi_{0}\rangle$ backward in time under the decoupled Hamiltonian $H_{0}=H_{L}+H_{R}$. Manipulating the coupling strength $\lambda$ 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. (50) remains true for sufficiently short times $t$ such that one can neglect the commutator term in the exponential of Eq. (49). This can be estimated from the condition

[TABLE]

If we assume that the relevant energy scales for $H_{0}$ and $H_{I}$ are $J$ and $\mu$, respectively, then $\langle H_{0}\rangle\sim J$ and $\langle i[H_{0},H_{I}]\rangle\sim J\mu$. Eq. (51) hence yields a constraint

[TABLE]

on the time duration over which the evolution backward in time following Eq. (50) can be achieved. For weakly coupled systems $(\mu\ll J)$, this constraint gives a sufficient window to probe OTOCs using the method described in Sec. II–B.

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\sigma^{z}$ and $H_{I}=-\mu\sigma^{x}$. The ground state $|\Psi_{0}\rangle$ of the combined Hamiltonian $H(\lambda)=H_{0}+\lambda H_{I}$ with $\lambda=1$ has the spin pointing along the direction parallel to the total magnetic field ${\bm{B}}=(\mu,0,J)$. If we switch off $H_{I}$ the spin will start precessing counterclockwise around the field direction ${\bm{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 $\lambda$, the spin will start precessing counterclockwise around the new field direction ${\bm{B}}_{\lambda}=(\lambda\mu,0,J)$, as illustrated in Fig. 7a. At short times $t$ we observe that the evolution for $\lambda>1$ approximates backward time evolution under $H_{0}$. This effect can be quantified by calculating the overlap

[TABLE]

between the ground state evolved backward in time according to $H_{0}$, i.e. $|\Psi_{0}(-t)\rangle=e^{+iH_{0}t}|\Psi_{0}\rangle$ and the ground state evolved forward in time according to $H(\lambda)$, $|\Psi_{\lambda}(t)\rangle=e^{-iH(\lambda)t}|\Psi_{0}\rangle$. Elementary but somewhat tedious calculation gives an explicit expression for $P(t)$ that we plot in Fig. 7b, for several values of $\mu$. We observe that for short times and $\mu\ll J$ the overlap remains very close to 1, confirming that the method indeed yields an excellent approximation to the state $|\Psi_{0}\rangle$ 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 $|\Psi_{0}(-t)\rangle$ in a wide variety of settings.

Appendix B Short-time expansion of LR two-point correlator

In this Appendix we derive Eq. (II.4), which we used in the main text to argue that the $LR$ two-point correlator $G_{LR}(t,t^{\prime})$, defined in Eq. (22), contains at short times information on the OTOC. We proceed by evaluating the two-point correlator

[TABLE]

Here $|\Psi_{0}\rangle$ 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 $|{\rm TFD}_{\beta}\rangle$. We work in the Heisenberg picture where $V_{\alpha}(t)$ is an arbitrary Hermitian operator evolving according to the full Hamiltonian $H$.

$G_{LR}(t,t^{\prime})$ 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 $\beta$ of a single, decoupled system. Mathematically, the goal is to express $G(t,t^{\prime})$ 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 Fetter (1971)

[TABLE]

where superscript $I$ denotes the interaction picture and $U(t,t^{\prime})=e^{iH_{0}t}e^{-iH(t-t^{\prime})}e^{-iH_{0}t^{\prime}}$ is the unitary operator that translates between the Heisenberg and interaction pictures. Using Eq. (55), the correlator becomes

[TABLE]

where we used the property $U(t,s)U(s,t^{\prime})=U(t,t^{\prime})$. For simplicity, we henceforth also assume that $t>t^{\prime}$.

We now employ a standard result of diagrammatic many-body theory Fetter (1971) that express $U(t,t^{\prime})$ as a series expansion in $H_{I}(t)$ of the form

[TABLE]

Here the time evolution of $H_{I}(t)$ is according to $H_{0}$. Substituting this into Eq. (56), and retaining only terms up to first order in $H_{I}(t)$, we find

[TABLE]

where $\langle\dots\rangle_{0}$ denotes the expectation value with respect to the ground state $|\Psi_{0}\rangle$. This expression is valid when one can neglect all higher order terms in the expansion (B) of $U(t,t^{\prime})$, $U(0,t)$ and $U(t^{\prime},0)$. This requires short time durations $|t-t^{\prime}|$ as well as individually small $|t|$ and $|t^{\prime}|$. In the following we focus on the symmetric case $t^{\prime}=-t$, which has a convenient property that short duration $|t-t^{\prime}|$ automatically assures that $|t|$ and $|t^{\prime}|$ are small.

As the final step we substitute for $H_{I}$ the Maldacena-Qi form given in Eq. (20), and approximate $|\Psi_{0}\rangle$ by $|{\rm TFD}_{\beta}\rangle$. With these choices, the expectation values in Eq. (B) 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

[TABLE]

where $\eta=+/-$ 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\rangle$ of the same Hamiltonian. Each individual trace in Eq. (B) is time-translation 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. (II.4).

Mathematically, Eq. (B) can be viewed as an expansion of the propagator $G$ in powers of a dimensionless time variable $\tilde{t}=\mu t$ to first order. Higher order contributions that would result from the omitted terms in Eq. (B) can be neglected when $\mu t\ll 1$ which constrains the expected domain of validity of Eq. (B) to short times or small values of coupling $\mu$.

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

[TABLE]

evaluated in the ground state $|\Psi_{0}\rangle$ of the system. Expanding the anticommutator and using the basic property of the inner product $\langle a|b\rangle=\langle b|a\rangle^{*}$, we can rewrite the average as $\langle A(t)B(t^{\prime})\rangle_{0}+\langle[B(t^{\prime})A(t)]^{\dagger}\rangle^{\ast}_{0}$. For Hermitian operators $A$ and $B$ this equals $2{\rm Re}\langle A(t)B(t^{\prime})\rangle_{0}$, and $G^{\rm ret}(t,t^{\prime})$ is therefore purely imaginary.

As we already mentioned, a tunneling measurement can be used to extract the electron spectral function $\rho_{x}(\omega)$ which is related to the electron propagator ${\cal G}_{x}(\omega)$. We therefore start by considering the corresponding quantity defined in the time domain as

[TABLE]

where $c_{x}(t)$ annihilates electron at time $t$ and at point $x$ of the wire. We want to relate ${\cal G}(t)$ to the Majorana propagators $G_{\alpha\beta}(t)$ discussed in the main text. To this end recall the relations

[TABLE]

which follow from the solution of the Bogoliubov-de Gennes (BdG) equations for the Majorana wire Alicea (2012); Beenakker (2013); Leijnse and Flensberg (2012); Stanescu and Tewari (2013); Elliott and Franz (2015). Here $\Phi_{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. (62) is constrained by the choice that $\chi_{L/R}$ transform as even/odd under time reversal.

We may invert Eqs.  (62) to express the electron operators in terms of the Majoranas,

[TABLE]

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. (61) and neglecting these terms, we find the Majorana fermion contribution to the electron Green’s function

[TABLE]

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 $i{\cal G}^{\rm ret}_{x}(t)=\theta(t)\langle\{c_{x}(t),c_{x}^{\dagger}(0)\}\rangle_{0}$, with a similar result:

[TABLE]

Given that $G_{LL}^{\rm ret}(t)$, $G_{LR}^{\rm ret}(t)$ are imaginary while $\Phi_{L/R}(x)$ are real, we find

[TABLE]

To complete this calculation, it remains to relate ${\rm Re}[{\cal G}_{x}^{\rm ret}(t)]$ to the electron spectral function which is the observable quantity. We use the spectral representation ${\cal G}_{x}^{\rm ret}(\omega)=\int_{-\infty}^{\infty}d\omega^{\prime}\rho_{x}(\omega^{\prime})/(\omega-\omega^{\prime}+i\delta)$, which, upon Fourier transforming, gives

[TABLE]

Taking the real part and recalling that $\rho_{x}(\omega)$ is strictly real, we obtain

[TABLE]

Finally combining with Eq. (66), 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 Chew et al. (2017). 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. Chew et al. (2017).

The system in Fig. 5 is described by $N$ sets of Majorana modes $\chi_{L}^{j},~{}\chi_{R}^{j}$ at the left and right ends of the nanowires, and by the $N_{L/R}\gg N$ complex fermions $c_{s,\alpha=L/R}$ hosted in disordered wave-functions of its quantum dots. We write the dot fermions in Majorana representation $c_{s,\alpha}=(\eta_{s,\alpha}+i\tilde{\eta}_{s,\alpha})/2$, where $\eta_{s,\alpha}$ is even under time-reversal (TR) while $\tilde{\eta}_{s,\alpha}$ is odd. Similarly we take nanowire Majoranas $\chi_{L}^{j}$ and $\chi_{R}^{j}$ to be even and odd under TR, respectively. Assuming that both quantum dots preserve the BDI symmetry class of the Majorana sector Chew et al. (2017), 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 $\sim i\chi_{L}^{j}\chi_{R}^{j}$ 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 $\vec{\chi}_{L/R}$, $\vec{\eta}_{L/R}$ and $\vec{\tilde{\eta}}_{L/R}$ for the respective left/right groups of Majorana fermions. Following Ref. Chew et al. (2017), we first describe the hybridization of Majoranas $\vec{\chi}_{L/R}$ into the left and right quantum dots as

[TABLE]

Here $M_{\lambda\alpha}$ are real rectangular matrices of couplings $\sim\lambda$ between the wire and dot Majoranas, cf. Fig. 5, and $M_{\epsilon\alpha}$ are real $N_{\alpha}\times N_{\alpha}$ matrices encoding the dots level structure $\sim\epsilon$. The coupling Hamiltonians $H_{0,\alpha=L/R}$ can be diagonalized by orthogonal rotations obtained from a singular-value decomposition of the coupling matrix encoded in $M_{\lambda\alpha}$ and $M_{\epsilon\alpha}$ Chew et al. (2017), giving new Majorana spinors

[TABLE]

and $\vec{\tilde{\eta}}_{L}^{\prime}={\cal O}_{\tilde{\eta},L}\vec{\tilde{\eta}}_{L}$, and similar for right Majorana operators with $L\to R$ and $\vec{\eta}\leftrightarrow\vec{\tilde{\eta}}$. In these new operators,

[TABLE]

where $\eta_{j,\alpha}^{\prime}$ and $\tilde{\eta}_{j,\alpha}^{\prime}$ correspond to the rotated spinors $\vec{\eta}_{\alpha}^{\prime}$ and $\vec{\tilde{\eta}}_{\alpha}^{\prime}$, and $E_{j,\alpha}$ 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 $\vec{\chi}_{\alpha}^{\prime}$ that do not appear in $H_{0,\alpha}$. Further, assuming strong wire-dot hybridizations $\lambda\gg N\delta\epsilon$, where $\delta\epsilon$ is the typical level spacing in the dot, all wire Majoranas $\vec{\chi}_{\alpha}$ are absorbed into the respective left and right dots. Finite-energy modes in Eq. (71) are split off to energies $\sim N\delta\epsilon$ by level repulsion.

We now add the $LR$ coupling between the two SYK-dots. In the original basis, sensible couplings are pairwise between Majoranas $\chi_{L}^{j}$ and $\chi_{R}^{j}$ on each wire, yielding $H_{\mathrm{int}}=i\sum_{j=1}^{N}\mu_{j}\chi_{L}^{j}\chi_{R}^{j}=i\vec{\chi}_{L}^{T}\mu_{LR}\vec{\chi}_{R}$ with diagonal coupling matrix $\mu_{LR}$. The specific form and possible tuning of couplings $\mu_{j}$ is discussed below. Inverting the orthogonal transformations in Eq. (70), one can represent the original wire Majoranas in the new variables as $\vec{\chi}_{L}={\cal O}_{\chi\chi,L}^{T}\vec{\chi}_{L}^{\prime}+{\cal O}_{\eta\chi,L}^{T}\vec{\eta}_{L}^{\prime}$ and $\vec{\chi}_{R}={\cal O}_{\chi\chi,R}^{T}\vec{\chi}_{R}^{\prime}+{\cal O}_{\tilde{\eta}\chi,R}^{T}\vec{\tilde{\eta}}_{R}^{\prime}$ and obtain

[TABLE]

Here the new $LR$ coupling matrices are given as $\mu_{\chi\chi}={\cal O}_{\chi\chi,L}\mu_{LR}{\cal O}_{\chi\chi,R}^{T}$ etc., connecting different species of the new Majorana operators. Assuming that terms $\sim\mu_{j}$ are small against wire-dot hybridizations and the on-site level repulsion scale, even in their presence the gapped Majorana-pairs in Eq. (71) will stay far-removed. The remaining effective Hamiltonian then reads

[TABLE]

In the absence of left-right symmetry of the underlying disordered quantum dots in Fig. 5, there is no reason to assume that $\mu_{\chi\chi}$ is diagonal. Nevertheless, it is always possible to re-diagonalize the $LR$ Hamiltonian by re-applying the orthogonal transformations to Majoranas $\vec{\chi}_{L/R}^{\prime}$. Given that the original coupling matrix $\mu_{LR}$ was diagonal, one finds $H_{\mathrm{int,eff}}=i\vec{\chi}_{L}^{\prime\prime T}\mu_{LR}\vec{\chi}_{R}^{\prime\prime}$ with

[TABLE]

and similarly $\vec{\chi}_{R}^{\prime\prime}={\cal O}_{\chi\chi,R}^{T}\vec{\chi}_{R}^{\prime}$. Note ${\cal O}_{\chi\chi,s}^{T}{\cal O}_{\chi\chi,s}\neq 1$, since only the full transformation in (70) is orthogonal.

We hence observe that under the approximations taken above, the interaction matrix $\mu_{LR}=\mathrm{diag}(\mu_{j})$ remains unchanged by the absorption of Majorana modes $\chi_{j,L/R}$ into the SYK dots, but now acts on the new modes $\vec{\chi}_{L/R}^{\prime\prime}$.

Finally, for ease of notation, we relabel $\vec{\chi}_{L/R}^{\prime\prime}\to\vec{\chi}_{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,\alpha=L/R}$ Chew et al. (2017). Generically one then obtains a Hamiltonian $H=H_{L}+H_{R}+H_{\mathrm{int,eff}}$ with site-dependent four-Majorana interactions as in Eq. (27)

[TABLE]

From here, assuming left-right symmetry of the SYK-dots in Fig. 5 such that $J_{ijkl}^{L}=J_{ijkl}^{R}$ identically, one obtains the Maldacena-Qi Hamiltonian Maldacena and Qi (2018) in Sec. III.

We now discuss how one may realize (tunable) one-to-one bilinear couplings across nanowires of the device in Fig. 5. While one may use residual Majorana hybridizations $\mu_{j}\sim\mu_{0}e^{-L/\xi}$ 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 $\chi_{L}^{j},~{}\chi_{R}^{j}$ van Heck et al. (2012). On the Hamiltonian level, this term can be incorporated as

[TABLE]

where $\mu_{0}$ depends on both $E_{c}$ and $E_{J}$ van Heck et al. (2012), 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 Majorana-coupling across all pairs of modes $\chi_{L/R}^{j}$. 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 coupled-wire 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 parities $\hat{q}_{j}=i\chi_{L}^{j}\chi_{R}^{j}$ 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. Plugge et al. (2017); Karzig et al. (2017).

To illustrate how this readout works, assume that a gate parameter $n$ in Eq. (76) is set such that the corresponding interaction is nearly switched off, $n(t)=\frac{1}{2}+v(t)$ with a small fluctuating gate voltage $|v(t)|\ll 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)\to[a(t)+a^{\dagger}(t)]$. The total nanowire-charge readout setup then is described by

[TABLE]

where $\hat{q}(t)$ is the time-evolving nanowire charge parity. $H_{\mathrm{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 $\langle a(t)\rangle$, and via Eq. (77) also on $\langle\hat{q}(t)\rangle$ Plugge et al. (2017). 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.

With similar but somewhat more complex measurements, monitoring e.g. the time-dependent resonator photon Green’s function $B_{a}(t,t^{\prime})=\langle a^{\dagger}(t)a(t^{\prime})\rangle$ (or any quantity quadratic in photon operators $a,~{}a^{\dagger}$), one finds

[TABLE]

with nanowire-charge correlator $D_{q}(t,t^{\prime})=\langle\hat{q}(t)\hat{q}(t^{\prime})\rangle$. Here $b_{0}(t,t^{\prime})=b_{0}(t-t^{\prime})$ is a bare photon Green’s function of the uncoupled resonator Hamiltonian $H_{\mathrm{res}}$. Note that a time-dependent Green’s functions as in Eq. (78) is encoded by time- and frequency-resolved resonator occupations $B_{a}(\tau,\omega)=\int d\tau^{\prime}e^{i\omega\tau^{\prime}}B_{a}(\tau+\frac{\tau^{\prime}}{2},\tau-\frac{\tau^{\prime}}{2})$, that can be measured in principle. Ignoring back-action of the resonator on the system (that generates an effective interaction as in Eq. (76)), and to lowest order in charge-resonator interaction $g$, this gives information about $d_{0}(t,t^{\prime})=\langle\hat{q}(t)\hat{q}(t^{\prime})\rangle_{0}$. The latter expression constitutes a four-Majorana two-sided correlator evaluated with respect to the bare SYK Hamiltonian $H=H^{\mathrm{SYK}}_{L}+H^{\mathrm{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.

Appendix F Numerical details

F.1 Exact diagonalization

The numerical exact diagonalization of the Maldacena-Qi model proceeds in a standard fashion. One first defines a complex fermion basis in order to build the Hilbert space of dimension $2^{N}$ (for $2N$ Majorana zero-modes). The most convenient choice, as discussed in the main text, is to define complex fermions delocalized accross the two subsystems,

[TABLE]

This basis has two key advantages: first, the anti-unitary time-reversal operator takes the simple form $\Theta=\mathcal{K}$, which makes the Hamiltonian $H_{S}$ explicitly real and thus saves computational resources; second, the infinite-temperature TFD takes the simple form

[TABLE]

reflecting the absence of fermions in all modes $j$. In practice, one can then generate the TFDs for finite $\beta$ through the relationship Maldacena and Qi (2018)

[TABLE]

where $H_{0}=H_{L}+H_{R}$ is the Hamiltonian of the system with zero coupling, $\mu=0$.

This procedure can be made more efficient by implementing the symmetries of the problem: fermion parity $P$ or – even better – fermion number modulo 4 ($Q_{4}$), cf. Ref. García-García et al. (2019), which includes fermion parity. The Hamiltonian then splits into four blocks of unequal dimensions. The ground state $|G\rangle$ is always found to be in the $Q_{4}=0$ sector, and any thermofield double state constructed using Eq. (81) also has $Q_{4}=0$ because $[Q_{4},H_{0}]=0$.

F.2 Solution of Schwinger-Dyson equations

Here we discuss how to solve the large-$N$ Schwinger-Dyson (SD) equations pertaining to the two coupled SYK models, cf. Sec. III-A Eqs. (31)-(33), in a numerically efficient way. To this end, we first perform analytical manipulations to implement as many of their subtle symmetries as possible.

Since we are interested in real-time dynamics of the SD equations, we work in real frequency $\omega$ and time $t$. To obtain the relevant retarded propagators we apply the standard analytical continuation $i\omega_{n}\to\omega+i\delta$ to Eqs. (31) and write

[TABLE]

with $D(\omega)=\left[\omega-\Sigma^{\rm ret}_{LL}\right]^{2}+\left[i\mu-\Sigma^{\rm ret}_{LR}\right]^{2}$. The retarded self-energies follow from Eq. (33) as

[TABLE]

We then employ the method introduced by Banerjee and Altman Banerjee and Altman (2017) (supplement S2 therein) to express the self energies as

[TABLE]

The factors $n_{ss^{\prime}}(t)$ are calculated from the spectral representation of the corresponding propagators $G^{\rm ret}_{ab}(\omega)$ as

[TABLE]

where $n_{F}(\omega)=1/(e^{\beta\omega}+1)$ is the Fermi function.

Given certain symmetries of the spectral functions $\rho_{LL}(\omega)$ and $\rho_{LR}(\omega)$, we now show how to deduce all occupations $n_{ss^{\prime}}(t)$ from just a single one. The spectral functions in our case read

[TABLE]

While the first line is standard, the $i$ factor on the second line is unconventional but comes about to give the correct result for $\rho_{LR}(\omega)$ within our conventions. Consider as an example the retarded Green’s functions for the non-interacting case

[TABLE]

which solve Eqs. (82) and (83) when $J=0$. The spectral functions then show Lorentzian peaks (with weight $\frac{1}{2}$)

[TABLE]

centered at $\omega=\pm\mu$, and are symmetric (anti-symmetric) around $\omega=0$, where $\sigma_{LL}=+1$ ($\sigma_{LR}=-1$).

In our numerical solution we begin from the ansatz (87) with small non-zero $\delta$ and iterate together with Eqs. (85) and (84). Given that the above symmetry properties persist throughout the iterations of the SD equations, one can relate factors $n_{ss^{\prime}}(t)$ as

[TABLE]

The first two of above equations hold since both spectral and Fermi functions are real, while the last two require the knowledge about the $\omega$-parity $\sigma_{ab}$ of the spectral function. First and third equations yield $n_{--}$ and $n_{-+}$ from $n_{++}$. Inserting the third equation into the second, one further obtains $n_{+-}(t)=\sigma_{ab}[n_{++}(t)]^{\ast}$.

Hence only $n(t)=n_{++}(t)$ is required to evaluate the retarded self-energy $\Sigma^{\rm ret}_{ab}(\omega)$. Henceforth we only consider this quantity and attach to it an index $n_{ab}(t)$ to denote its relation to the specific spectral function $\rho_{ab}(\omega)$. Note that $n_{ab}(t)$ is not real in general; it is the Fourier transform of a real function $\rho_{ab}(\omega)n_{F}(\omega)$, so it is Hermitian on the time-axis: $n_{ab}(t)=n_{ab}^{\ast}(-t)$.

Since we are working with imaginary $G^{\rm ret}_{LR}(\omega)$ (leading also to the unconventional spectral function equation), one has to put back a factor of $-i$ in the self-energy Eq. (84) above, $\Sigma^{\rm ret}_{LR}\to\tilde{\Sigma}^{\rm ret}_{LR}=-i\Sigma^{\rm ret}_{LR}$. With this additional factor, one can also rephrase the self-energies as

[TABLE]

This form makes apparent the explicit anti-Hermiticity of $\Sigma^{\rm ret}_{ab}(\omega)$ on the frequency axis: the argument of the Fourier transform in Eq. (90) is real, hence the transformed function is Hermitian. With the prefactors $\mp 2iJ^{2}$, $\Sigma^{\rm ret}_{LL}(\omega)$ and $\tilde{\Sigma}^{\rm ret}_{LR}(\omega)$ are manifestly anti-Hermitian: $\Sigma^{\rm ret}_{LL}(\omega)^{\ast}=-\Sigma^{\rm ret}_{LL}(-\omega)$ and $\tilde{\Sigma}^{\rm ret}_{LR}(\omega)^{\ast}=-\tilde{\Sigma}^{\rm ret}_{LR}(-\omega)$. Since both the bare non-interacting propagator (terms $\omega$ and $i\mu$ in Eq. (82)) and $\Sigma^{\rm ret}_{ab}(\omega)$ are anti-Hermitian, the same holds true for any odd-power product of the two emerging from the Dyson series. Hence also the full propagator $[G^{\rm ret}_{ab}(\omega)]^{\ast}=-G^{\rm ret}_{ab}(-\omega)$ is anti-Hermitian on the frequency axis, and $G^{\rm ret}_{ab}(t)$ is strictly imaginary.

After performing the above manipulations, we further find it convenient to switch to Green’s functions and self-energies that diagonalize the SD equation explicitly:

[TABLE]

The spectral densities $\rho_{\pm}=\rho_{LL}\pm\rho_{LR}$ here follow from the standard relation $\rho_{\pm}(\omega)=-\frac{1}{\pi}\mathrm{Im}G^{\rm ret}_{\pm}(\omega)$. From our ansatz in Eq. (87) one sees that $g_{+}(\omega)=-g_{-}^{\ast}(-\omega)$ which implies $\rho_{+}(\omega)=\rho_{-}(-\omega)$ initially. One then finds $\Sigma_{+}(\omega)=-\Sigma_{-}(-\omega)^{\ast}$, and consequently also $G^{\rm ret}_{+}(\omega)=-G^{\rm ret}_{-}(-\omega)^{\ast}$ which a-priori is not obvious. Hence the symmetry $\rho_{+}(\omega)=\rho_{-}(-\omega)$ is kept throughout, and the full solution of the SD equations is encoded in

[TABLE]

where broadening and shift of resonances in $G^{\rm ret}_{+}(\omega)$ is directly linked to the real/imaginary part of the self-energy

[TABLE]

Here we used

[TABLE]

The physical propagators $G^{\rm ret}_{LL,LR}$ follow from Eq. (92) and using symmetries, e.g., as $G^{\rm ret}_{LR}(\omega)=-\frac{i}{2}[G^{\rm ret}_{+}(\omega)+G^{\rm ret}_{+}(-\omega)^{\ast}]$.

We made substantial progress in understanding the structure of the SD equations for the Maldacena-Qi model and, using symmetries, reduced them to a single propagator $G^{\rm ret}_{+}(\omega)$ and self-energy $\Sigma^{\rm ret}_{+}(\omega)$. Equation (94) is however highly nonlinear and Eq. (93) contains an energy scale $\mu$ which makes the solution difficult to access analytically.

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. (88) with small nonzero $\delta$. The self-energy $\Sigma^{\rm ret}_{+}(\omega)$ is calculated using Eqs. (94) and (95) with help of fast Fourier transform algorithms. $G^{\rm ret}_{+}$ is then computed from Eq. (93) and used to reconstruct the full retarded propagator. New spectral functions are then extracted from Eqs. (86) and the procedure is restarted from these. To improve the convergence properties we follow Ref. Maldacena and Stanford (2016) 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 $\delta$. Results of our numerics are discussed in Sec. III of the main text.