Out-of-time-ordered correlators of the Hubbard model: SYK strange metal in the spin freezing crossover region

# Out-of-time-ordered correlators of the Hubbard model: SYK strange metal in the spin freezing crossover region

Naoto Tsuji RIKEN Center for Emergent Matter Science (CEMS), Wako 351-0198, Japan    Philipp Werner Department of Physics, University of Fribourg, 1700 Fribourg, Switzerland
July 20, 2019
###### Abstract

The Sachdev-Ye-Kitaev (SYK) model describes a strange metal that shows peculiar non-Fermi liquid properties without quasiparticles. It exhibits a maximally chaotic behavior characterized by out-of-time-ordered correlators (OTOCs), and is expected to be a holographic dual to black holes. While a faithful realization of the SYK model in condensed matter systems may be involved, a striking similarity between the SYK model and the Hund-coupling induced spin-freezing crossover in multi-orbital Hubbard models has recently been pointed out. To further explore this connection, we study OTOCs for fermionic single-orbital and multi-orbital Hubbard models, which are prototypical models for strongly correlated electrons in solids. We introduce an imaginary-time four-point correlation function with an appropriate time ordering, which by means of the spectral representation and the out-of-time-order fluctuation-dissipation theorem can be analytically continued to real-time OTOCs. Based on this approach, we numerically evaluate real-time OTOCs for Hubbard models in the thermodynamic limit, using the dynamical mean-field theory in combination with a numerically exact continuous-time Monte Carlo impurity solver. The results for the single-orbital model show that a certain spin-related OTOC captures local moment formation in the vicinity of the metal-insulator transition, while the self-energy does not show SYK-like non-Fermi liquid behavior. On the other hand, for the two- and three-orbital models with nonzero Hund coupling we find that the OTOC exhibits a rapid damping at short times and an approximate power-law decay at longer times in the spin-freezing crossover regime characterized by fluctuating local moments and a non-Fermi liquid self-energy . These results are in a good agreement with the behavior of the SYK model, providing firm evidence for the close relation between the spin-freezing crossover physics of multi-orbital Hubbard models and the SYK strange metal.

## I Introduction

The strange-metal state found by Sachdev and Ye Sachdev and Ye (1993) in a random all-to-all interacting Heisenberg spin model shows a peculiar non-Fermi liquid behavior without quasiparticle excitations, which is reminiscent of the unconventional electronic and magnetic properties of high-temperature superconductors above the superconducting dome. The model exhibits various unusual properties Sachdev and Ye (1993); Georges et al. (2000, 2001), including the absence of magnetic ordering (in the fermionic representation) and a residual entropy down to zero temperature, and an approximate scale invariance at low energies. In particular, the resulting self-energy has a characteristic frequency dependence of at low frequencies. The Sachdev-Ye model has been generalized to a --like model with spins coupled to fermions Parcollet et al. (1998); Parcollet and Georges (1999), which exhibits a wide doping range with similar properties.

Recently, a related fermionic model with random all-to-all interactions and no single-particle hopping has been introduced by Kitaev Kit (); Sachdev (2015); Polchinski and Rosenhaus (2016); Maldacena and Stanford (2016). This model, dubbed Sachdev-Ye-Kitaev (SYK) model, not only retains the non-Fermi liquid properties of the Sachdev-Ye model but also allows for semi-analytic calculations of out-of-time-ordered correlators (OTOCs) Larkin and Ovchinnikov (1969). OTOCs are a novel type of four-point correlation functions such as that ignore the usual time-ordering rule, and are used to describe quantum chaotic properties and information scrambling of quantum many-body systems Shenker and Stanford (2014a, b, 2015); Maldacena et al. (2016); Hosur et al. (2016); Swingle et al. (2016). In particular, the SYK model has been shown to be maximally chaotic in the large- limit Kit (); Maldacena and Stanford (2016), in the sense that OTOCs grow exponentially () at early time with the growth rate Maldacena et al. (2016); Tsuji et al. (2018a) ( is the Boltzmann constant, is the temperature, and is the Planck constant). This property is shared with black holes in Einstein gravity Shenker and Stanford (2014a, b, 2015); Maldacena et al. (2016), which supports expectations that the SYK model may be the holographic dual to gravitational theories. As shown in Ref. Bagrets et al. (2017), the initial exponential growth of OTOCs in the SYK model crosses over to an exponential decay at intermediate times, and eventually to a power-law relaxation in the long-time limit. In the mean time, the SYK model has been generalized to lattice models in which each site represents an SYK atom Gu et al. (2017); Davison et al. (2017); Song et al. (2017); Chowdhury et al. (2018).

While the SYK model shows intriguing universal properties, a question that is of particular interest here is: where can we find a concrete realization of the SYK strange-metal state in condensed matter systems? A faithful implementation of the random all-to-all interaction without single-particle hopping in fermion systems may be involved. Despite its difficulty, there have been several proposals for the realization of the SYK model in condensed-matter systems Danshita et al. (2017); García-Álvarez et al. (2017); Pikulin and Franz (2017); Chen et al. (2018). On the other hand, as pointed out in Ref. Wer (), there is a striking similarity between the SYK strange metal and the spin-freezing crossover regime of multi-orbital Hubbard models with nonzero Hund coupling Werner et al. (2008); Ishida and Liebsch (2010); de’ Medici et al. (2011); Hoshino and Werner (2015), which are prototypical models of strongly correlated electron materials. In these multi-orbital lattice systems, there is a competition between the Hund effect that favors the formation of local magnetic moments and the Kondo effect that screens local magnetic moments. When the two effects are balanced, there emerges a fluctuating moment state with non-Fermi liquid properties. This so-called spin-freezing crossover regime separates a Fermi liquid metal from a spin-moment-frozen metal. A schematic phase diagram of the two-orbital Hubbard model with Hund coupling in the space of the Coulomb interaction and the filling per orbital and spin is shown in Fig. 1, where the ratio between the Hund coupling and is fixed. The crossover from the Fermi liquid to an incoherent metal state with frozen magnetic moments occurs in the region of the doped half-filled Mott insulator. Near the crossover line (red curve), the self-energy shows a non-Fermi-liquid frequency dependence over a significant energy range, and a spin-spin correlation function which decays on the imaginary-time axis as Werner et al. (2008). This is exactly the same non-Fermi liquid behavior as realized in the Sachdev-Ye and SYK models. At low enough temperature, there is a crossover to the Fermi liquid scaling () Stadler et al. (2015), similar to what is found in lattice generalizations of the SYK model Chowdhury et al. (2018).

There are similarities and differences at the level of the Hamiltonian. The local interaction term of multi-orbital Hubbard models with nonzero Hund coupling can be written in the form , where denotes the interaction matrix element, labels the spin and orbital, and () is the electron creation (annihilation) operator. Although the pattern of in realistic systems may be complicated, does not in general resemble a Gaussian random distribution. Furthermore, there are additional single-particle hopping terms in the Hubbard models, which can be a relevant or irrelevant perturbation depending on the strength García-García et al. (2018). Hence it is a nontrivial question whether or not the spin-freezing crossover regime in multi-orbital Hubbard models can be regarded as a realization of the SYK strange metal.

In this work, we calculate out-of-time-ordered correlators for single-orbital and multi-orbital Hubbard models in the thermodynamic limit. Since OTOCs are dynamical (real-time) four-point correlation functions with an unusual ordering of operator sequences, which often requires a formidable effort of summing all the relevant diagrams and solving a Bethe-Salpeter equation, it is numerically challenging to evaluate OTOCs for correlated many-body systems. Previously, the exponential growth of OTOCs has been calculated analytically by summing ladder diagrams for the SYK and several related models Kit (); Maldacena and Stanford (2016); Banerjee and Altman (2017); García-García et al. (2018). For single-particle problems, OTOCs have been numerically calculated for the quantum kicked rotor model Rozenbaum et al. (2017) and the quantum stadium billiard Roz (). OTOCs have also been numerically evaluated for relatively small-size systems using exact diagonalization Che (); Fan et al. (2017); Fu and Sachdev (2016); He and Lu (2017); Huang et al. (2016); Shen et al. (2017); Yao (); Bohrdt et al. (2017); Dóra and Moessner (2017); Dóra et al. (2017). For interacting fermion lattice models, which are typically studied in condensed matter physics, a useful approach is provided by the dynamical mean-field theory (DMFT) Georges et al. (1996), which can be directly applied to the thermodynamic limit, and becomes an exact treatment in the limit of large lattice dimension Metzner and Vollhardt (1989). It involves a self-consistent mapping of the lattice model to an effective quantum impurity model. DMFT has been generalized to deal with nonequilibrium states by switching from the imaginary-time Matsubara formalism to the real-time Kadanoff-Baym or Keldysh formalism defined on the singly-folded time contour Aoki et al. (2014). To calculate OTOCs, nonequilibrium DMFT can be further extended to a doubly folded time-contour formalism, with which OTOCs have been evaluated for the Falicov-Kimball model Tsuji et al. (2017), an integrable lattice fermion model that can be exactly solved within DMFT Freericks and Zlatić (2003). The application of this method to Hubbard models, however, is difficult due to the lack of reliable and versatile impurity solvers that can be used within the doubly folded time-contour formalism.

As an alternative approach, we develop a general method to derive real-time OTOCs from imaginary-time four-point correlation functions through an analytic continuation and the use of the recently formulated out-of-time-order fluctuation-dissipation theorem Tsuji et al. (2018b), in an analogous way as the retarded Green’s function is obtained by analytic continuation from the imaginary-time Matsubara Green’s function. The advantage of this method is that the imaginary-time four-point function can be accurately evaluated numerically for general quantum many-body systems by using an appropriate continuous-time quantum Monte Carlo (QMC) method Gull et al. (2011). This is in contrast to the direct application of QMC methods to the calculation of real-time correlation functions, which usually suffers from a sign problem Werner et al. (2009), that is exacerbated in the OTOC case by the doubling of the real-time contour. We use the proposed method to evaluate OTOCs of single-orbital and multi-orbital Hubbard models within the framework of DMFT using a QMC impurity solver.

The results show that accurate real-time OTOCs can be obtained up to a few hopping times, and also the general long-time behavior can be approximately captured, while the details of the oscillations at intermediate and long times cannot be resolved. For the half-filled single-orbital Hubbard model, we evaluate several different types of OTOCs in the vicinity of the metal-insulator transition. We show that although these functions capture nontrivial correlations in the strongly interacting metallic regime, there is no evidence for the non-Fermi liquid behavior with SYK-like exponents due to the dominance of the Kondo effect. We then turn to the doped Mott insulating phase of the two- and three-orbital Hubbard models with nonzero Hund coupling, for which we focus on a spin-related OTOC that has a counterpart in the SYK model and is sensitive to fluctuating magnetic moments. We find that the OTOC in the spin-freezing crossover regime damps quickly (roughly exponentially) at short times, and decays as a power law at longer times. We confirm that this behavior agrees qualitatively with that of the OTOC for the SYK model with a finite number of orbitals obtained from exact diagonalization. The power-law exponents agree almost quantitatively if we identify the variance of the SYK interaction with the square of the Hund coupling, as suggested in Ref. Wer (). Our analysis of OTOCs thus establishes a close connection between the spin-freezing crossover regime of multi-orbital Hubbard models and the SYK strange metal.

The paper is organized as follows: In Sec. II, we define the imaginary-time four-point correlation functions and discuss their general properties. In Sec. III, we prove that these imaginary-time four-point functions can be analytically continued to real-time OTOCs by introducing the spectral representation and using the out-of-time-order fluctuation-dissipation theorem. In Sec. IV, we present numerical results for OTOCs of the single-, two- and three-orbital Hubbard models, and compare them with the SYK model. Section V contains a summary and conclusions.

## Ii Imaginary-time four-point functions

In this section, we define imaginary-time four-point correlation functions, which we prove in the next section to be analytically continuable to real-time OTOCs, and discuss their general properties. For arbitrary operators and , we define

 CM(AB)2(τ) ≡⎧⎨⎩−⟨^A(τ+βℏ2)^B(βℏ2)^A(τ)^B(0)⟩0≤τ≤βℏ2,−⟨^B(βℏ2)^A(τ+βℏ2)^B(0)^A(τ)⟩−βℏ2≤τ<0. (1)

Here is the inverse temperature, represents the statistical average, is the Hamiltonian of the system, is the partition function, and is the Heisenberg representation for the imaginary-time evolution. We use the label ‘’ for the imaginary-time four-point function because of the analogy with the Matsubara Green’s function defined by

 CMAB(τ) ={−⟨^A(τ)^B(0)⟩0≤τ≤βℏ,∓⟨^B(0)^A(τ)⟩−βℏ≤τ<0, (2)

where the sign is taken when both and are fermionic (i.e., Grassmann odd) and the sign is taken when either or is bosonic (i.e., Grassmann even). In Eq. (2) the function is defined for , while in Eq. (1) the function is defined for . Note that the definition (1) does not have a sign change, independent of the statistical nature (bosonic or fermionic) of and , since the order of and is exchanged twice for in Eq. (1).

An important property of the imaginary-time four-point function is its time periodicity,

 CM(AB)2(τ+βℏ2) =CM(AB)2(τ)(−βℏ2≤τ<0). (3)

The proof for (3) follows straightforwardly from the definition (1):

 CM(AB)2(τ+βℏ2) =−⟨^A(τ+βℏ)^B(βℏ2)^A(τ+βℏ2)^B(0)⟩ =−1ZTr[e−β^H^A(τ+βℏ)^B(βℏ2)^A(τ+βℏ2)^B(0)] =−1ZTr[^A(τ)e−β^H^B(βℏ2)^A(τ+βℏ2)^B(0)] =−⟨^B(βℏ2)^A(τ+βℏ2)^B(0)^A(τ)⟩=CM(AB)2(τ)(−βℏ2≤τ<0). (4)

We can extend the region of the definition of from to by repeatedly applying (3) for (). In this way, can be considered as a periodic function of with period . This allows us to Fourier transform into

 CM(AB)2(iϖn) =∫βℏ20dτeiϖnτCM(AB)2(τ), (5)

where takes the discrete values

 ϖn =4nπβℏ(n∈Z). (6)

Let us compare this situation with the one for the usual imaginary-time two-point function , which is (anti)periodic,

 CMAB(τ+βℏ) =±CMAB(τ), (7)

with period . In Eq. (7), the minus sign is taken when both and are fermionic, and plus otherwise. Due to the (anti)periodicity, one can Fourier transform into

 CMAB(iωn) =∫βℏ0dτeiωnτCMAB(τ) (8)

with the Matsubara frequency given by

 ωn =⎧⎪ ⎪ ⎪⎨⎪ ⎪ ⎪⎩2nπβℏeither % ^A or ^B is bosonic,(2n+1)πβℏboth ^A and ^B % are fermionic. (9)

The period for (1) is half of that for (2). Due to this difference, the frequency step changes between (5) and (8).

## Iii Analytic continuation to real-time OTOCs

In the previous section, we have introduced the imaginary-frequency four-point function (5). Let us recall that the Matsubara Green’s function (8) can be analytically continued to the retarded Green’s function via the replacement . Here is the Fourier transform of

 CRAB(t,t′) ≡−iθ(t−t′)⟨[^A(t),^B(t′)]∓⟩ (10)

with representing the anticommutator () when both and are fermionic and the commutator () otherwise. It is thus natural to ask what kind of function corresponds to the analytic continuation of .

Below we show that the analytic continuation of through is given by what we call the retarded OTOC , which is defined by the Fourier transform of

 CR(AB)2(t,t′) ≡−iθ(t−t′)[⟨^A(t)^B(t′),^A(t)^B(t′)⟩−⟨^B(t′)^A(t),^B(t′)^A(t)⟩]. (11)

Here is the step function defined by () and (), and we used the notation of the bipartite statistical average

 ⟨^X,^Y⟩ ≡Tr(^ρ12^X^ρ12^Y) (12)

(with being the density matrix), which has previously appeared in the study of OTOCs. The original motivation to employ this form was that the squared commutator might be ill-defined in the context of quantum field theory, because two operators can approach each other arbitrarily close in time, which may cause divergences. In this situation, one usually needs to regularize the squared commutator. One prescription to regularize it is to take the bipartite statistical average, , with which the two commutators are separated in the imaginary-time direction. There is an information-theoretic meaning of the difference between the usual and bipartite statistical averages, which is given by the Wigner-Yanase (WY) skew information Wigner and Yanase (1963),

 I12(^ρ,^O) ≡−12Tr([^ρ12,^O]2)=⟨^O2⟩−⟨^O,^O⟩, (13)

for a quantum state and an observable (which is a hermitian operator). It represents the information content of quantum fluctuations of the observable contained in the quantum state (for further details on the WY skew information in the present context, we refer to Refs. Tsuji et al. (2018b); Luo (2005)). Underlying this regularization procedure is the expectation that OTOCs in the form of the usual and bipartite statistical averages would show qualitatively similar characteristic behavior such as the chaotic exponential growth in the short-time regime (butterfly effect). In terms of the bipartite statistical average, the imaginary-time four-point function introduced in the previous section can be written as

 CM(AB)2(τ) =⎧⎨⎩−⟨^A(τ)^B(0),^A(τ)^B(0)⟩0≤τ≤βℏ2,−⟨^B(0)^A(τ),^B(0)^A(τ)⟩−βℏ2≤τ<0. (14)

If we introduce the commutator-anticommutator representation of OTOCs,

 C[A,B]α1[A,B]α2(t,t′) ≡⟨[^A(t),^B(t′)]α1,[^A(t),^B(t′)]α2⟩(α1,α2=±), (15)

can be written in the form

 CR(AB)2(t,t′) =−iθ(t−t′)C{A,B}[A,B](t,t′). (16)

The relation between and is most clearly seen in the spectral representation. To obtain the spectral representation for , we expand it in the basis of eigenstates of denoted by with eigenenergies ,

 CM(AB)2(τ) =−1Z∑klmne−β2(Ek+Em)e1ℏ(Ek−El+Em−En)τ⟨k|^A|l⟩⟨l|^B|m⟩⟨m|^A|n⟩⟨n|^B|k⟩ =−1Z∫∞−∞dω′e−ω′τ∑klmne−β2(Ek+Em)δ(ω′+1ℏ(Ek−El+Em−En)) ×⟨k|^A|l⟩⟨l|^B|m⟩⟨m|^A|n⟩⟨n|^B|k⟩. (17)

From the first to the second line, we inserted , where is the delta function. By Fourier transforming , we obtain

 CM(AB)2(iϖn) =1Z∫∞−∞dω′1−e−βℏω′2iϖn−ω′∑klmne−β2(Ek+Em)δ(ω′+1ℏ(Ek−El+Em−En)) ×⟨k|^A|l⟩⟨l|^B|m⟩⟨m|^A|n⟩⟨n|^B|k⟩ =1Z∫∞−∞dω′1iϖn−ω′∑klmn(e−β2(Ek+Em)−e−β2(El+En))δ(ω′+1ℏ(Ek−El+Em−En)) ×⟨k|^A|l⟩⟨l|^B|m⟩⟨m|^A|n⟩⟨n|^B|k⟩. (18)

Motivated by the above expression, let us define the spectral function for the OTOC by

 A(AB)2(ω) ≡1Z∑klmn(e−β2(Ek+Em)−e−β2(El+En))δ(ω+1ℏ(Ek−El+Em−En)) ×⟨k|^A|l⟩⟨l|^B|m⟩⟨m|^A|n⟩⟨n|^B|k⟩. (19)

Note that takes real values when , since . However, in this case is not necessarily positive semidefinite for . One exception is the low-temperature limit, where becomes positive semidefinite for . To see this, let us denote the ground state as with the eigenenergy . In the zero-temperature limit, the spectral function approaches

 A(AA†)2(ω) →1Z∑lne−βEgδ(ω+1ℏ(2Eg−El−En))⟨g|^A|l⟩⟨l|^A†|g⟩⟨g|^A|n⟩⟨n|^A†|g⟩ −1Z∑kme−βEgδ(ω+1ℏ(Ek+Em−2Eg))⟨k|^A|g⟩⟨g|^A†|m⟩⟨m|^A|g⟩⟨g|^A†|k⟩ =∑km[δ(ω−1ℏ(Ek+Em−2Eg))−δ(ω+1ℏ(Ek+Em−2Eg))]|⟨g|^A|k⟩|2|⟨g|^A|m⟩|2 ≥0(ω≥0). (20)

The spectral sum is given by

 ∫∞−∞dωA(AB)2(ω) =⟨{^A,^B},[^A,^B]⟩=:cAB. (21)

Using the spectral function , the imaginary-frequency function can be written as

 CM(AB)2(iϖn) =∫∞−∞dω′A(AB)2(ω′)iϖn−ω′. (22)

This is analogous to the Lehmann representation for the Matsubara Green’s function,

 CMAB(iωn) =∫∞−∞dω′AAB(ω′)iωn−ω′, (23)

where is the spectral function for the Matsubara Green’s function defined by

 AAB(ω) ≡1Z∑kl(e−βEk∓e−βEl)δ(ω+1ℏ(Ek−El))⟨k|^A|l⟩⟨l|^B|k⟩. (24)

Here the sign is taken when both and are fermionic and the sign is taken otherwise.

In a similar manner, we can obtain the spectral representation of the retarded OTOC, which is expanded in the eigenbasis of the Hamiltonian as

 CR(AB)2(t,t′) =−iθ(t−t′)1Z∑klmne−β2(Ek+Em)[eiℏ(Ek−El+Em−En)(t−t′)⟨k|^A|l⟩⟨l|^B|m⟩⟨m|^A|n⟩⟨n|^B|k⟩ −e−iℏ(Ek−El+Em−En)(t−t′)⟨k|^B|l⟩⟨l|^A|m⟩⟨m|^B|n⟩⟨n|^A|k⟩]. (25)

We permute the summation labels for the second term in Eq. (25) as to obtain

 CR(AB)2(t,t′) =−iθ(t−t′)1Z∑klmn[e−β2(Ek+Em)−e−β2(El+En)]eiℏ(Ek−El+Em−En)(t−t′) ×⟨k|^A|l⟩⟨l|^B|m⟩⟨m|^A|n⟩⟨n|^B|k⟩. (26)

By using the expression for the Fourier transformation of the step function

 θ(t) =i2π∫∞−∞dω′e−iω′tω′+iδ (27)

with a positive infinitesimal constant , we can Fourier transform the retarded OTOC as

 CR(AB)2(ω) =∫∞−∞dω′1ω′+iδ1Z∑klmn[e−β2(Ek+Em)−e−β2(El+En)] ×δ(ω−ω′+1ℏ(Ek−El+Em−En))⟨k|^A|l⟩⟨l|^B|m⟩⟨m|^A|n⟩⟨n|^B|k⟩. (28)

One notices that the same form of the spectral function (19) has appeared in the above expression. Thus, we find that the retarded OTOC has a spectral representation

 CR(AB)2(ω) =∫∞−∞dω′A(AB)2(ω′)ω−ω′+iδ. (29)

One can see that is analytic in the upper half of the complex plane. In the limit of , it behaves as

 CR(AB)2(ω) ∼cABω. (30)

By comparing Eq. (22) and (29), we prove that the imaginary-frequency function can be analytically continued to the retarded OTOC through ,

 CM(AB)2(iϖn) iϖn→ω+iδ−−−−−−→CR(AB)2(ω). (31)

Since is analytic in the upper half plane and uniformly decays to zero as in Eq. (30) for , it should satisfy the Kramers-Kronig relation,

 ReCR(AB)2(ω) =−1πP∫∞−∞dω′ImCR(AB)2(ω′)ω−ω′, (32) ImCR(AB)2(ω) =1πP∫∞−∞dω′ReCR(AB)2(ω′)ω−ω′. (33)

We also define the advanced OTOC as

 CA(AB)2(t,t′) ≡iθ(t′−t)C{A,B},[A,B](t,t′) =iθ(t′−t)[⟨^A(t)^B(t′),^A(t)^B(t′)⟩−⟨^B(t′)^A(t),^B(t′)^A(t)⟩]. (34)

In the same way as for the retarded OTOC, the advanced OTOC has the spectral representation

 CA(AB)2(ω) =∫∞−∞dω′A(AB)2(ω′)ω−ω′−iδ. (35)

Hence the advanced OTOC is analytic in the lower half plane. By comparing Eq. (22) and Eq. (35), we can see that is obtained by analytic continuation from via . The retarded and advanced OTOCs are related via

 CR(AB)2(ω)∗ =CA(B†A†)2(ω). (36)

In the case of , the spectral function (which is real in this case) is given by the imaginary part of the retarded OTOC,

 A(AA†)2(ω) =−1πImCR(AA†)2(ω). (37)

So far, we have explained how to obtain the retarded and advanced OTOCs by analytic continuation of the imaginary-time four-point function . This allows us to access [see Eqs. (11) and (34)]. In order to get the full information on OTOCs, we also need to calculate the complementary part, . This can be done by using the out-of-time-order fluctuation-dissipation theorem, which is the out-of-time-order extension of the conventional fluctuation-dissipation theorem, expressed as

 CKAB(ω) =⎧⎪ ⎪ ⎪⎨⎪ ⎪ ⎪⎩coth(βℏω2)[CRAB(ω)−CAAB(ω)]either ^A or ^B is bosonic,tanh(βℏω2)[CRAB(ω)−CAAB(ω)]both ^A and ^B are fermionic. (38)

Here we have defined the Keldysh Green’s function

 CKAB(ω) =−i⟨[^A(t),^B(t′)]±⟩ (39)

with the sign taken if either or are bosonic and the sign taken if both and are fermionic. Following the analogy between the Green’s functions and OTOCs, let us define the “Keldysh” component of OTOCs as

 CK(AB)2(t,t′) ≡−i2[C{A,B}2(t,t′)+C[A,B]2(t,t′)] =−i[⟨^A(t)^B(t′),^A(t)^B(t′)⟩+⟨^B(t′)^A(t),^B(t′)^A(t)⟩]. (40)

One can see that is exactly the complementary part that we needed to reconstruct OTOCs from the imaginary-time data. The out-of-time-order fluctuation-dissipation theorem has an analogous form to the conventional one,

 CK(AB)2(ω) =coth(βℏω4)[CR(AB)2(ω)−CA(AB)2(ω)]. (41)

Note that the argument of the cotangent factor () is just half of that for the conventional fluctuation-dissipation theorem (38). The out-of-time-order fluctuation-dissipation theorem takes the same form for arbitrary statistics (bosonic or fermionic) for the operators and . In Table 1, we list the definitions and properties of the Green’s function and OTOC. One can see a clear parallelism between the two types of correlation functions.