Quantum chaos in the Brownian SYK model with large finite $N$: OTOCs and tripartite information

TL;DR

This paper investigates quantum chaos and information scrambling in the Brownian SYK model with large finite N, using OTOCs and tripartite information, revealing a logarithmic scrambling time and exponential decay of correlators.

Contribution

It introduces a method to analyze the dynamics of OTOCs and tripartite information in the Brownian SYK model for very large N using a bosonic collective mode approach.

Findings

Scrambling time scales as ln N.

OTOCs and tripartite information decay exponentially after t*.

Numerical calculations performed up to N=10^6.

Abstract

We consider the Brownian SYK model of $N$ interacting Majorana fermions, with random couplings that are taken to vary independently at each time. We study the out-of-time-ordered correlators (OTOCs) of arbitrary observables and the R\'enyi-$2$ tripartite information of the unitary evolution operator, which were proposed as diagnostic tools for quantum chaos and scrambling, respectively. We show that their averaged dynamics can be studied as a quench problem at imaginary times in a model of $N$ qudits, where the Hamiltonian displays site-permutational symmetry. By exploiting a description in terms of bosonic collective modes, we show that for the quantities of interest the dynamics takes place in a subspace of the effective Hilbert space whose dimension grows either linearly or quadratically with $N$, allowing us to perform numerically exact calculations up to $N = 10^6$. We analyze in detail the interesting features of the OTOCs, including their dependence on the chosen observables, and of the tripartite information. We observe explicitly the emergence of a scrambling time $t^\ast\sim \ln N$ controlling the onset of both chaotic and scrambling behavior, after which we characterize the exponential decay of the quantities of interest to the corresponding Haar scrambled values.

Equations

$$H_{\mathrm{SYK}}(t)=i^{q/2}\sum _{1\le i_1<i_2<\dots <i_q\le N}{J}_{i_1,\dots ,i_q}(t){\psi }_{i_1}{\psi }_{i_2}\dots {\psi }_{i_q}.$$

$$\overline{J_{i_1\dots i_q}(t)J_{i_1^{\prime }\dots i_q^{\prime }}\left(t^{\prime }\right)}={\delta }_{i_1{i}_1^{\prime }}\cdots {\delta }_{i_q{i}_q^{\prime }}\delta \left(t-t^{\prime }\right){\sigma }_J\frac{(q-1)!}{N^{q-1}},$$

$$\overline{J_{i_1\dots ,i_4}(t)J_{i_1^{\prime }\dots ,i_4^{\prime }}\left(t^{\prime }\right)}={\delta }_{i_1{i}_1^{\prime }}\cdots {\delta }_{i_4{i}_4^{\prime }}\delta \left(t-t^{\prime }\right)\frac{1}{N^3}.$$

$${\mathcal{F}}_{\mathcal{O},{\mathcal{O}}^{\prime }}(t)=\mathrm{tr}\left[\rho \mathcal{O}(t){\mathcal{O}}^{\prime }(0)\mathcal{O}(t){\mathcal{O}}^{\prime }(0)\right],$$

$$\mathcal{C}(t)=\frac{1}{2}\mathrm{tr}\left[\rho {\left(\left\{{\psi }_j(t),{\psi }_k(0)\right\}\right)}^{†}\left(\left\{{\psi }_j(t),{\psi }_k(0)\right\}\right)\right],$$

$$\mathcal{C}(t)=1+\mathrm{Re}\left[\mathrm{tr}\left(\rho {\psi }_j(t){\psi }_k(0){\psi }_j(t){\psi }_k(0)\right)\right].$$

$$|I⟩=\frac{1}{D^{N/2}}\sum _{j=1}^{D^N}|j⟩\otimes |j^{\prime }⟩,$$

$$U\mapsto |U⟩⟩=\left(1\otimes U\right)|I⟩.$$

$$I_3(A:C:D)=I(A:C)+I(A:D)-I(A:CD),$$

$$I(X:Y)=S_X+S_Y-S_{XY},$$

$$S_{AC}=-\mathrm{tr}\left[{\rho }_{AC}\ln {\rho }_{AC}\right],$$

$$\left({\psi }_j^a-i{\psi }_j^b\right)|I^{ab}⟩=0,\forall j.$$

$$U(t)\mapsto |U(t)⟩⟩=U^a(t)|I^{ab}⟩.$$

$$I_3^{(2)}(A:C:D)=I^{(2)}(A:C)+I^{(2)}(A:D)-I^{(2)}(A:CD),$$

$$I^{(2)}(X:Y)=S_X^{(2)}+S_Y^{(2)}-S_{XY}^{(2)},$$

$$S_X^{(2)}=-\ln \left[\overline{\mathrm{tr}\left({\rho }_X^2\right)}\right].$$

$$I_3^{(2)}(A:C:D)=\frac{N}{2}\ln (2)-S_{AC}^{(2)}-S_{AD}^{(2)},$$

$${\mathcal{F}}_{\mathcal{O},{\mathcal{O}}^{\prime }}(t)=\frac{1}{2^{N/2}}\mathrm{tr}\left\{\mathcal{O}U(t){\mathcal{O}}^{\prime }{U}^{†}(t)\mathcal{O}U(t){\mathcal{O}}^{\prime }{U}^{†}(t)\right\}.$$

$$U(t)=e^{-i\Delta t{H}_{\mathrm{SYK}}(t_n)}{e}^{-i\Delta t{H}_{\mathrm{SYK}}(t_{n-1})}\dots e^{-i\Delta t{H}_{\mathrm{SYK}}(t_1)},$$

$$\overline{J_{i_1\dots ,i_4}(t_r)J_{i_1^{\prime }\dots ,i_4^{\prime }}\left(t_s\right)}={\delta }_{i_1{i}_1^{\prime }}\cdots {\delta }_{i_4{i}_4^{\prime }}\frac{1}{N^3}\frac{{\delta }_{r,s}}{\Delta t}.$$

$${\mathcal{F}}_{\mathcal{O},{\mathcal{O}}^{\prime }}$$

$$\times$$

$${\mathcal{F}}_{\mathcal{O},{\mathcal{O}}^{\prime }}=⟨L|\left(U\otimes U^{\ast }\otimes U\otimes U^{\ast }\right)|R⟩,$$

$$|L⟩=$$

$$|R⟩=$$

$${\chi }_j^a:=$$

$${\chi }_j^c:=$$

$${\mathcal{Q}}^{\alpha }=\prod _{k=1}^N{\chi }_k^{\alpha },\alpha =a,b,c,d.$$

$${\psi }_j^a=$$

$${\psi }_j^c=$$

$$\mathcal{U}(t)=U_{+}^a(t)U_{-}^b(t)U_{+}^c(t)U_{-}^d(t),$$

$$U_{\pm }^{\alpha }(t)=e^{\mp i\Delta t{H}_{\mathrm{SYK}}^{\alpha }(t_n)}{e}^{\mp i\Delta t{H}_{\mathrm{SYK}}^{\alpha }(t_{n-1})}\dots e^{\mp i\Delta t{H}_{\mathrm{SYK}}^{\alpha }(t_1)},$$

Full text

1. I Introduction

2. II The model and the chaos quantifiers

3. II.1 The OTOCs and the operator spreading

4. II.2 Diagnostic of scrambling: the tripartite information in fermionic systems

5. III Exact approach from emergent permutational symmetry

6. III.1 Decomposing the dynamical problem

7. III.2 The generator of the dynamics: mapping to a bosonic system

8. III.3 The OTOC and the tripartite information

9. IV The physical results

10. IV.1 The OTOCs: numerical results

11. IV.2 The Rényi-$2$ tripartite information

12. V Deriving the key formulas

13. V.1 The Hamiltonian

14. V.2 Extracting OTOCs and Entropies

15. V.2.1 The OTOCs

16. V.2.2 The Rényi-$2$ entanglement entropy $S^{(2)}_{AC}({\bar{\ell}})$

17. V.3 Some large-$N$ limits

18. VI Conclusions

19. A Non-interacting case: $q=2$

20. B Relation between OTOCs and Rényi-$2$ entropies

21. B.1 The case of Pauli matrices

22. B.2 The case of Majorana operators

23. C Details on the numerical implementation

24. D The Rényi-$2$ entanglement entropy $S^{(2)}_{AD}({\bar{\ell}})$

I Introduction

The study of many body quantum chaos is currently experiencing a golden age, also due to its implications on important aspects in many-body physics such as the thermalization Deut91 ; Sred94 of isolated systems RiDO08 ; HoQi16 ; DKPR16 ; LeQi18 ; HaMa13 ; LiSu14 , or the scrambling of quantum information HQRY16 ; LFSL19 ; BeGL19 . In fact, the field has already enjoyed intense research activity more than thirty years ago Gutzwiller_book ; stockmann_book , when the relations between chaotic many-body systems and random matrix theory were first explored. Recently, a renewed interest came from the study of black hole physics and concepts such as scrambling of quantum information, and computational complexity HaPr07 ; LSHO13 ; SeSu08 ; ShSt14 ; KitaevTalk2 ; Suss16 ; BRSS16 .

An important milestone in the recent literature is the Sachdev-Ye-Kitaev model SaYe93 ; Kitaev_talk , which was originally proposed by Sachdev and Ye as a model of strongly correlated electron systems, and generalized by Kitaev in 2015 who pointed out its connection to holographic duality. This model describes $N$ Majorana fermions or complex fermions with random all-to-all interactions. In the work we will focus on the Majorana fermion version. The SYK model has already drawn enormous attention from different communities, ranging from quantum gravity MaSt16 ; MaSY16 to condensed-matter and many-body physics PoRo16 ; BaAK16 ; DFGG17 ; KlTa17 ; GuQS17 ; SoJB17 ; CWBS18 , due to the concomitance of several unique features. Among these, the model has been shown to be maximally chaotic and yet amenable to exact analysis in the large-$N$ limit Kitaev_talk ; PoRo16 ; MaSS16 ; BaAK16 ; BaAK17 , making it an ideal playground for the study of chaos and scrambling of quantum information.

In the same years, the effort to better characterize quantum chaos led to the systematic development of reliable indicators for its diagnosis. In particular, out-of-time-ordered correlation (OTOC) functions, historically introduced in the context of disordered superconductors LaOv69 , were naturally selected as ideal probes to detect the “scrambling” of local observables ShSt14_bh ; ShSt14 ; Kitaev_talk ; RoSS15 ; RoSS18 , namely the spreading of their spatial support in the operator basis. It is important to mention that these ideas had far reaching ramifications, motivating the study of OTOCs also in many-body systems with short-range interactions RoSw16 ; TsWU17 ; AlFI16 ; SwCh17 ; Yung17 ; PaSa17 ; KuGP17 ; GuQS17 ; DoMo17 ; LiMo18 ; LiMo18 ; SKMK18 ; NIDS19 ; McNK19 ; HuBZ19 ; CLBS19 and in spatially local “quantum unitary circuits” NRVH17 ; SPHS18 ; NaVH18 ; KRPS18 ; ChDC18 ; RaPK18 ; KhVH18 ; KoLP18 ; BeKP18 ; Hunt19 ; GuHu19 ; ZhNa19 ; FCDC19 , which provide minimally structured models for chaotic quantum dynamics. In fact, related studies on information scrambling in a class of random, and in general non-local, circuits (the so-called approximate $t$-designs) were already carried out within quantum information theory EmLL05 ; EWSL03 ; DaOP07 ; GrAE07 ; OlDP07 ; Znid07 ; Znid08 ; ArBr08 ; HaLo09 ; BrVi10 ; DiJo11 ; BrFa12 ; BrHH16 ; NHKW17 ; SoCBQA19 , where the latter were used to provide rapid approximations to random unitary operators. Finally, we note that OTOCs were also shown to be directly related to the growth of the operator size, i.e. the size of its support RoSS18 ; QiSt18 .

So far, computations of OTOCs in the SYK model have been carried out through field-theoretical approaches in the large-$N$ limit Kitaev_talk ; MaSt16 ; PoRo16 . On the other hand, despite the many works devoted to this topic, results for finite values of $N$ are difficult to obtain, and remain scarce GaVe16 ; GaVe17 ; GaJV18 . This is also true for numerical computations: the exponential growth of the Hilbert space dimension, and the presence of disorder averages yield strong limitations on the sizes of the systems that can be simulated FuSa16 ; RoSS18 ; CGHP17 ; GuMV18 . Still, it would be highly desirable to develop a systematic approach to investigate the properties of the SYK model at finite $N$, even numerically. Indeed, not only would this allow for inspection of finite-size corrections to the large-$N$ results, but also to compute quantities beyond multi-point correlation functions, for which field-theoretical approaches might be difficult to apply. A notable example is given by the (negative) tripartite information of the evolution operator introduced in Ref. HQRY16 in the context of unitary circuits. This was suggested as a valuable tool to quantify the scrambling power of a quantum channel, namely its ability to delocalize information provided as an input. We note that, so far, this quantity has been computed only numerically for small system sizes HQRY16 ; SPMK18 (see also Refs. IySa18 ; PRZI18 ; SeML18 , where the tripartite information of given states, and not of the channel, was studied).

Motivated by the above picture, we consider a simpler, but closely related, Brownian SYK model, and address the problem of its exact analysis at finite sizes. The model was introduced in Ref. SaSS18 , and differs from the traditional SYK in that the random couplings are chosen to vary independently at each time. The simplification arising in this case is similar to the one we have in unitary circuits by choosing random gates independently in each layer. Experience from the latter framework suggests that the main features of the chaotic dynamics remain qualitatively unaltered by introducing an additional time-dependence to the spatial disorder, except that random circuits and Brownian models behave like infinite-temperature systems since they do not display energy conservation.

In this work, we focus on the development of a systematic approach to the chaotic dynamics in the Brownian SYK model, which could also be applied, more generally, to other time-dependent, disordered Hamiltonians with infinite-range interactions. In particular, we aim to compute OTOCs of arbitrary local observables, and other dynamical quantities which can be extracted from disordered averages involving up to four unitary evolution operators. These include a Rényi-$2$ version of the tripartite information introduced in HQRY16 , which has been shown to encode information about all possible OTOCs  HQRY16 .

As a main result of our work, we show that the averaged dynamics of the OTOCs and of the Rényi tripartite information can be studied as a quench problem at imaginary times in a model of $N$ qudits, where the Hamiltonian displays full site-permutational symmetry. We analyze this problem by means of a description in terms of bosonic collective modes, and prove that for the quantities of interest the dynamics takes place in a subspace of the Hilbert space whose dimension grows either linearly or quadratically with $N$. This allows us to perform numerically exact calculations up to one million particles, and, consequently, analyze in great detail the behavior of OTOCs and of the Rényi tripartite information, highlighting their most interesting features. While some of our results depend on simplifications arising in the special case of the SYK model, we expect that suitable generalizations of our method could be successfully applied also to the study of other disordered time-dependent Hamiltonians with all-to-all interactions.

It is useful to compare our method with that of existing studies, as some of the ideas used in our work are related to other approaches in the literature. First, Ref. SaSS18 proposed the Brownian SYK model as a simplified version of the original SYK, and mainly focused on the computation of the spectral form factor Haake_book . For this specific quantity, it was shown that in the Brownian SYK model an exact solution could be achieved, by means of an elementary mapping to a classical partition function. Our results on OTOCs and tripartite information cannot be obtained using the same approach.

Next, we discuss Refs. Znid07 ; Znid08 ; BrVi10 , where a class of random quantum circuits was considered, in which at each layer a single unitary gate is applied to a pair of qudits randomly chosen. There, it was shown that the moments of the evolution operator associated with a time step could be mapped onto a permutational invariant Hamiltonian which generalizes the Lipkin-Meshkov-Glick model RiVM08 . Even though the idea underlying our method is similar, both our mapping and the quantities studied in this paper are different.

We also note that the computation of OTOCs in models with a continuous-time evolution in the presence of Brownian disorder and infinite-range interactions have been already addressed in ShSt15 ; ZhCh19 (see also ChZh18 ; XuSw18 ). The system studied in these works, consisting of $N$ qudits driven by an Hamiltonian which is bilinear in the Pauli operators, was introduced as a chaotic toy model in Ref. LSHO13 , where its scrambling time was first estimated to be logarithmic in $N$ (see also Ref. GHST18 , where the spectral form factor was analyzed for the same system). The approach of ShSt15 ; ZhCh19 relies on the derivation, based on an application of Itô calculus Parthasarathy_book , of a system of differential equations for the OTOCs of interest. Solving the latter, numerical results were given in Ref. LSHO13 for sizes comparable to those that can be reached with our method, while an analytical solution was found in  ZhCh19 for a particular average of OTOCs. As we will see, our approach differs from that of ShSt15 ; ZhCh19 , as we tackle directly the computation of the averaged moments of the evolution operator. This allows us to use the same formalism to also analyze the tripartite information discussed above, which was not addressed in these studies. Finally, we note that rigorous results, relevant to the present paper, for the scrambling properties of continuous-time evolution generated by random Hamiltonians were recently presented in Refs. BaBK17 ; OBKB17 .

The organization of the rest of this paper is as follows. In Sec. II we introduce the Brownian SYK model and the quantities which will be investigated in this work. We proceed to present the key features of our method in Sec. III, while our physical results are reported in Sec. IV. The most technical aspects of our study are consigned to Sec. V and to a few appendices. Finally, our conclusions are discussed in Sec. VI.

II The model and the chaos quantifiers

The object of study of this work will be the Brownian SYK model, describing a set of $N$ Majorana fermions with $q$-local, all-to-all random interactions, cf. Fig. 1. It is defined on a Hilbert space $\mathcal{H}_{N}$ of dimension $D=2^{N/2}$, with $N$ operators $\psi_{j}$ acting on $\mathcal{H}_{N}$. They are the representation of standard Majorana fermions, and thus satisfy $\{\psi_{j},\psi_{k}\}=2\delta_{j,k}$ and $\psi_{j}^{\dagger}=\psi_{j}$ (the quantities of interest in this work will not depend on the representation chosen for the $N$ Majorana fermions). Its time-dependent Hamiltonian reads

[TABLE]

Here, the couplings $J_{i_{1},\ldots\,,i_{q}}(t)$ are random variables, which we assume to be Gaussian distributed with vanishing mean and variance

[TABLE]

where we denoted by $\overline{[\ldots]}$ the average over disorder realizations. While our method could be applied for arbitrary integer values of $q$, we will focus for concreteness on the case $q=4$. Furthermore, we will choose the constant $\sigma_{J}$ in such a way that

[TABLE]

In comparison, the original SYK Hamiltonian shares the same form of (1), but with time-independent couplings. In appendix A we additionally discuss the case $q=2$, which lacks chaotic behavior as each disorder realization is non-interacting.

II.1 The OTOCs and the operator spreading

As we have already discussed in Sec. I, we will be mainly interested in two quantifiers of quantum chaos and scrambling. The first one is given by OTOCs of local observables: explicitly, given two operators $\mathcal{O}$, $\mathcal{O}^{\prime}$, we define their OTOC on a state $\rho$ as

[TABLE]

where $\mathcal{O}(t)=U^{\dagger}(t)\mathcal{O}U(t)$, and $U(t)$ is the unitary evolution operator. In this work we will choose the infinite-temperature Gibbs density matrix $\rho=\mathbb{1}/2^{N/2}$, which represents a stationary state for the time-dependent Hamiltonian (1).

Importantly, we recall that the OTOC (4) can be related to an intuitive notion of the spreading of localized operators under unitary evolution. To this end, we choose for simplicity $\mathcal{O}=\psi_{j}$, $\mathcal{O}^{\prime}=\psi_{k}$ with $j\neq k$, and consider the quantity

[TABLE]

which measures the magnitude of the anticommutator between $\psi_{j}(t)$ and $\psi_{k}(0)$. At time $t=0$, one simply has $\mathcal{C}(t)=0$. On the other hand, as time increases, the spacial support of $\psi_{j}(t)$ will also increase; namely $\psi_{j}(t)$ will evolve into a complicated sum of strings of Majorana operators. Then, we see that deviations of $\mathcal{C}(t)$ from zero signal that the support of $\psi_{j}(t)$ has grown to include site $k$. Accordingly, $\mathcal{C}(t)$ can be understood as a measure of the spatial spreading of the local operator $\psi_{j}(t)$. The connection between the latter and OTOCs is finally established by the simple relation

[TABLE]

In conclusion, the discussion above allows one to view the OTOCs as a measure of chaos: chaotic dynamics corresponds to OTOCs that vanish sufficiently rapidly with time. On the other hand, for a non-chaotic Hamiltonian one expects information to spread coherently: for large system sizes this results in either a slow decay or a non-vanishing asymptotics of OTOCs LiMo18 , while for small ones this causes revivals, consisting in OTOCs frequently returning close to their original value GoYD19 .

II.2 Diagnostic of scrambling: the tripartite information in fermionic systems

The OTOCs provide a physically clear definition of quantum chaos in terms of correlation functions between local operators. Other measures probing different features intuitively associated with chaos exist. Among these, the notion of scrambling of information, originally introduced in the study of black hole physics HaPr07 ; SeSu08 , is particularly clear: a quantum system is a good scrambler if a localized perturbation in the initial state spreads over all its degrees of freedom, in such a way that it can no longer be detected by local measurements at large times. In this context, it is useful to think of the unitary evolution as a quantum channel, taking an initial state as the input, and returning the evolved state as the output. In this logic, it was proposed in Ref. HQRY16 that the scrambling power of a channel could be conveniently measured by the tripartite information between bipartitions of its input and output, as we review in the following.

For simplicity, let us first consider a system of $N$ qudits, associated with a Hilbert space $\mathcal{H}_{N}=h_{1}\otimes\ldots\otimes h_{N}$, where $h_{j}\simeq\mathbb{C}^{D}$, and a unitary operator $U:\mathcal{H}_{N}\to\mathcal{H}_{N}$. In order to study the scrambling properties of $U$, we wish to interpret it as a state in a suitable space. To this end, we introduce a copy of the original Hilbert space $\mathcal{H}^{\prime}_{N}$, and define the maximally entangled state $\ket{I}\in\mathcal{H}_{N}\otimes\mathcal{H}^{\prime}_{N}$ as

[TABLE]

where $\{\ket{j}\}_{j=1}^{D^{N}}$, $\{\ket{j^{\prime}}\}_{j^{\prime}=1}^{D^{N}}$ are orthonormal bases for $\mathcal{H}_{N}$ and $\mathcal{H}^{\prime}_{N}$, respectively. Note that we choose the basis such that $\ket{I}$ is a direct product of EPR pairs between qudits in the two systems, as is illustrated in Fig. 2. Then, the operator $U$ can be interpreted as a state in $\mathcal{H}_{N}\otimes\mathcal{H}_{N}^{\prime}$ through the Choi-Jamiolkowski mapping

[TABLE]

Here the operator $U$ is depicted as a box, whose legs correspond to the local Hilbert spaces $h_{j}$; we see that one could intuitively think of the state $\ket{U}\rangle$ as obtained by “bending” the output legs, so as to treat input and output, associated with $\mathcal{H}_{N}$ and $\mathcal{H}^{\prime}_{N}$ respectively, on an equal footing. It should be noted that the mapping from $U$ to $\ket{U}\rangle$ is not unique, as it depends on the choice of state $\ket{I}$. However, different $\ket{I}$ are related by a local unitary transformation, which does not affect the entropy-related quantities we discuss in the following.

Given $\ket{U}\rangle\in\mathcal{H}_{N}\otimes\mathcal{H}^{\prime}_{N}$, one can compute the entanglement entropy between different spatial regions in $\mathcal{H}_{N}$ and $\mathcal{H}^{\prime}_{N}$. We consider in particular bipartitions of $\mathcal{H}_{N}$ and $\mathcal{H}^{\prime}_{N}$ into the complementary subsystems $A$, $B$ and $C$,$D$ respectively; in Fig. 2 a special choice for these regions is shown. Given a pair of bipartitions $(A,B)$ and $(C,D)$, we define the tripartite information as HQRY16

[TABLE]

where $CD$ denotes the union of the regions $C$ and $D$. Here $I(X:Y)$ is the mutual information between the regions $X$ and $Y$

[TABLE]

where $S_{X}$ is the von Neumann entropy of the reduced density matrix $\rho_{X}$. For instance, we have

[TABLE]

where $\rho_{AC}=\operatorname{tr}_{BD}[\rho]$.

The tripartite information in Eq. (9) was suggested in Ref. HQRY16 as a natural and convenient diagnostic for scrambling. In fact, as in the case of OTOCs, its underlying physical meaning is easy to grasp. From Eq. (9), we see that $-I_{3}(A:C:D)$ quantifies the amount of information on the input region $A$ that can be recovered by global measurements in $C\cup D$, but can not be obtained by probing $C$ and $D$ individually. Recalling that $\mathcal{H}^{\prime}_{N}=C\cup D$ corresponds to the output, this is exactly a measure of scrambling: if $-I_{3}(A:C:D)$ is large, it means that the information localized in a subsystem $A$ of the input state can be recovered only by global measurements in the output state, and information has been scrambled. Accordingly, if for any bipartition of $\mathcal{H}_{N}$ and $\mathcal{H}^{\prime}_{N}$, $I_{3}(A:C:D)$ is negative with an absolute value close to the maximum possible value, the channel $U$ has large scrambling power. Finally, a close connection was established in Ref. HQRY16 between the tripartite information (9) and the OTOCs, which further corroborated the appeal of the former as a valuable diagnostic of scrambling and, more generally, of chaotic dynamics. This connection is reviewed in Appendix B, where we also discuss its generalization to the fermionic setting.

The above discussion is carried out in terms of qudits, whereas in our work we are interested in a fermionic system. At this point, one could employ a Jordan-Wigner representation of the Majorana operators in the Hamiltonian (1), interpret the resulting evolution operator as a unitary channel acting on a system of $N/2$ qubits, and define the tripartite information for the latter according to the discussion above. However, given a correspondence between Majorana and Pauli operators via the Jordan-Wigner transformation, it is known that the reduced density matrix of disjoint intervals written in terms of the two is not the same, leading to different results for the corresponding von Neumann entanglement entropy IgPe10 ; FaCa10 . In our case, we stress that the physical degrees of freedom are represented by the Majorana operators and, accordingly, the tripartite information in Eq. (9) should be computed in terms of the latter. In this respect, we find it useful to discuss explicitly the generalization of the above construction for Majorana operators, without making direct reference to the tensor-product structure of the doubled Hilbert space associated with the input and output of the channel.

As a first ingredient, we wish to interpret the evolution operator generated by the Hamiltonian (1) as a state. To this end, we introduce a system of $2N$ Majorana operators $\psi^{\alpha}_{j}$, where $j=1\,,\ldots\,,N$, while $\alpha$ is an index labeling two different species which we denote by $a$ and $b$. The maximally entangled state $\ket{I^{ab}}$ is then defined as the vacuum state for the complex fermions $c_{j}=\psi_{j}^{a}-i\psi_{j}^{b}$ GuLQ17 , namely

[TABLE]

The operator $U$ can now be interpreted as a state in the doubled system through the mapping

[TABLE]

Here the superscript $a$ indicates that the Hamiltonian generating the unitary evolution operator $U^{a}(t)$ is written in terms of the fermions $\psi^{a}_{j}$. A pictorial representation of this construction is shown in Fig. 3. One can now proceed to compute the fermionic reduced density matrices for the evolved state $\ket{U(t)}\rangle$, and consequently the corresponding tripartite information as in Eq. (9). We refer the reader to Sec. V for further details.

Unfortunately, despite its great interest, the computation of the tripartite information (9) is a very difficult task, which so far has been carried out only numerically for qudit systems of small sizes HoQi16 ; SPMK18 . For this reason, we study a simpler but closely related quantity, which is obtained from $I_{3}(A:C:D)$ by considering Rényi, rather than von Neumann, entropies. Specifically, we will compute the following Rényi-$2$ tripartite information

[TABLE]

where

[TABLE]

and

[TABLE]

We note that, strictly speaking, $S^{(2)}_{X}$ is not the averaged Rényi entropy of order $2$, as the disorder average is taken inside the logarithm. However, Ref. HoQi16 showed that the OTOC for a pair of operators in $A$ and $C$ averaged over all operator choices is determined by ${\rm tr}\left(\rho_{AD}^{2}\right)$. Therefore the averaged purity $\overline{{\rm tr}\left(\rho_{X}^{2}\right)}$ is a meaningful physical quantity to consider. Also, for $N$ not too small, one expects the effect of fluctuations in the disorder to be small, so that $S^{(2)}_{X}$ remains a good approximation for the Rényi-$2$ entropy HoQi16 ; KRPS18 .

It is worth to notice that Eq. (14) can be simplified in general. Indeed, it is easy to show HQRY16

[TABLE]

where we used that the dimension of the Hilbert space associated with $N$ Majorana fermions is $D=2^{N/2}$. Eq. (17) tells us that, in order to obtain the tripartite information, it is sufficient to compute the entropies $S^{(2)}_{AC}$ and $S^{(2)}_{AD}$ between different regions of the input and the output.

We conclude this section by stressing that while the Rényi tripartite information (14) differs quantitatively from $I_{3}(A:C:D)$, based on previous studies GuLQ17 , we can still expect it to display the same qualitative features of the latter, and thus to be a suitable measure for scrambling.

III Exact approach from emergent permutational symmetry

Having introduced the model and the quantities of interest, we proceed by presenting the general ideas of the method developed in this work. The physical results will be then discussed in Sec. IV, while we postpone the most technical details of our calculations to Sec. V.

III.1 Decomposing the dynamical problem

We will begin our discussion with the concrete problem of computing the OTOC (4), which we rewrite as

[TABLE]

We recall that the time-dependent, disordered Hamiltonian (1) gives rise to a dynamics which can be interpreted as the continuous limit of the discrete process defined by

[TABLE]

where $\Delta t=t/n$ and $t_{j}=j\Delta t$, while the delta function in Eq. (3) is regularized as

[TABLE]

In practice, one can work with the discrete form (19) of the evolution operator, and take the continuum limit at the end of the calculations.

In order to compute $\mathcal{F}_{\mathcal{O},\mathcal{O}^{\prime}}(t)$, we first introduce a resolution of the identity between each pair of operators in  (18), yielding

[TABLE]

Here $\{\ket{j}\}$ denotes a basis for the Hilbert space $\mathcal{H}_{N}$ [introduced before Eq. (1)] on which the operators $\psi_{j}$ act. Rearranging the above sum, we obtain

[TABLE]

where

[TABLE]

Here $U^{\ast}(t)$ denotes the complex conjugate of $U(t)$ (which is well defined, once a basis $\{\ket{j}\}$ of $\mathcal{H}_{N}$ is given) and we introduced the vectors $\ket{i,j,k,l}=\ket{i}\otimes\ket{j}\otimes\ket{k}\otimes\ket{l}\in\mathcal{H}_{N}^{\otimes 4}$. According to Eq. (22), the dynamical information about the OTOC is uniquely encoded in the operator $\mathcal{U}(t)\equiv U\otimes U^{\ast}\otimes U\otimes U^{\ast}$, while $\mathcal{O}$, $\mathcal{O}^{\prime}$ only affect the “left” and “right” states $\ket{L}$, $\ket{R}$, cf. Fig. 4 .

From Eq. (19), we see immediately that $\mathcal{U}(t)$ is written in terms of the operators

[TABLE]

which provide a basis for all the operators in $\mathcal{H}^{\otimes 4}_{N}$. Note that, as we already stressed, $\psi_{j}$ is the representation of a Majorana fermion, and thus is an operator acting on $\mathcal{H}_{N}$, for which the tensor product is defined in the usual way. Due to the tensor-product structure of $\mathcal{H}_{N}^{\otimes 4}$, the operator $\chi_{j}^{\alpha}$ satisfy mixed commutation and anticommutation relations, namely $\left[\chi^{\alpha}_{j},\chi^{\beta}_{k}\right]=0$ if $\alpha\neq\beta$, while $\left\{\chi^{\alpha}_{j},\chi^{\alpha}_{k}\right\}=2\delta_{j,k}$. On the other hand, it is possible to introduce related operators in $\mathcal{H}_{N}^{\otimes 4}$ which are all anti-commuting with one another, realizing a truly fermionic algebra. We consider for concreteness the case $N\equiv 0\ ({\rm mod}\ 4)$ [if $N\equiv 2\ ({\rm mod}\ 4)$, one has a similar treatment], and introduce

[TABLE]

Then, we can define

[TABLE]

One can easily verify that $\{\psi^{\alpha}_{j}\}_{j,\alpha}$ satisfy fermionic anticommutation relations, namely $\{\psi_{j}^{\alpha},\psi_{k}^{\beta}\}=2\delta_{\alpha,\beta}\delta_{j,k}$, and that $\left(\psi^{\alpha}_{j}\right)^{\dagger}=\psi^{\alpha}_{j}$. Furthermore, since $\left(\mathcal{Q^{\alpha}}\right)^{2}=\mathbb{1}$, we have $\prod_{k=1}^{M}\chi_{j_{k}}^{\alpha}=\prod_{k=1}^{M}\psi_{j_{k}}^{\alpha}$ for any even integer $M$. Since the Hamiltonian (1) contains a sum of products of Majorana operators with an even number of particles, it is then straightforward to show

[TABLE]

where

[TABLE]

while $H_{\rm SYK}^{\alpha}$ is the Hamiltonian (1) written in terms of the fermions $\psi^{\alpha}_{j}$. We see that $\mathcal{U}(t)$ can be viewed as an evolution operator on the space of four “replica” Majorana fermions $\psi^{\alpha}_{j}$, labeled by $\alpha=a$, $b$, $c$, $d$. Eq. (30) represents the starting point for our subsequent calculations.

The above discussion allows us to decompose the problem of computing the OTOC (18) into two logically separated steps:

• •

compute the disorder average of the generator of the dynamics $\mathcal{U}(t)$, defined in Eq. (30) (cf. Sec. III.2);

• •

given the operator $\overline{\mathcal{U}}(t)$, evaluate the matrix element $\braket{L}{\overline{\mathcal{U}}(t)}{R}$, where $\ket{L}$, $\ket{R}$ were introduced in Eq. (22) and pictorially represented in Fig. 4 (cf. Sec. III.3).

Importantly, it is possible to show that the same procedure can be employed for the Rényi-$2$ tripartite information (14): one can express also this quantity in the form of a matrix element $\braket{L}{\overline{\mathcal{U}}(t)}{R}$, for an appropriate choice of the vectors $\ket{L}$ and $\ket{R}$, cf. Sec. III.3. We will address the two points above separately in the following subsections, for both the OTOCs and the tripartite information, providing a complete overview of the approach developed in this work.

III.2 The generator of the dynamics:

mapping to a bosonic system

We start by addressing the computation of the average evolution operator defined in Eq. (30). Using that even numbers of different Majorana operators commute, and that one can factor disorder averages at different times, we note that Eqs. (30), (31) imply

[TABLE]

This allows us to write down a linear differential equation for $\overline{\mathcal{U}}(t)$, as follows.

First, from Eq. (20), we see that, in order to expand the first line at the first order in $\Delta t$, each exponential factor has to be expanded up to the second order. By doing this, and carefully taking into account the correlations between the couplings, one obtains an equation of the form

[TABLE]

namely, taking the limit $\Delta\to 0$

[TABLE]

where

[TABLE]

Here, the indexes $a,b,c,d$ are ordered as $a<b<c<d$, while we introduced

[TABLE]

We note that, since the disorder average has been already taken, the operator $L$ is time- and disorder-independent. Eq. (34) can thus be seen as a Schrodinger equation (at imaginary times) for $\overline{\mathcal{U}}(t)$ in the space ${\rm End}(\mathcal{H}_{N}^{\otimes 4})$ of the linear endomorphisms acting on $\mathcal{H}_{N}^{\otimes 4}$, where the left matrix multiplication by $L$ is interpreted as a superoperator. In the following, it will be useful to denote by $\ket{\mathcal{O}}\rangle$ the state in ${\rm End}(\mathcal{H}_{N}^{\otimes 4})$ associated with the operator $\mathcal{O}$.

In order to proceed, we note that at every time $t$, the operator $\overline{\mathcal{U}}(t)$ can be written as a linear superposition of operators of the form

[TABLE]

where $\mathcal{O}_{j}^{\alpha_{j}}$ is chosen within the set of operators

[TABLE]

Indeed, due to the anticommutation relations of the Majorana operators and the form of the Hamiltonian $H$, it is easy to see that $\overline{\mathcal{U}}(t)$ can not contain terms with an odd number of fermions at site $\psi^{\alpha}_{j}$. Hence, the dynamics of $\ket{\overline{\mathcal{U}}(t)}\rangle$ takes places in the Hilbert space generated by the vectors

[TABLE]

Here, $\alpha_{j}\in\{1,ab,ac,ad,bc,bd,cd,abcd\}$, with the convention $\mathcal{O}_{j}^{1}=\mathbb{1}$, $\mathcal{O}_{j}^{ab}=\psi^{a}_{j}\psi^{b}_{j}$, $\ldots$, $\mathcal{O}_{j}^{abcd}=\psi^{a}_{j}\psi^{b}_{j}\psi^{c}_{j}\psi^{d}_{j}$, i.e the ordered set $\{\mathcal{O}_{j}^{\alpha}\}_{\alpha=1}^{abcd}$ coincides with $I_{j}$ in Eq. (38).

Eq. (39) defines the previously announced mapping to a system of $N$ qudits, as one can interpret

[TABLE]

where $\mathcal{K}_{N}=h_{1}\otimes\ldots\otimes h_{N}$ and $h_{j}\simeq\mathbb{C}^{8}$ is the space generated by $\{\ket{1}\,,\ket{ab}\,,\ldots\,,\ket{abcd}\}$. In this picture, the differential equation (34) is equivalent to a quench problem in $\mathcal{K}_{N}$: the system is prepared in the initial product state

[TABLE]

and left to evolve according to the differential equation

[TABLE]

Here, $H$ [not be confused with $H_{\rm SYK}$ in (1)] is an operator acting on $\mathcal{K}_{N}$ which plays the role of the Hamiltonian driving the imaginary-time dynamics. The precise form $H$ in terms of local operators can be derived by computing the action on the basis operators (37) of the left multiplication by $L$ in (35); however, even without doing this explicitly, it is easy to show that $H$ is invariant under any permutation of the sites in $\mathcal{K}_{N}$. This comes from the fact that the operator $L$ in (35) is left unchanged under the exchange of the pairs $\psi^{\alpha}_{i}\psi^{\beta}_{i}$ and $\psi^{\alpha}_{j}\psi^{\beta}_{j}$ for any choice of $i$ and $j$. Since the initial state (41) also enjoys the same symmetry, we can conclude that the dynamics of $\ket{\overline{\mathcal{U}}(t)}\rangle$ takes place in the subspace $\mathcal{S}_{N}\subset\mathcal{K}_{N}$ which is invariant under arbitrary permutations of the sites. This is of course a great simplification for our problem. The permutational symmetry of the Hamiltonian $H$ is “emergent” in the sense that it manifests itself only after taking averages over the Brownian disorder, while the Hamiltonian $H_{\rm SYK}$ in Eq. (1) does not display this symmetry for individual random realizations.

In order to study the dynamics in this subspace, we introduce the basis vectors $\ket{\mathcal{O}_{\vec{n}}}\rangle$ for $\mathcal{S}_{N}$

[TABLE]

where we used the same notations as in Eqs. (39), (40). Here $\pi$ is the unitary operator associated with a generic element in the symmetric group $S_{N}$, whose action permutes different sites in $\mathcal{K}_{N}$. Note that, since the sum runs over all the permutations, not all the elements in the sum are linearly independent.

The basis vectors (43) of the permutation symmetric space $\mathcal{S}_{N}$ are labeled by sets of $8$ integers $\{n_{j}\}$, satisfying $\sum_{k}n_{k}=N$, where each integer $n_{k}$ [with $k=1\,,ab\,,\ldots\,,abcd]$ “counts” the number of qudits in the level associated with $k$. In fact, it is possible to employ a more convenient representation, by viewing the state (43) as an $8$-mode Fock state generated by bosonic creation operators acting on a vacuum $\ket{\Omega}$. In particular, we have the identification

[TABLE]

Here, each operator $a_{k}^{\dagger}$ creates a collective mode corresponding to the level associated with $k$. In this language, the initial state (41) is written as $\ket{\overline{\mathcal{U}}(0)}\rangle=\frac{1}{\sqrt{N!}}\left(a_{1}^{\dagger}\right)^{N}\ket{\Omega}$.

This representation is particularly convenient, due to the fact that the Hamiltonian $H$ in Eq. (42) can be written in terms of the same bosonic operators appearing in Eq. (44):

[TABLE]

where $X_{r}$ is a bilinear operator of bosons. The explicit form of $X_{r}$ is derived in Sec. V.1, cf. Eqs. (81)-(86). A formal solution to the problem of computing $\overline{\mathcal{U}}(t)$ is then obtained as

[TABLE]

From its explicit form, one can see that the Hamiltonian $H$ commutes with the operator $\sum_{j=1}^{8}a^{\dagger}_{j}a_{j}$, which “counts” the total number of bosonic modes; accordingly, the evolved state (46) always belongs to the finite-dimensional Hilbert space generated by the basis vectors (44). However, the dimension of the latter is $D=\binom{N+7}{7}$, which grows as $N^{7}$, strongly limiting any numerical computation based on a brute force implementation of Eq. (46). Luckily, it is possible to show that the Hamiltonian $H$ has additional symmetries, which are unveiled by means of an appropriate Bogoliubov transformation

[TABLE]

further reducing the dimension of the effective Hilbert space explored by the dynamics. The matrix element $C_{n,m}$ are reported in Sec. V.1 [cf. Eq. (87)]. Then, from the form of the Hamiltonian $H$ in terms of the modes $b_{j}$ [cf. Eqs. (88)–(93)], one obtains that the number operators

[TABLE]

are conserved, namely they commute with $H$. Of course, the initial state can also be expressed in terms of the modes introduced in Eq. (47). Using the explicit form of $C_{n,m}$, we obtain

[TABLE]

and find that the total conserved number $n_{1,2}+n_{3,4}+n_{5,6}+n_{7,8}$ is $N$.

As we will see in the next section, these formulas allow us to work with effective Hilbert spaces whose dimensions grow either linearly or quadratically with $N$, and hence to provide numerically exact results for very large system sizes.

III.3 The OTOC and the tripartite information

We now discuss the last step of our method, namely the computation of the matrix elements of the form (22). Let us consider the most general OTOC

[TABLE]

where we introduced

[TABLE]

and where all indices are different, i.e. the operators have only $p$ Majorana fermions in common. Considering Eq. (22), we can expand $\overline{\mathcal{U}}(t)=\overline{U\otimes U^{\ast}\otimes U\otimes U^{\ast}}$ into the basis of operators $\mathcal{O}_{\vec{n}}$ corresponding to the vector (43) in ${\rm End}(\mathcal{H}_{N}^{\otimes 4})$. We obtain

[TABLE]

where the sum runs over all the sets $\vec{n}=\{n_{j}\}$ with $j=1,ab,\ldots,abcd$ and $\sum_{j}n_{j}=N$, while $c_{\vec{n}}(t)$ are the coefficients of $\overline{\mathcal{U}}(t)$ in the basis of the operators $\mathcal{O}_{\vec{n}}$. One can now interpret the sum (54) as the scalar product between an appropriate state $\ket{W_{(p,n,m)}}\rangle\in{\rm End}(\mathcal{H}_{N}^{\otimes 4})$ and $\ket{\overline{\mathcal{U}}(t)}\rangle$; namely we can write

[TABLE]

The whole problem of extracting the numerical value of the OTOC from the knowledge of $\overline{\mathcal{U}}(t)$ then boils down to writing down explicitly $\ket{W_{(p,n,m)}}\rangle$. After this is done, one can straightforwardly compute the overlap (55).

The derivation of the explicit form of $\ket{W_{(p,n,m)}}\rangle$ is however rather technical, and for this reason we postpone it to Sec V.2. The final result, instead, is extremely simple, and reads

[TABLE]

where $\ket{\Omega}$ and $b_{j}$ were introduced in Eqs. (44) and (47) respectively.

Surprisingly, one can also express the Rényi-$2$ entropies entering in the definition of the tripartite information (14) in the same form. More precisely, choosing the same conventions as Fig. 2 for the bipartitions of input and output of the evolution operator, one can write

[TABLE]

Here $\bar{\ell}$ is the length of $B$ and $D$ (chosen to be of the same size), while we will use $\ell$ for the length of the regions $A$ and $C$.

Once again, we refer the reader to Sec. V, where this is explicitly shown, while in the following we report the final result of this analysis, which gives

[TABLE]

and

[TABLE]

Similar formulas could be in principle derived also for more general choices of the bipartitions of input and output. This, however, would introduce additional technical difficulties, so we don’t derive them here.

It is now very important to comment on the form of the formulas presented above, as it allows us to reduce the computational cost required to obtain the physical quantities of interest. Let us first consider the case of the OTOC (51), which is conveniently rewritten as

[TABLE]

Namely, in order to compute $\mathcal{F}_{(p,n,m)}(t)$, one can evolve $|W_{(p,n,m)}\rangle\rangle$ and then take the overlap with the state $\ket{\overline{\mathcal{U}}(0)}$. This is important, as is best appreciated by looking at the simplest instance $p=0,n=m=1$. In this case, Eq. (56) implies that the state $|W_{(0,1,1)}\rangle\rangle$ belongs to the sector of the Hilbert space labeled by the quantum numbers $n_{1,2}=N-1,n_{3,4}=1,n_{5,6}=n_{7,8}=0$, where $n_{i,i+1}$ were introduced in Eq. (48)–(49). Since $n_{j,j+1}$ are conserved by the Hamiltonian $H$, the dynamics takes place in this sector of the Hilbert space, whose dimension can be easily seen to be $D=N$. Accordingly, one can conveniently represent the restricted Hamiltonian in a basis consisting of $N$ elements, and compute $e^{Ht}|W_{(0,1,1)}\rangle\rangle$ in this basis, which allows us to go to system size one million.

Similar considerations hold for the generic OTOC $\ket{W_{(p,n,m)}}\rangle$ (which belongs to the sector $n_{1,2}=N-p-m$, $n_{3,4}=p+m$, $n_{5,6}=n_{7,8}=0$) and for the Rényi-$2$ entropies corresponding to (59), (60). In the latter cases, expanding

[TABLE]

one is left with a sum of terms, each of which requires a simulation within Hilbert spaces up to dimensions $N\bar{\ell}\sim N^{2}$. Putting all together, we see that the computation of the quantities of interest requires us to simulate the dynamics in a Hilbert space whose dimension grows either linearly (for the OTOCs) or quadratically (for the tripartite information) with $N$.

IV The physical results

In this section we present the main physical results of our work. We begin with the analysis of the OTOCs, and continue with the Rényi-$2$ tripartite information introduced in Eq. (14).

IV.1 The OTOCs: numerical results

We start by presenting our numerical results for the simplest OTOC

[TABLE]

Due to the infinite range of the interactions and the disorder averages, $\mathcal{F}_{x,y}(t)$ does not depend on the precise choice of $x$ and $y$, but only on whether $x=y$ or $x\neq y$. Both cases are displayed in Fig. 5, where we report data for increasing values of the system size $N$. We see that $\mathcal{F}_{x,x}(t)$ and $\mathcal{F}_{x,y}(t)$ (with $x\neq y$) behave qualitatively differently at short times: the former displays an initial exponential decay, while the latter appears to remain approximately constant. In fact, based on the formulas of the previous section, one can make these statements more precise and show

[TABLE]

where the convergence is point-wise in $t$. This is proven in Sec. (V.3). In both OTOCs $\mathcal{F}_{x,x}(t)$ and $\mathcal{F}_{x,y}(t)$, we see the emergence of a characteristic time $t^{\ast}(N)$, increasing with $N$, which is required before they begin to decay towards zero at large times. One naturally interprets $t^{\ast}(N)$ as a scrambling time, which is also consistent with our subsequent analysis of the tripartite information. Finally, in Fig. 5$(c)$ we plot together the OTOCs for $x=y$ and $x\neq y$, for different systems sizes. We see that after an initial time window, the two OTOCs become indistinguishable, meaning that the information on the initial operators chosen has been completely washed out by the chaotic dynamics.

In order to quantitatively characterize the dependence of the scrambling time $t^{\ast}(N)$ on the system size, we test the short-time behavior of $\mathcal{F}_{x,y}(t)$ against the analytical ansatz MaSS16

[TABLE]

where $c_{x,y}$ is a constant (independent of $N$). In particular, we compute

[TABLE]

and compare the numerical data against a linear behavior. The results are shown in Fig. 6$(a)$. We clearly see that as the system size is increased, the curves for $\ell_{N}(t)$ approach the linear fit $y(t)=\frac{4}{3}t+\ln(2)$, within an initial time interval that is also increasing with $N$. In turn, this means that the ansatz (67) is valid, with the free parameters fixed as

[TABLE]

From this result, we can identify the scrambling time with $t^{\ast}(N)=3\ln(N)/4$.

The initial behavior in Eq. (67), together Fig. 5$(b)$, suggests that a data collapse should take place if we consider the shifted functions

[TABLE]

where we assumed that the parameters (69) are exact. This is plotted in Fig. 6$(b)$, where we see a remarkable data collapse at all times. In particular, the data appear to be perfectly collapsed already for $N\gtrsim 800$.

Next, we have tested how robust the above predictions are, against different choices of the local observables. We have considered in particular

[TABLE]

We have verified that at short times the ansatz (67) is always valid, and that a data collapse always takes place using the shift in Eq. (70). Furthermore the exponent is universal, namely it is independent of the observables chosen (while the prefactor is not). However, the OTOCs corresponding to distinct choices of local operators are quantitatively different, also after the scrambling time $t^{\ast}(N)$, as it can be appreciated from Fig. 6$(c)$. This can be interpreted by saying that, at finite times, the system retains some information on the initial observable chosen.

In order to investigate this point further, we plot in Fig. (7)$(a)$ different OTOCs, corresponding to distinct choices of local observables, which are labeled according to the convention of in Eq. (51). The curves correspond to initial operators that all have the same length, namely that are product of the same number of fermions. In this case, we see that all the OTOCs converge to the same function (up to small corrections) after the scrambling time. Comparing with the results displayed in Fig. 6$(c)$, we can conclude the following: after the scrambling time, information regarding the specific initial observables is lost, whereas OTOCs corresponding to operators with different initial length can still be distinguished.

Finally, we have investigated the large-time exponential decay of the OTOCs. The data in Fig. 5 suggest to consider an ansatz of the form

[TABLE]

where $\tau_{N}$ should be asymptotically independent of $N$. In Fig. 7$(b)$, we plot $\ln(-\mathcal{F}_{x,y}(t))$ for large values of $t$, and we see that the data are indeed consistent with an exponential decay of $\mathcal{F}_{x,y}(t)$. To be quantitative, we have performed a fit of $\ln(-\mathcal{F}_{x,y}(t))$ using $r_{N}(t)=a-t/\tau_{N}-b/t$. For the values of time $t$ available, we have found that the fitted $\tau_{N}$ has a weak dependence on $N$, with $\tau_{N}\simeq 1.53\pm 0.04$ for $N\simeq 10^{5}$. The fitted value appears to be independent of the choice of the local observables, up to the inaccuracy of the extrapolation method.

IV.2 The Rényi-$2$ tripartite information

We finally present our results for the Rényi-$2$ tripartite information introduced in Eq. (14). As we discussed in Sec. III.3, for this quantity the effective dynamics to be computed takes place in a Hilbert space whose dimension grows quadratically with $N$, so that we are restricted to smaller system sizes than in the case of OTOCs. Furthermore, for large subsystems the value of the entropy becomes very large, so that we also have to deal with issues of numerical precision. Overall, for the computationally worst case of bipartitions of equal size, we are able to provide data up to $N\simeq 400$. More details on the numerical implementations are reported in Appendix C.

In Fig. 8 we present data for the time evolution of the Rényi entropies of the subsystems $A\cup C$ and $A\cup D$, where we used the same partitions as Fig. 2. The plots correspond to fixed subsystem size and increasing $N$. Based on the formulas of Sec. III, in this limit we are able to compute (cf. Sec. V.3)

[TABLE]

We see from Fig. 8$(a)$ that the numerical results are in perfect agreement with this prediction. For finite $N$, the entropy $S^{(2)}_{AC}(t)$ displays an initial linear increase, eventually reaching a saturation value, as expected from the traditional picture of quantum quenches. The behavior of the Rényi entropy $S^{(2)}_{AD}(t)$ is instead not monotonic, as displayed in Fig. 8$(b)$. Indeed, one has $S^{(2)}_{AD}(0)=(N/2)\ln 2$, which is the maximum entropy possible, so that at small times $S^{(2)}_{AD}(t)$ has to decrease. Its dynamics is then non-trivial during the initial scrambling time $t^{\ast}(N)$, after which it begins an exponential decay towards its large-time stationary value.

Figs. 9 and 10 show the same quantities for all the possible values of the subsystems $\bar{\ell}$, at different times and system sizes. First, we notice that the entropies and the tripartite information are symmetric under exchange $\bar{\ell}\leftrightarrow\ell=N-\bar{\ell}$, as they should. Furthermore, we see that for different values of $\bar{\ell}$ we have the same qualitative behavior, where at large times an asymptotic value is always reached. In fact, it is possible to compute the latter exactly, as it is known that unitary evolutions driven by Brownian Hamiltonians converge in the infinite-time limit to unitary $k$-designs, for arbitrary positive integers $k$ BaBK17 ; OBKB17 . As a consequence, the asymptotic properties can be computed using Haar averages. The latter, which were already computed in Ref. HQRY16 , are reported as dashed lines in Fig. 9 and 10, towards which convergence is apparent. We note that, while their infinite-time limit could be expected, the entropies undergo non-trivial dynamics at short and intermediate times. This is best appreciated by looking at the entropy $S^{(2)}_{AD}(\bar{\ell})$ in Fig. 10. We see that up to the scrambling time $t^{\ast}(N)$ it appears to be decreasing (precisely, its average over $\bar{\ell}$), while at later times it increases again. This results in the non-trivial dynamics of the tripartite information, which can become positive at short times [cf. Fig. 8$(c)$].

V Deriving the key formulas

In this last section, we finally address the most technical aspects of our calculations, including several details of the method outlined in Sec. III. We start by presenting the explicit form of the Hamiltonian driving the dynamics in the four-replica space in Sec. V.1. Next, we derive the key formulas (56), (59) and (60) in Sec. V.2. Finally, in Sec. V.3 we report the proof of Eqs. (65) and (74) for the large-$N$ limit of the OTOC $\mathcal{F}_{x,y}(t)$, and of the Rényi entropy $S^{(2)}_{AC}(\bar{\ell})$.

V.1 The Hamiltonian

In this section we show how to derive the explicit form (45) of the Hamiltonian driving the imaginary-time evolution in Eq. (42), from Eq. (35). We start with the identity

[TABLE]

with $X=\sum_{i=1}^{N}x_{i}$, for commuting operators $x_{i}$ satisfying $x_{i}^{2}=-1$. This can be easily derived as follows (see e.g Ref. SaSS18 ). First, define

[TABLE]

Then, using $x_{j}^{2}=-1$, it is straightforward to show

[TABLE]

which immediately yields the desired identity. Eq. (75) allows us to write the Hamiltonian in terms of global sums of pairs of single-site Majorana operators.

Next, suppose that for a single-site operator $x_{i}$ we have

[TABLE]

where $\mathcal{O}^{\alpha}_{j}$ have been defined after Eq. (39). Then one can make the following identification

[TABLE]

namely the action of $X$ on the permutation symmetric basis (43) is the same as the r.h.s. of Eq. (79), as can be checked directly. From this, the final form of the Hamiltonian in terms of bosonic modes $a_{j}$ is readily obtained, and reads

[TABLE]

where $(-1)^{\gamma_{r}}$ is given in (36), while the operators $X_{r}$ are defined as

[TABLE]

Inspection of Eq. (80) reveals that the Hamiltonian displays several conservation laws. It is natural to look for a Bogoliubov transformation of the modes which makes some of the symmetries apparent. In addition, one would also like this transformation to simplify the convoluted $a$-mode expression $\ket{\mathcal{W}_{(p,n,m)}}\rangle$ for the OTOCs (111). Motivated by this, we look for a transformation where the first boson $b_{1}$ is associated with the macroscopically occupied mode in Eq.  (111), and choose the other modes $b_{j}$ to satisfy canonical commutation relations. While this can be done in different ways, it turns out that a particularly convenient transformation is the one defined by Eq. (47), where $C_{n,m}$ is the element in the line $n$ and in the column $m$ of the matrix

[TABLE]

Indeed, after this Bogoliubov transformation the form of the Hamiltonian immediately reveals the presence of additional symmetries which can be directly exploited for our computations.

It is straightforward to rewrite the operators (81)–(86), and hence the Hamiltonian (80), in terms of the new modes $b_{j}$. In particular, we have

[TABLE]

V.2 Extracting OTOCs and Entropies

We now wish to show how to derive an explicit expression for the vectors $|W_{(p,n,m)}\rangle\rangle$, $|W_{S^{(2)}_{AC}(\bar{\ell})}\rangle\rangle$ and $|W_{S^{(2)}_{AD}(\bar{\ell})}\rangle\rangle$ in Eqs. (56), (59) and (60), respectively. In order to simplify this task, we start by proving the following lemma. Let

[TABLE]

where $z_{x}\in\{1,ab,ac,ad,bc,bd,cd,abcd\}$ is the operator at site $x$ for the permutation $\pi$, c.f. (43) and $\alpha_{x}(z_{x})$ constants. Then

[TABLE]

The equivalence between Eqs. (95) and (96) is best established by directly expanding the product in Eq. (96), and regrouping the different terms.

Next, we introduce some notations to handle our subsequent calculations in a compact way. In particular, let us rewrite the basis operator $\mathcal{O}_{\vec{n}}$ in Eq. (43) as

[TABLE]

Here we introduced the notations $\Psi^{ab}_{ab}=\prod_{p\in I_{ab}}\psi_{p}^{a}\psi_{p}^{b}$, $\Psi^{ac}_{ac}=\prod_{p\in I_{ac}}\psi_{p}^{a}\psi_{p}^{c}$, $\ldots$, $\Psi^{abcd}_{abcd}=\prod_{p\in I_{abcd}}\psi_{p}^{a}\psi_{p}^{b}\psi_{p}^{c}\psi_{p}^{d}$, where $I_{ab}$, $I_{ac}$, $\ldots$ $I_{abcd}$ are ordered, pairwise disjoint, subsets of $\{1,2,\ldots N\}$, such that $|I_{ab}|=n_{ab}$, $\ldots$ $|I_{abcd}|=n_{abcd}$. In this notation, upper indexes in $\Psi^{\alpha}_{\beta}$ indicate the type of single-site operators, while lower indices specify which subset of $\{1,2,\ldots,N\}$ the product of such operators runs over. Consistent with this convention, we also introduce

[TABLE]

where $I_{a}$ is the ordered set defined by

[TABLE]

and whose elements are ordered as they appear in Eq. (99) (namely, the first $n_{ab}$ elements of $I_{a}$ are those of $I_{ab}$ with the same order, followed by those of $I_{ac}$, and so on). Analogously, one can define $\Psi^{b}_{b}$, $\Psi^{c}_{c}$ and $\Psi^{d}_{d}$, where

[TABLE]

with the elements of $I^{a}$, $I^{b}$ and $I^{c}$ ordered as they appear in Eqs. (100), (101) and (102). Finally, let us consider two disjoint subsets $A\cup B=\{1\ldots N\}$. Then, we define

[TABLE]

and analogously for the other cases. Using these notations, we can rewrite

[TABLE]

where $\mathcal{N}=(N!n_{1}!\cdots n_{abcd}!)^{-1/2}$ is the normalization. In order to write down the first line, we used that even string of different fermions commute, while sorting the Majorana operators in the second line resulted in the phases $(-1)^{\gamma_{A}}$, $(-1)^{\gamma_{B}}$. We will not write $\gamma_{A}$, $\gamma_{B}$ explicitly, as they will cancel at the end of the calculations. Conversely, one can easily compute the phase $(-1)^{\delta}$ appearing in the last line of Eq. (105):

[TABLE]

Here we have used that $n_{a}=n_{aB}+n_{aA}$ and $n_{aB}+n_{bB}+n_{cB}+n_{dB}$ are even, which can be seen by writing explicitly $n_{aB}=n_{abB}+n_{acB}+n_{adB}+n_{abcdB}$ etc.

Eq. (105) is the starting point to derive the explicit form of the vectors $|W_{(p,n,m)}\rangle\rangle$, $|W_{S^{(2)}_{AC}(\bar{\ell})}\rangle\rangle$ and $|W_{S^{(2)}_{AD}(\bar{\ell})}\rangle\rangle$ for OTOCs and Rényi entropies respectively. These are treated in the following, in dedicated subsections.

V.2.1 The OTOCs

We wish to calculate the OTOC (51) of the initial operators (52), (53). Starting from (22)-(24), we can insert (105) for the correlated time evolution operator, where simply $A=\{1\ldots N\},B=\{\}$ such that there is only one $(-1)^{\gamma}$ and no $(-1)^{\delta}$ and we omit the labels $A,B$. Since this expression involves products of even numbers of Majorana fermions acting on the same “replica” space, we may switch to the operators $\chi_{j}^{\alpha}$ in Eqs. (25), (26), and perform backwards the steps to derive Eq. (22) from Eq. (18). This gives

[TABLE]

where $\Phi^{(p,n)}$, $\Phi^{(p,m)}$ are defined in Eqs. (52), (53). Here, we simply wrote $\Psi_{a}$, $\Psi_{b}$, $\Psi_{c}$ and $\Psi_{d}$ without superscript, as we only have a single copy of the fermionic space. More explicitly, we have, for instance

[TABLE]

where $I_{a}$ is defined in (99).

Next we move the operator pairs $\Phi^{(p,n)}$ and $\Phi^{(p,m)}$ together such that they cancel. Of course, this generates phases through the anti-commutation relations of the fermions; we obtain

[TABLE]

where $n_{a,c}\left(\{i_{\alpha}\}\right)$ is the number of indeces in $\{i_{\alpha}\}_{\alpha=1}^{p}$ which also belong to $I_{a}\cup I_{c}$ [as defined in Eqs. (99), (101)]. Analogously, $n_{a,b}\left(\{j_{\alpha}\}\right)$ and $n_{b,c}\left(\{k_{\alpha}\}\right)$ are, respectively, the numbers of indexes in $\{j_{\alpha}\}_{\alpha=1}^{n}$ and in $\{k_{\alpha}\}_{\alpha=1}^{m}$ which also belong to $I_{a}\cup I_{b}$ and $I_{b}\cup I_{c}$. Finally, noticing

[TABLE]

we see that the factor $(-1)^{\gamma}$ in (109) is exactly canceled. We are left with an equation of the form (95), setting $\alpha$’s appropriately. Thus we may apply the lemma proved before [cf. Eq. (96)], which directly gives

[TABLE]

The result (56) is finally obtained after expressing the operators $a_{j}^{\dagger}$ in terms of the $b$-modes, introduced in Eq. 47.

V.2.2 The Rényi-$2$ entanglement entropy $S^{(2)}_{AC}({\bar{\ell}})$

Next, we turn to the task of deriving the vector $|W_{S^{(2)}_{AC}(\bar{\ell})}\rangle\rangle$ introduced in Eq. (57), corresponding to the exponential of the second Rényi entropy $S^{(2)}_{AC}({\bar{\ell}})$.

As we have explained in detail in Sec. II.2, in order to compute $S_{AC}^{(2)}(\bar{\ell})$, we need to consider the evolved state $U^{a}\ket{I^{ab}}$, where $\ket{I^{ab}}$ was introduced in Eq. (12), while the time evolution operator $U^{a}$ acts only on the “replica” space $a$ (cf. Fig. 3). For our fermionic system, the reduced density matrix of the union of the disjoint sets $A$ and $C$ can then be written as (see e.g. FaCa10 )

[TABLE]

Here we denoted by $\{F^{a}_{A}\}$, and $\{F^{b}_{C}\}$ a complete basis of operators in $A$, $C$ respectively; namely $F^{a}_{A}$ and $F^{b}_{C}$ take value in all the possible strings of Majorana operators supported in $A$ and $C$. Here, as before, we followed the convention that upper indexes indicate the type of single-site operators, while lower indices specify which subset of $\{1,2,\ldots,N\}$ the product of such operators runs over. Through simple manipulations, we have

[TABLE]

where $U^{\alpha}_{\pm}(t)$ are defined in Eq. (31). From this equation we clearly see that, in complete analogy with the case of the OTOCs, we can write also $\operatorname{tr}\left[\overline{\rho_{AC}^{2}}\right]$ in the form $\bra{L}\overline{\mathcal{U}}(t)\ket{R}$. As anticipated, this allows us to apply a procedure similar to the one employed for the OTOCs, and derive the vector $|W_{S^{(2)}_{AC}(\bar{\ell})}\rangle\rangle$. In particular, we can use the notations introduced in Sec. V.2, and exploit directly Eq. (105). This yields straightforwardly

[TABLE]

Next, we compute

[TABLE]

In (V.2.2), all Majorana operators are in the same system, so we leave away the doubled system label. To extract the traces, note that those are the only operators acting on the region $B$ of the system. Thus the Majorana fermions on those sites $B$ already have to cancel in pairs, or the expectation value $\braket{I}{\cdot}{I}$ will be zero. Similarly, the only value of $F_{C}$ with non-zero contribution has

[TABLE]

which can be inserted into the second expectation value, canceling the $\pm 1$ and giving

[TABLE]

Here, in order to go from the first to the second line, we made use of the identity (151) in Appendix B.2. The last line of (117) is non-zero only for

[TABLE]

With this we see that the traces evaluate as $(-1)^{\gamma_{B}+\gamma_{A}}$. Also, it shows that $n_{aB}\equiv n_{bB}\equiv n_{bA}\equiv n_{cA}\equiv n_{cB}\equiv n_{dB}\ ({\rm mod}\ 2)$, such that $(-1)^{\delta}=+1$. Putting all together, we get

[TABLE]

where the Kronecker delta $\delta(\pi)$ enforces the constraints (118). This expression can also be cast into the form (95) by setting some $\alpha$’s to zero. Then Eq. (96) gives us

[TABLE]

Transformation to $b$-modes (47) finally yields the result anticipated in Eq. (59).

An analogous treatment can be carried out for the case of the entropy $S^{(2)}_{AD}$. Since the technical steps are very similar, we report them in Appendix D.

V.3 Some large-$N$ limits

In this section, we finally show how one can compute the limit $N\to\infty$, while keeping time $t$ fixed, for the OTOCs and the Rényi-$2$ entropies, and derive in particular Eqs. (65) and (74).

We start with the case of OTOCs, and consider Eq. (61). As a first simplification, we only need to keep modes $b_{1}$ through $b_{4}$ in the Hamiltonian and initial state $\langle\bra{\overline{\mathcal{U}}(0)}$ as the others are not present in $\ket{W_{(p,n,m)}}\rangle$. Next, we switch to ladder operators with an unusual normalization, specifically

[TABLE]

This now allows us to take the leading order in $N$ for each term of the exponential $e^{Ht}$, using that $\tilde{b}_{1}\sim N$ as $p,n,m\ll N$. We obtain

[TABLE]

with

[TABLE]

As $\langle\bra{\overline{\mathcal{U}}(0)}(\tilde{b}_{2}^{{\dagger}2}\tilde{b}_{1}^{2}/N^{2}-1)=0$ at the highest order in $N$, terms with $H_{A}$ do not contribute at the leading order. For the OTOC $\mathcal{F}_{x,y}(t)$, also $H_{B}$ and $H_{C}$ cannot occur, because $\ket{W_{x,y}}\rangle$ (56) does not contain any $b_{3}$-modes. The asymptotic result is then the constant $\mathcal{F}_{x,y}(t)\to\langle\braket{\overline{\mathcal{U}}(0)}{1}{W_{x,y}}\rangle=-1$, as reported in (66). In contrast, for the OTOC $\mathcal{F}_{x,x}(t)$, the state $\ket{W_{x,x}}\rangle$ does contain one $b_{3}$-mode such that $H_{B}$ can appear at most once. The remaining Hamiltonian is still simple enough to finally derive the exponential decay (65). We stress that we can only expect these limits to be point-wise in $t$ due to the exchange of limits in (123); in fact, convergence is clearly not uniform, as can be seen from the exact numerical results.

The case of the entropy $S_{AC}^{(2)}(t)$ is treated along similar lines. We first perform a further mode transformation

[TABLE]

such that

[TABLE]

We may now follow the same procedure as for the OTOCs. In fact, the Hamiltonian has the exact same form in terms of the modes $b_{j}$ and $c_{j}$. Taking $\bar{\ell}\ll N$, Eq. (124) is therefore valid, after substituting the modes $\tilde{b}_{j}$ with $\tilde{c}_{j}$. Now, the initial state

[TABLE]

annihilates both $H_{A}$ and $H_{B}$ ensuing in a very simple (quadratic) asymptotic Hamiltonian $H_{C}$. From this, Eq. (74) follows straightforwardly.

VI Conclusions

In this work, we have developed an approach to analyze the chaotic dynamics in the Brownian SYK model, a system of $N$ Majorana fermions coupled together via random, time-dependent interactions. We have shown that the OTOCs and the tripartite information of the unitary evolution can be studied as a quench problem (at imaginary times) in a system of $N$ qudits, which can be conveniently investigated in terms of bosonic modes, due to an emergent permutational symmetry. Exploiting the latter, we were able to produce numerically exact results up to $N=10^{6}$, and to study several features of the chaotic dynamics at finite size.

We have analyzed in detail the dependence of the OTOCs on the observables chosen, highlighting the pieces of information on the initial operators which are not washed out by the chaotic dynamics. In particular, after the scrambling time $t^{\ast}(N)\sim\ln N$, the OTOCs of distinct operators converge to the same curve if they have the same length, namely if they are written as products of the same number of Majorana fermions, whereas the curves of different OTOCs can be distinguished after the scrambling time $t^{\ast}(N)$ if the length is different. Furthermore, we have verified that the exponent of the initial exponential growth of the OTOCs is universal and performed a data collapse for increasing system sizes. Regarding the tripartite information, we have shown that its evolution is non-trivial during the initial scrambling time, while at large times it always decays exponentially to the corresponding Haar-scrambled value; this result is consistent with the rigorous recent findings of Refs. BaBK17 ; OBKB17

The approach developed in this paper can be generalized to other models where the Hamiltonian displays all-to-all random interactions, with time-dependent Brownian disorder. Indeed, one can straightforwardly follow the steps outlined in Sec. III, and study the dynamics of OTOCs and tripartite information as a quench problem in a qudit system with site permutational symmetry. In turn, this implies that the effective imaginary-time dynamics takes place in a Hilbert space whose dimension grows as a polynomial in $N$. Of course, one would need to investigate for each case whether a further reduction of the effective dimension takes place, as for the Brownian SYK model studied in this paper.

It is possible that the final formulas obtained with our method (which have been used in this work mainly for efficient numerical computations) could be simplified further and evaluated to exact analytic expressions in the large-$N$ limit. In fact, by means of a different approach, an exact result for a suitable average of OTOCs was found in Ref. ZhCh19 for the Brownian dynamics generated by a disordered Hamiltonian in a qudit system. It would be interesting to see whether ideas related to the work ZhCh19 could be used here, to obtain analytic expressions for the OTOCs of arbitrary observables and for the tripartite information, in the large-$N$ limit.

Finally, the approach presented in this paper could also be applied to compute quantities involving higher moments of the evolution operator $U(t)$, such as Rényi entropies of higher order, or the Rényi-$2$ operator entanglement entropy of local observables PrPi07 ; Duba17 . In these cases, however, the application of our method would be inevitably more complicated. More importantly, it is not granted that a reduction of the Hilbert-space dimension could be achieved by means of a transformation analogous to (47). In any case, it would be certainly interesting to investigate these points further.

Acknowledgements.

LP thanks Bruno Bertini and Tomaz Prosen for discussions related to this work. XLQ thanks the helpful discussions with David Huse and Alex Streicher. LP acknowledges support from the Alexander von Humboldt foundation. JIC and NS acknowledge support by the EU Horizon 2020 program through the ERC Advanced Grant QENOCOBA (No. 742102, JIC) and the ERC Starting Grant WASCOSYS (No. 636201, NS), and from the DFG (German Research Foundation) under Germany’s Excellence Strategy - EXC-2111 - 390814868. XLQ is supported by the National Science Foundation under grant No. 1720504, and in part by the Department of Energy under grant No. DE-SC0019380.

Appendix A Non-interacting case: $q=2$

In this section, we study the Brownian SYK model (1) for $q=2$. We choose the constant $\sigma_{J}$ in (2) such that the disorder’s correlations are given by

[TABLE]

Each disorder realization is governed by a free Hamiltonian, therefore we do not expect any scrambling of operators or decay of OTOCs.

The method developed in this article can be applied to arbitrary $q$ and we may study the non-interacting case within its framework. The states

[TABLE]

representing the OTOC (56), Rényi-2 entanglement entropies (57) and (58), and initial time evolution operator (50) are independent of $q$ as long as $q$ is even. However, the effective Hamiltonian reflects the change of $q$ and is simpler. Along the same lines as for $q=4$ (see section III.2), we can compute

[TABLE]

The corresponding representation of the effective Hamiltonian after operator-state mapping in bosonic modes is

[TABLE]

where the six $X_{r}$ operators are the same as in the corresponding expression (45) for $q=4$.

For the two simple OTOCs $\mathcal{F}_{x,y}(t)$ and $\mathcal{F}_{x,x}(t)$ the dynamics only explores the two-level subspace spanned by $\ket{N-1,0,1,0}$ and $\ket{N-2,1,0,1}$. Therefore we can compute these OTOCs analytically and obtain

[TABLE]

The curves are plotted in Fig. 11. As expected, the OTOCs do not decay to zero at long times, as the non-interacting model is not chaotic. Since the tripartite information can be written as an average of OTOCs HQRY16 , it too will lack the characteristics of scrambling.

We now make a comment on the so-called “length” of the operator $\psi_{j}(t)$, see e.g. RoSS18 . At any time $t$, we can always write

[TABLE]

and define the average length $L(t)$ as

[TABLE]

It can be shown that the average length is related to an appropriate average over OTOCs, namely RoSS18

[TABLE]

Using this relation, if follows from our results (138)–(139) that the length is constant 1. This is expected because the Gaussian dynamics preserves the length of products of Majorana operators.

Next, we can calculate the entanglement entropies numerically, just like in the interacting case. We present the results in Fig. 12. In the limit $N\to\infty,{\bar{\ell}},t$ fixed, we can derive

[TABLE]

along the same lines as in section V.3. While the entropy $S_{AC}^{(2)}({\bar{\ell}})$ saturates to its maximal Haar value at large $N$ and $t$, the behavior of $S_{AD}^{(2)}({\bar{\ell}})$ is qualitatively different from the interacting case (Fig. 8). This leads to the tripartite information being positive at all times and system sizes, indicating the absence of scrambling.

As for the interacting case, we can also study the entanglement entropies’ dependence on the subsystem size, see Fig. 13. Comparing against Figs. 9 and 10, we see that in the free case, the entanglement entropies do not reach the maximal Haar scrambled values at finite ratios $\bar{\ell}/N$.

Appendix B Relation between OTOCs and Rényi-$2$ entropies

In this appendix we review the relation between OTOCs and Rényi-$2$ entropies in the case of unitary evolution operators defined on qubit systems, and generalize the latter for fermionic (Majorana) systems. We focus on the configuration displayed in Fig. 2, taking the regions $A$ and $C$ of the same size and position (and analogously for $B$ and $D$).

B.1 The case of Pauli matrices

We start by giving a derivation of the aforementioned relation for a system of qubits, along the lines of the one in Ref. HQRY16 . First, we can write the reduced density matrix $\rho_{AC}$ as

[TABLE]

where the sum is over the complete bases $\{O_{j}^{A}\}$ and $\{O_{k}^{C}\}$ of strings of Pauli operators in $A$ and $C$, while $a$ and $c$ are equal to the number of sites in $A$ and $C$. The state $\ket{I}$ is the maximally entangled state connecting $A\cup B$ and $C\cup D$, satisfying

[TABLE]

while $U_{AB}$ is the evolution operator acting non-trivially only on the system $A\cup B$. Then, using the orthogonality of the Pauli operators, and after simple simplifications, we have

[TABLE]

Consider now the sum

[TABLE]

Using the identity HQRY16

[TABLE]

(here $A_{j}$ are a complete basis of operators for the Hilbert space associated with $A$, while $|A|$ is its dimension), one immediately obtains that the r.h.s. of Eq. (146) is equal to (147). Therefore

[TABLE]

In the last step, we summed over $k$, used once again the identity (148), and finally renamed the indexes $\tilde{k}=l$. Putting all together, we find

[TABLE]

This is exactly the same result as in HQRY16 . An analogous derivation holds for the case of $S^{(2)}_{AD}$. This equation encodes a close connection between the tripartite information and the OTOCs, and allows one to establish that chaos, as measure by small values of all OTOCs, implies scrambling HQRY16 . In the next section we show that a similar relation, with the addition of proper signs, holds in the case of fermionic systems.

B.2 The case of Majorana operators

For Majorana operators one needs a different treatment. Indeed, the identity (148) is no longer valid, but should be modified as follows. Let $\mathcal{O}$ be an operator with a well defined parity of Majorana operators, i.e. $\mathcal{O}$ is the sum of strings of operators that are either all even or all odd. Then, by expanding in the operator basis of Majorana operators, one can prove

[TABLE]

Here $\ell_{j}$ is the length of the operator $A_{j}$. For example, if $A_{j}=\psi_{1}\psi_{4}$ then $\ell_{j}=2$. Analogously, $\ell_{\mathcal{O}}$ is the length of one of the terms in $\mathcal{O}$. Since all these terms have the same parity, it does not matter which one we choose. We can then proceed as in the previous sections, now paying attention to the order of the operators involved in the calculations. First, we have

[TABLE]

where now $a$ and $c$ are half the number of sites in $A$ and $C$. Proceeding as before, we have

[TABLE]

where $n=N/2$. Consider now the sum

[TABLE]

Noticing now that the evolution operator can always be written as sum of even strings of Majorana operators, one can directly apply the identity (151) to prove that the r.h.s. of (153) is equal to (154). On the other hand, using

[TABLE]

we can sum over $k$ by employing once again the identity (151), and finally rename $\tilde{k}=r$. Putting all together, we obtain

[TABLE]

We see that additional signs appear in the sum over the OTOCs with respect to the case of Pauli matrices. However, the Rényi-$2$ entropy still encodes global information about the sum of the OTOCs over extended regions of the system.

Appendix C Details on the numerical implementation

In this section, we explain a few details for the numerical computation. For the OTOCs, we implement (61) in terms of the modes $b_{1},b_{2},b_{3},b_{4}$, as explained in the main text. For the entropies however, we use different modes. While we could likewise implement (57) and (58) in $b$-modes, this introduces large numerical error as $N$ increases, due to cancellations of large numbers. Instead we have found that for the entropy $S_{AC}^{(2)}$, a numerical calculation in terms of modes $c_{1},c_{2},c_{3},c_{4}$ (128) is most stable. For the entropy $S_{AD}^{(2)}$, we found that using the modes $b_{1},b_{2},c_{3},c_{4}$ gives the most stable results. For our numerical calculation we have therefore transformed initial state $\ket{\overline{\mathcal{U}}(0)}\rangle$ (50), effective Hamiltonian $H$ (80) and $\ket{W_{S_{AC,BD}^{(2)}({\bar{\ell}})}}\rangle$ (59)-(60) into these modes.

Appendix D The Rényi-$2$ entanglement entropy $S^{(2)}_{AD}({\bar{\ell}})$

In this appendix, we turn to the task of deriving the vector $|W_{S^{(2)}_{AD}(\bar{\ell})}\rangle\rangle$ introduced in Eq. (58), corresponding to the exponential of the second Rényi entropy $S^{(2)}_{AD}({\bar{\ell}})$. The discussion goes along the same lines of the one presented in Sec. V.2.2 for the entropy $S^{(2)}_{AC}({\bar{\ell}})$. Writing out the partial trace as done for the other entropy, we get

[TABLE]

where we denoted by $\{F^{a}_{A}\}$ and $\{F^{b}_{D}\}$ a complete basis for the operators in $A$, $D$ respectively; namely $F^{a}_{A}$ and $F^{b}_{D}$ take value in all the possible strings of Majorana operators supported in $A$ and $D$. We can continue along the same lines as for the other entropy, giving

[TABLE]

Here the term $(*)$ can be evaluated as for the other entropy up until (V.2.2). Then, however, $\Psi_{bB}^{a{\dagger}}$ and $\Psi_{aB}^{a}$ are not the only parts in the first expression acting on this subspace, now $F_{D}^{a}$ also does. So, continuing from $(*)$ and dropping the doubled system label as all operators are in the same system, we have

[TABLE]

The left side is only non-zero for

[TABLE]

such that we can evaluate the sum $\sum_{F_{A},F_{D}}$, inserting these in the right side. We get

[TABLE]

The traces now give $(-1)^{\gamma_{A}}$ and $(-1)^{\gamma_{B}}$, respectively. Inserted back into (V.2.2), these cancel, and the $\delta$ partially cancels the other phases,

[TABLE]

Again, this has the form (95) with suitable choices of $\alpha$’s. Then, according to Eq. (96) we obtain

[TABLE]

Finally, the transformation to $b$-modes in (47) results in Eq. (60).