Probing quantum chaos in multipartite systems

Understanding the emergence of quantum chaos in multipartite systems is challenging in the presence of interactions. We show that the contribution of the subsystems to the global behavior can be revealed by probing the full counting statistics of the local, total, and interaction energies. As in the spectral form factor, signatures of quantum chaos in the time domain dictate a dip-ramp-plateau structure in the characteristic function, i.e., the Fourier transform of the eigenvalue distribution. With this approach, we explore the fate of chaos in interacting subsystems that are locally maximally chaotic. Global quantum chaos can be suppressed at strong coupling, as illustrated with coupled copies of random-matrix Hamiltonians and of the Sachdev-Ye-Kitaev model. Our method is amenable to experimental implementation using single-qubit interferometry.

A prominent signature of quantum chaos is the repulsion among energy levels. For instance, the spacing between nearest-neighbor levels follows the Wigner-Dyson distribution in quantum chaotic systems, while it is described by Poisson statistics in the presence of conserved quantities (e.g., in integrable systems) [38]. The SFF is proportional to |Z(β + it)/Z(β)| 2 , where Z(·) is the partition function and β = 1/k B T . This quantity probes the level statistics of both close and far-separated energy eigenvalues, providing a tool to detect the ergodic nature of the system [1]. For a generic chaotic system, the SFF exhibits a dip-ramp-plateau structure [see e.g., Fig. 1]. Its short-time decay forms a slope. The physical origin of the subsequent ramp is the long-range repulsion between energy levels [11]. The transition from the slope to the ramp forms the dip. The final plateau originates * zhenyuxu@suda.edu.cn † adolfo.delcampo@uni.lu from the finite Hilbert space dimension and approaches a constant value Z(2β)/Z(β) 2 in the absence of degeneracies in the energy spectrum. The SFF has been widely employed in the study of the discrete energy spectrum of quantum chaotic systems [11,[39][40][41][42][43][44][45]. Quantum chaotic systems composed of multipartite subsystems subject to generic interactions typically have a complicated energy spectrum [45][46][47][48]. We shall focus on a global subsystem composed of strongly-chaotic subsystems, interacting with each other. In this setting, any subsystem can be seen as an open quantum system embedded in an environment, composed of the remaining subsystems. The subsystem dynamics is thus governed by dissipative quantum chaos [1], which is currently under exhaustive study [29,30,33,34,37,[49][50][51][52][53]. We shall depart from the standard practice of assuming an effective open quantum dynamics, as ubiquitously done in the literature. Instead, we will account for the exact unitary dynamics of the global composite system, with no approximations (e.g., without invoking the Markovian description or an effective master equation). The above diagnostics can be employed to detect global quantum chaos in multipartite systems [48]. However, apart from proposals like the fidelity-based SFF [34,37] and the related partial SFF [45], they are not suited to directly detect how chaotic behavior stems from the subsystems and their interactions. In this work, we provide an experimentally realizable approach to this end by considering the measurement of an energy observable X, which can be the Hamiltonian of a subsystem or the interaction energy. As measurement outcomes are stochastic, we propose to study the full counting statistics, characterized by the eigenvalue distribution of X at thermal equilibrium. Its Fourier transform, the characteristic function, reveals chaotic behavior through the dip-ramp-plateau structure. Its analysis shows that strong interactions among the different subsystems can suppress the global chaotic behavior of the multipartite system, even when the subsystems are maximally chaotic, as revealed by the study of global and local observables. This scheme not only provides a convenient theoretical tool to diagnose quantum chaos in complex multipartite quantum systems, but it can be experimentally realized by using single-qubit interferometry with an ancillary qubit. Our paper is organized as follows. In Sec. II, we introduce the characteristic function of an energy observable, which is then used as a tool to detect the chaotic features stemming from the subsystems and their interactions in Sec. III. Then, we employ this method to analyze the chaotic behavior in multipartite systems sampled from the Gaussian orthogonal ensemble (GOE) in Random Matrix Theory (RMT) [2,54] in Sec. IV and the coupled Sachdev-Ye-Kitaev (cSYK) models [55][56][57][58][59][60][61][62][63][64] in Sec. V. Finally, we summarize in Sec. VI with concluding remarks and a brief discussion of potential applications.

II. THE CHARACTERISTIC FUNCTION OF THE ENERGY OBSERVABLE X
Let H = l H l be the Hilbert space of a multipartite system and X be an energy observable of a local subsystem in the subspace k H k ⊆ H. We focus on the Hamiltonian of the subsystems (local energy) and the interaction energy, as choices of the observable X. The probability distribution of the observable X with eigenvalues {x} averaged over an initial thermal equilibrium state ρ th = e −βH /Z(β) is [65] P (x) = tr [ρ th δ(X − x)] . (1) The eigenvalue probability distribution P (x) encodes the full counting statistics of the observable X, that is, the probability to find the system in an eigenstate with eigenvalue x when prepared in the state ρ th . In terms of the integral representation of the Dirac delta function, the probability distribution can be expressed as the Fourier transform of the characteristic (moment generating) function that captures the statistical properties of the spectrum of the observable X. While in the following we refer to t as a time variable, it is to be understood as the Fourier conjugate to x. Experimentally, the characteristic function Eq.

III. A PROBE FOR QUANTUM CHAOS IN MULTIPARTITE SYSTEMS
In what follows, we consider the absolute square value of the generating function in Eq. (3), i.e., as a tool for probing quantum chaos in multipartite systems, identifying contributions from the subsystems and their interactions. The choice X = H j represents the local energy that can be employed to diagnose the chaotic behavior contributed by the jth subsystem in an N -partite system. By contrast, the observable X = H j ⊗ H k ⊗ H l can be used to detect signature of quantum chaos attributed to the interactions among subsystems j, k, and l.
with no obvious dip-ramp-plateau structure plays an important role in attenuating the chaoticity.
can be interpreted as the fidelity between the initial coherent Gibbs state (or a thermofield double state) and the state resulting from its evolution [34,44]. In addition, when X is a small perturbation of the Hamiltonian H, i.e., H = H 0 + X, and commutes with H (or H 0 ), Eq.
which captures the overlap between two identical initial states (|ψ 0 ) evolving under slightly different Hamiltonians H and H 0 [17]. Note that according to the Baker-Campbell-Hausdorff formula, Eq. (5) and the LE differ in the general case, when [X, H] = 0. We emphasize that the observable X in Eq. (5) is not required to represent a small perturbation or to commute with H. When it describes the local energy of a subsystem k H k ⊆ H, Eq. (5) provides the possibility to directly detect the chaotic behavior contributed by the subsystems or interactions to the global multipartite system. It can be used to either diagnose quantum chaos of a one-partite system (as done by the SFF) or a structured multipartite system. To support this observation, we illustrate its use in the following examples involving coupled random-matrix Hamiltonians and the coupled SYK model.

IV. PROBING THE CHAOTICITY IN COUPLED RANDOM-MATRIX HAMILTONIANS
Consider a N −partite system with a general Hamiltonian of the form For the sake of illustration, let us sample the Hamiltonians from the Gaussian orthogonal ensemble (GOE) [2, 54], which is a paradigmatic random matrix ensemble for physical applications involving systems with timereversal symmetry and exhibiting quantum chaos. GOE is the ensemble of real symmetric matrices, whose elements are chosen at random from a Gaussian distribution. The joint probability density of H j ∈ GOE(d j ) (d j denotes the dimension of the Hilbert space H j ) is proportional to exp(− 1 2σ 2 trH 2 j ), where σ is the standard deviation of the random matrix elements of H j .
The first example we consider is N = 1 and X = H. In this scenario, we focus on H ∈ GOE(d) and H ∈ GOE(d 1 )⊗GOE(d 2 ) (d 1 d 2 = d) as an example. Equation (5) averaged over the full GOEs [i.e., GOE(d)] yields where · represents the ensemble average and . = denotes the annealing approximation [11]. In Eq. (7), the GOE averaged partition function is given by (see Appendix A) where I n (·) is the modified Bessel function of first kind and order n, and the coefficient reads As shown by the red dotted curve (or the dark blue curve by numerical simulations) in Fig. 1, Eq. (7) exhibits a typical feature of quantum chaos, namely a dipramp-plateau structure. The early decay from unit value comes to an end, forming a dip (also known as correlation hole) with the onset of a ramp. The latter extends until it saturates at a plateau value at the characteristic plateau time t plateau = 2 √ d/σ [see Eq. (9)]. The existence of the ramp, a period of linear growth of G(t) GOE , is a consequence of the repulsion between long-range energy levels [11]. This long-range repulsion causes the energy levels to be anticorrelated. The plateau stems from the discrete energy spectrum, whose height is Z(2β) GOE / Z(β) 2 GOE . Similar chaotic features have been studied in the Gaussian Unitary Ensemble (GUE) [34,44,80], Gaussian ensembles under infinite temperature [81], and Sachdev-Ye-Kitaev (SYK) models [11,82,83].
For systems with less chaoticity than the full GOE, the span of the ramp will shrink or even disappear (see light blue curves in Fig. 1).
Without loss of generality, we then consider a bipartite system with H (0,1) 1,2 independently sampled from GOEs. As the total system is composed of two partitions each described by a random-matrix Hamiltonian, it is not surprising that in the absence of (or weak) interactions the full system exhibits visible dip-ramp-plateau structure when choosing X = H [see the dark blue curve in Fig.  2 (a)]. However, when the coupling strength 1 between subsystems 1 and 2 is enhanced, the dip-ramp-plateau structure gradually washes out [light blue curves in Fig.  2 (a)].
To account for this phenomenon, we look at the characteristic function for different choices of the observable X = H (0) 1 in Fig. 2 (b) and X = H Fig.  2 (c) and show how these choices identify the contributions to quantum chaos by the first subsystem and the interactions between subsystems 1 and 2. Obviously, H (1) plays an important role in diminishing the chaotic behavior, since the characteristic function reflects no obvious ramp structure. Indeed, from the perspective of the nearest-neighbor level distribution, the Kronecker product of random matrices will tend to break the Wigner-Dyson distribution under certain conditions [84]. When the coupling is enhanced, the interaction term H (1) dominates, and the chaoticity of the whole system is gradually suppressed.
Similar phenomena exist in more structured systems, as shown in Fig. 3  (2) 1,2,3 and omit the superscript therein. Both bipartite interactions (e.g., H 1 ⊗ H 2 ), and the tripartite interaction H 1 ⊗ H 2 ⊗ H 3 play an import role in decreasing the chaoticity of the composite global system, while the local chaotic nature of each subsystem can be detected by choosing a local observable (e.g. X = H 1 ) even at strong coupling.
It is worthwhile to note that the presence of interactions could also induce chaos if the interaction tend to mix the subsystems, just as shown in the coupled kicked rotors [46].

V. COUPLED SACHDEV-YE-KITAEV MODEL
The second example we consider is the coupled Sachdev-Ye-Kitaev (cSYK) model [55][56][57][58][59][60][61][62][63][64]. A system composed of 2N Majorana fermions is divided into two separated sides and each subsystem is described by the SYK model [85,86]. We consider the left and right SYK Hamiltonians H L,R with a bilinear coupling H b . The total Hamiltonian of the cSYK model reads where µ controls the strength of the the bilinear coupling and χ k denote Majorana fermion operators satisfying the anticommutation relations {χ k , χ l } = δ kl . J klmn and K kk are random coupling constants independently sampled from Gaussian distributions with zero expectation values and J 2 klmn = 3!J 2 N 3 , K 2 kk = K 2 N 2 . A similar model has been used to study the holographic duality of an eternal traversable wormhole: By preparing the SYK model in a thermofield double state and turning on the coupling between the two sides, the wormhole is traversable in the context of gravity [55]. Figure 4 shows the numerical result for the disorderaveraged G(t) , in agreement with the random matrix theory. From Fig. 4(b), we see that the bilinear coupling plays an important role in controlling the chaotic behavior of the whole system, as the generating function has no obvious dip-ramp-plateau structure. The chaotic character of the system is robust when the bilinear coupling is weak (µ 0.1). When enhancing the coupling, the dipramp-plateau structure of the entire system is gradually washed out [ Fig. 4(a)]. This implies that the entire system becomes less chaotic when the bilinear coupling is strong, even when the chaoticity of the subsystems remain consistent. As in the GOE example, the chaotic nature of each subsystem can be detected by choosing a local observable (e.g. X = H L ), as shown in Fig. 4(b).

VI. DISCUSSION AND CONCLUSION
We have introduced a protocol to directly detect quantum chaos in interacting multipartite systems by measuring the statistical distribution of an energy observable X at thermal equilibrium. Specifically, we make use of the absolute square value G(t) of the generating function of the eigenvalue distribution associated with the observable X. When the observable equals the total Hamiltonian, G(t) reduces to the SFF. For local observables, chaotic features give rise to a dip-ramp-plateau structure in G(t), which is similar to that in SFF. G(t) directly detects the contributions to quantum chaos in a composite system from different subsystems by choosing the observable for k-partite interactions. We have shown that the coupling of chaotic systems can give rise to the suppression of quantum chaos in the composite system, as the interaction strength among the subsystems is increased. From the perspective of decoherence, sampling the eigenvalue statistics of a subsystem is like sampling the local energy. Quantum chaos is generally expected to be suppressed as a result of decoherence [34,87]; see however [37]. In addition, even at strong coupling, the chaotic character of the subsystems can be unveiled by choosing X as a local observable, as demonstrated by considering the multipartite GOEs and the coupled SYK models.
Our scheme can be implemented in quantum devices, such as NMR systems [69,72,74] and trapped ions [75] by introducing an auxiliary qubit coupled to the systems, as both the random spin and SYK models are realizable in the laboratory [88][89][90][91][92][93][94]. Our approach for diagnosing the chaos in multipartite quantum systems may thus find broad applications in interdisciplinary studies in quantum information, quantum matter, and AdS/CFT duality, especially in analyzing quantum chaos in structured quantum many-body systems.
In Appendix A, we intend to briefly introduce the calculation of the generating function G(t) averaged over ensembles. The averaged G(t) in terms of annealing approximation is given by [34,44] This annealing average is in agreement with the quenched average |Z (β + it)| 2 /Z(β) 2 in high temperature region [11,34]. Then the denominator and numerator of Eq. (A1) can be written as and where ρ(E) is the spectral density and ρ(E, E ) is the two-point probability density function.

Gaussian orthogonal ensemble statistics
For GOE, the spectral density and two-point probability density function are given by [34] ρ and (A5) respectively, and K d (x, y) is the kernel [2]. Note that GOE and GUE share the same form of N -point probability density function. The only difference is that the kernel K d (x, y) in Eqs. (A4) and (A5) for GOE is a quaternion. Note that σ is selected as σ = 1/ √ 2 in Ref.
According to Dyson's theorem [2], and (A8) The above method can be straightforwardly extended to higher point correlation functions. Similar work has been done in Ref. [81] for evaluating spectral form factors under infinite temperature.

G(t) averaged by GOEs
With Eq. (A7), the partition function [Eq. (A2)] averaged over the GOEs is approximated by where I n (·) is the modified Bessel function of first kind and order n (Eq. (8) in the main text). According to Eqs. (A8), the imaginary time partition function [Eq. (A3)] averaged over the GOEs reads In the main text, we consider fixed number of realizations and temperature for illustration. In this section, we append Fig. 5 and Fig. 6 to illustrate G(t) with the dependence of the number of realizations and temperature.   Fig. 2(a). Data for G(t) is averaged over 1, 10, and 100 realizations of GOE.