Random Matrix Theory of the Isospectral twirling

In this paper, we present a systematic construction of probes into the dynamics of isospectral ensembles of Hamiltonians by the notion of Isospectral twirling, expanding the scopes and methods of ref.[1]. The relevant ensembles of Hamiltonians are those defined by salient spectral probability distributions. The Gaussian Unitary Ensembles (GUE) describes a class of quantum chaotic Hamiltonians, while spectra corresponding to the Poisson and Gaussian Diagonal Ensemble (GDE) describe non chaotic, integrable dynamics. We compute the Isospectral twirling of several classes of important quantities in the analysis of quantum many-body systems: Frame potentials, Loschmidt Echos, OTOCs, Entanglement, Tripartite mutual information, coherence, distance to equilibrium states, work in quantum batteries and extension to CP-maps. Moreover, we perform averages in these ensembles by random matrix theory and show how these quantities clearly separate chaotic quantum dynamics from non chaotic ones.


Introduction
Few concepts are more fascinating than the one expressed by the word Chaos (χάος). Originally meaning 'the Abyss', but soon opposed to the notion of an ordered universe [2], the concept of chaos lives up to its etymology by being a very challenging notion to define precisely. In classical mechanics, chaotic systems are those displaying high sensitivity to the initial conditions. A related notion is that of the butterfly effect, which relates the amplification of small local perturbations to the spreading of a large effect involving the whole system. In quantum unitary dynamics, different initial states cannot lead to exponentially separated trajectories, as an elementary result from the fact that the inner product between different states is conserved by unitary dynamics. To circumvent this elusiveness, quantum chaos has traditionally been associated to salient properties of those systems which were quantized from classical chaotic systems. Such systems were found to possess statistics of energy levels spacings corresponding to the predictions of Random Matrix Theory (RMT) [3][4][5][6][7][8][9][10][11][12][13][14]. Such a measure of chaos is impractical for a quantum many-body system, as it requires a perfect knowledge of the spectrum of the Hamiltonian, and is quite indirect. After all, we still would not know what quantum chaos means [15][16][17][18][19].
In recent years, more direct measures of chaos have been proposed. The one that is most closely related to the butterfly effect is the expectation value (in the thermal state) of commutators of spa-tially separated local operators. A quantum version of the butterfly effect means that dynamics spreads operators around the system in a way that measuring one would have non trivial effect on the other [20]; this behavior in turns corresponds to the decay of four-point out-of-time order correlation functions (OTOCs) . This behavior has been confirmed for several theories believed to be chaotic, for instance, theories holographically dual to gravity [45][46][47][48] or random unitary evolutions modelling the scrambling behavior of black holes [49][50][51][52][53], to non integrable quantum many-body systems [16,[54][55][56][57][58][59][60][61][62][63]. The scrambling of information in a quantum many-body systems consists in the delocalization of quantum information by a quantum channel over the entire system in the sense that no information about the input can be gained by any local measurement of the output. It has been a very insightful and remarkable result [64] the realization that the butterfly effect as measured by rapid decay of the OTOCs implies the scrambling of quantum information. When different properties proposed to define a notion come to a unification, this means that one is onto something.
Quantum chaos is also invoked as an explanation of thermalizing behavior in closed quantum many-body systems [17,[65][66][67][68][69][70][71][72][73][74][75][76][77][78][79][80][81][82][83][84] and the emergence of irreversibility [85] resulting from the establishment of a universal pattern of entanglement corresponding to that of a random state in the Hilbert space. The sensitivity of the dynamics can also be captured by the Loschmidt Echo (LE), and it has been shown that for bipartite systems LE and OTOCs coincide [86]. Entanglement dynamics is also invoked as an explanation for the onset of quantum chaos, from the aforementioned typicality of entanglement to complexity of the entanglement itself [87], and the fact that scrambling itself can be measured through quantum correlations between different parts of the system [40,64,88,89]. With a very insightful concept, the authors in [90] posit that the web of quantum chaos diagnostics should be unified and study further connections of these probes with complexity.
The need for a unifying framework is pressing, and if all the probes to quantum chaos can be unified in a single notion this would constitute a decisive step towards a theory of quantum chaos. An attempt towards unification has been made in [1], where it was shown that all the probes of quantum chaos can be cast in the form of the 2k−Isospectral twirling (IT). Consider a unitary evolution U = exp −iHt generated by the Hamiltonian H. All the moments of polynomials of order 2k that depend on such evolution can be constructed starting from U ⊗k,k ≡ U ⊗k ⊗ U †⊗k . The Isospectral twirling consists in averaging over all the possible U with a given spectrum, that is, the Haar average of a k−fold unitary channel. After the averaging [1], one obtains the probe as the expectation value in the thermal state of the IT with an appropriate permutation of the desired operator. After this average, this quantity only depends on the spectrum Sp(H) of the Hamiltonian generating the dynamics.
All the proposed probes to quantum chaos share a similar heuristics, that is, one has a desired characteristic that quantum chaos should have, say, the butterfly effect, and then defines a quantity that is able to characterize that effect. Moreover, one tries to show that physical systems thought to be chaotic do indeed possess that characteristic. In this very enlightening framework, though, one piece is missing. How does one know that these quantities and properties are not possessed also by systems that are far from being chaotic? The probability of return [91] was one of the first quantities conjectured to describe the high sensitivity of the dynamics. However, also integrable systems like diagonalizable quantum spin chains feature a rapid decay of such probability, see [92]. In the domain of quantum many-body physics, it has been shown that some OTOCs may indeed behave similarly also for such systems [41,93]. In order to complete the program of characterizing quantum chaos, one needs a tool that at the very least distinguishes clearly between chaotic and integrable systems. The Isospectral twirling does exactly this. Since it is a quantity that is completely characterized by the spectrum of the Hamiltonian generating the dynamics, one can evaluate it for Hamiltonians with spectra given by random matrix theory, e.g., the Gaussian Unitary Ensemble (GUE) and for integrable Hamiltonians with spectra belonging to the ensembles Poisson or the Gaussian Diagonal Ensemble (GDE) and compare. If the probes behave distinctly in the different ensembles, then one has truly characterized quantum chaos. In [94] the OTOCs are averaged over the GUE. Comparison between different spectra requires the use of random matrix theory for all the ensembles of Hamiltonians.
In this paper, we show that random matrix theory applied to the Isospectral twirling of a unitary quantum channel shows that all the proposed probes to quantum chaos clearly separate dynamics generated by GUE Hamiltonians, that are chaotic, from integrable ones, that certainly are not chaotic. The distinction can be seen in the salient values and corresponding timescales of certain spectral functions that are averaged over the corresponding ensembles of isospectral Hamiltonians. These values and time scales clearly separate the ensemble of GUE Hamiltonians that are chaotic from the Poisson and GDE ensembles in the way they depend on the dimension d of the Hilbert space. A complete table of the results is shown in Table 2. The paper is organized as follows: in Sec. 2 we introduce the steps involved in the calculation of the Isospectral twirling. In Sec. 3 we study the Isospectral twirling for the frame potential, OTOCs and Loschmidt echo, coherence and Wigner-Yanase-Dyson skew information, the mutual and tripartite mutual information, and quantum batteries. In Sec. 4 and the Appendices the details of the calculation of the Isospectral twirling and the random matrix theory is presented. Conclusions follow.

Definitions and strategy 2.1 General approach and notation
We now briefly describe the structure of this paper and provide a birdeye view to its scopes, techniques and results. The first goal of this paper is to show that many of the proposed probes to quantum chaos can be cast in the form of the expectation value of the Isospectral twirling times a suitable operator characterizing the probe. More specifically, we consider: • measures of randomness and adherence to Haar measure, like the frame potential (see Sec. 3.1).
• measures of the butterfly effect as OTOC and Loschmidt Echos(LE) (in Sec. 3.2), • measures of quantum correlations and direct measures of scrambling as the tripartite mutual information (in Sec. 3.3), • measure of coherence (in Sec. 3.4), • quantum thermodynamics measures of chaos (in Sec. 3.5).
We shall denote these measures by the notation P O (U ). Here P is a linear functional of an operator O whose expectation value may be of interest as a probe to quantum chaos, and U is the unitary evolution describing the dynamics of interest. We will show that, after averaging over all the possible dynamics with a given spectrum, all the probes to quantum chaos assume the general form The average · G is the average over all the unitary channels with a given spectrum. The quantitŷ R (2k) (U ) is the Isospectral twirling of U , that is, the average over its 2k−fold channel [1,95,96]. TheT ) is a rescaled permutation operator corresponding to the element π of the permutation group S 2k (in the cycle representation, see App. 4.2 for details). More precisely, the 2k−Isospectral twirling of a unitary channel is the average of (G † U G) ⊗k,k over G sampled uniformly from the unitary group. Let us define it formally. Be H C d a d−dimensional Hilbert space and let U ∈ U(H) be a unitary quantum channel. The 2k−Isospectral twirling of U is defined aŝ Above, dG represents the Haar measure over the unitary group U(d) and we recall that U ⊗k,k ≡ U ⊗k ⊗ U †⊗k . The Haar average will be computed by the Weingarten functions method, discussed in Sec. 4 resulting in a weighted sum over permutation operators: where the matrix with components Ω πσ = tr (T π T σ ) is the inverse of the Weingarten functions, see Sec. 4.1 for details. The connection between the unitary channel and the dynamics is given by the Hamiltonian of the system, that is, U = exp −iHt. Through this connection, the Isospectral twirling becomes a function of timeR (2k) (t). Effectively, the Haar average performed by the Isospectral twirling corresponds to averaging over all the possible eigenstates of H. The result of the twirl only depends on the spectrum of H.
From Eq. (3) we see that the dependence of the Isospectral twirling on the spectrum of the unitary channel U is stored in the functions known as spectral form factors [8] (tr T (2k) and therefore we haveR

Probes to Chaos
Let us now show the connection between the typical probes of quantum chaos and Eq. (1). The definitions of the quantities on the left hand side are given in Sec. 3. For every quantity we show the Isospectral twirling in the form of Eq. (1) and its asymptotic evaluation for large dimension d. The temporal dependence is contained in the coefficients c (2k) π (U ) which also contain all the information about the spectrum of the Hamiltonian generating the dynamics. The large d expression for some of the quantities above is too cumbersome to be reported in this section. We systematically compute and analyze all these quantities in Sec. 3.
The meaning of the above formulae is the following. We pick a random matrix within an ensemble of Hamiltonians defined by some spectral properties: for example, for a qubit system, we can consider all the Hamiltonians with the spectrum of a Pauli string. We ask: how will the dynamics generated by a Hamiltonian in this ensemble behave as described by quantities like mutual information, OTOCs, entanglement, coherence, etc. Knowing the spectrum of the Hamiltonian H, one can compute the coefficients c (2k) π (U ) that depend on the time t through U . The dynamics of the probes is on average specified solely in terms of the functions c (2k) π (t) in the Table 1. In the following, we will be interested by ensembles of Hamiltonians whose spectrum is fixed by some given probability distribution.

Probe
Isospectral twirling Large d limit

Spectral average of the Isospectral twirling: Random Matrix Theory
The second important goal of this paper is to show that the aforementioned probes do demonstrate peculiar characteristics of quantum chaos. To this end, we should show that a behavior that is possessed by an ensemble describing chaotic Hamiltonians is not possessed by an ensemble describing nonchaotic, e.g., integrable Hamiltonians. We can obtain this insight by looking at the average behavior of the probes P O (U ) G within that ensemble and showing that that behavior is typical. The Isospectral twirling is capable of being a tell-tale quantity for chaotic behavior because it only depends on the spectrum of the Hamiltonian, and we can define ensembles of chaotic or integrable Hamiltonians by spectral properties. For instance, for integrable systems typical energy distributions are known, the Gaussian Diagonal Ensemble (GDE) or the Poisson distribution (P). Similarly, for chaotic systems a typical spectrum is given by the Gaussian Unitary Ensemble (GUE) or the Gaussian Orthogonal Ensemble (GOE).
In the following we will denote by E the ensemble of Hamiltonians of interest and by R (2k) (t) E the average behavior of the Isospectral twirling in that ensemble. In turn, this will become an average over the ensemble E for the spectral functions c (2k) π (t). We will use random matrix theory for the ensembles E ≡ GUE, P, GDE. These calculations are presented in detail in Sec. 4.3. The spectral functions c (2k) π (t) depend on the permutation element. We will show they take the form where α, β, γ, δ ∈ N and depend on k and π and C d is a constant depending on the Hilbert space dimension d, see Table 3 and Table 4 for c (2) π (U ) and c (4) π (U ) respectively. Since we perform averages of these functions with respect to spectrum statistics of certain isospectral classes of Hamiltonians H, it is worth introducing some definitions. Let H = i E i |i i| be a spectral decomposition of the Hamiltonian, with the assumption E i+1 ≥ E i . Throughout the paper we assume H to be timeindependent and use a decomposition of the form U = j e −iE j t Π j , with Π 2 j = Π j orthogonal projectors. The spectral functions Eq. (4) one obtains, up to k = 2 are the following: In the above expression we have adopted a convenient notation for these functions which appear all around in the paper: c e (U ) = c 4 (t) and c (12) , it makes sense to introduce the tilded quantities for a = 2, 3, 4 defined abovẽ As the expressions Eq. (7) explicitly show, these coefficients are spectral functions of U ; let us introduce the spectral average we use. Let E be an isospectral family of Hamiltonians H with spectral distribution P ({E k }) ≡ P (E 1 , . . . , E d ), the spectral average of the Isospectral twirling over E is defined as: where d{E k } ≡ dE 1 · · · dE d ; since the Isospectral twirling depends linearly from the coefficients Eq. (7), see Eq. (3), its spectral average with respect to E is a linear combination of the spectral average of c π (U )s. One can immediately see that we can also express the coefficients in Eq. (7) in terms of (nearest) level spacings, s i ≡ E i+1 − E i ; for this reason we also consider isospectral families E with given distribution of nearest level spacings P ({s i }). Notice thatR (2) (U ) depends only on c 2 (t), whileR (4) (U ) depends on c 2 (t), c 3 (t) and c 4 (t). As such, the higher the order of the Isospectral twirling, the higher the amount of information that it is contained and its ability to distinguish different types of dynamics cfr. Table 2 and Fig. 2. A technical but important definition that we use later is the following. The spectrum of a Hamiltonian Sp (H) on H C d is said to be generic if for any l s.t. d ≥ l ≥ 1: unless E n i = E m j , ∀i, j = 1, . . . , l and for some permutation of the indices n i , m j . 1 Also, in some occasions also c2(2t) appears.

Eigenvalues distribution: GDE and GUE
We consider two types of eigenvalues distributions P ({E k }), corresponding to the following classes of matrices: GUE(d, 0, 1/d 1/2 ), random unitarily invariant d × d matrices whose entries have null mean and d −1/2 as standard deviation and GDE(d, 0, 1/2) random diagonal matrices whose d diagonal elements are stochastic independent variables with null mean and 1/2 as standard deviation. The eigenvalue distributions read [8][9][10]97]: the first distribution refers to a chaotic [9,98], while the second is an integrable spectrum [99]. The difference between these eigenvalue distributions is the presence of level-repulsion term, namely This term inhibits the presence of degeneracies and pops up also in the Wigner-Dyson distribution, as follows.

Nearest level spacings distributions
We have seen that the Isospectral twirling can be also regarded as a function of the nearest level spacings, P E ({s i }). Because of the dependence of the Isospectral twirling for generic Hamiltonian on the level spacings, this paper focuses on performing averages with respect Poisson (P), Wigner-Dyson from the Gaussian Orthogonal Ensemble (WD-GOE), and Wigner-Dyson from Gaussian Unitary Ensemble (WD-GUE). The nearest level distributions for these distributions are given by [3,8,14,[100][101][102][103][104]: Wigner-Dyson from GUE, which we use in the remainder of the paper, and whose distributions are shown Fig. 1. An important remark, which is cleared during the explicit calculations of random matrix theory (Sec. 4), is that for the distributions in Eq. (14) we are considering the nearest level spacings s i to be independent from each other; this assumption works well for the Poisson distribution.

Results and plots
In Sec. 4 we perform the computation of the spectral functionsc a (t). Here we discuss the results and present some plots for the temporal evolution of these functions. Since the goal of this paper is to infer from the Isospectral twirling the dynamical properties of integrable vs. non-integrable spectra, the properties ofc a (t) directly enter the Isospectral twirling. A striking feature ofc a (t) is that both the magnitude of the oscillations with respect to the mean, the location of the first minimum and the typical timescale with which the asymptotic state is reached depend on the dimension d of the Hilbert space. We present the salient features of the spectral functionc 2 (t) andc 4 (t) in Table 2. One can observe that typically the GDE distribution reaches the asymptotic limit in a √ log d time, while for  14). . the case of GUE and Poisson such timescale is a scaling law d γ , with exponents γ = 1 and γ = 1/2 respectively. While we are able to distinguish different ensembles during the dynamics, the behavior in t → ∞ is universal and does not depend on the ensemble (cfr. Table 2), see Proposition in Sec. 4.5. A detailed case by case analysis is presented below, while the spectral average of c 2 (t) and c 4 (t) are plotted in Fig. 2.

Typicality
In this section we show how typical is the behavior in the two ensembles of dynamical systems on which we perform the averages, namely, the typicality in the ensemble of isospectral Hamiltonians (that is, in the ensemble {G † HG, G ∈ U(d)}, and in the ensemble of Hamiltonians with a given spectrum, namely E.
The Isospectral twirling is a tool to study quantum mechanical quantities depending only on the spectrum of the unitary generating the dynamic, and thus the details, such as the basis of the Hamiltonian, are washed away in the calculation [68,75,86,105,106]. As we show in this section, however, the Isospectral twirling exhibits typicality, e.g. a typical dynamical curve is close to the one obtained from the Isospectral twirling. In this paper, given the operator P O , we analyze its Isospectral twirling P O G and then study the average behavior of such quantity over some well-known distribution of spectra and nearest level spacings distributions. Since we are taking two averages, there are two typicalities to talk about. With the use of Lévy lemma, we prove the typicality of the average with respect to the randomization of the eigenvectors, and from which the typicality of subsequent average with respect to the spectrum follows. Taking P O its Isospectral twirling can be written as P O G = tr (T (2k) OR (2k) (U )). Then: SciPost Physics Submissioñ   Table 2 (Left). Right: Log-Log plot of the spectral functionc 4 (t) E for different ensemble E = P , E = GDE and E = GUE for d = 2 12 . The starting value is 1, while the asymptotic value is (2d−1)d −3 . From these plots and from scalings in Table 2(Right) we conclude that c 4 (t) distinguishes chaotic from non-chaotic dynamics more efficiently than c 2 .
This bound holding for any t ≥ 0, e.g. taking t → ∞: Once again, in the large d limit we can set the probability to be arbitrary small. As a result of these calculations, the isospectral twirling presents typicality in both the averages over G and over E.

Isospectral twirling as a universal quantity for quantum chaos
We have discussed in the Introduction how challenging is a direct characterization of quantum chaos. In classical systems, chaos can emerge when phase space trajectories diverge exponentially due to small differences in their initial states [108]. Such exponential growth leads to fast mixing in phase space, and thus the information on the initial state is quickly lost in time. In classical dynamical systems, chaotic behavior is thus a possible pathway towards ergodic mixing [109] and thermalization. In quantum systems, the analog is not as obvious because of the linearity of the dynamics, which leads to unitary and seemingly information preserving evolution. Nonlinearity is an essential element of the theory of classical chaos, and thus there seems to be a certain tension between the two notions. However, the difference becomes thinner if the linear system is infinite dimensional. In quantum systems one can then intuitively see that in the case of many-body dynamics hints of rapid mixing between trajectories can occur. Such mixing in quantum systems is often referred to as scrambling [33,51,55,[110][111][112][113][114][115].
As we discussed earlier, many quantities have been proposed towards a diagnostic of quantum chaos [17,80,90,116,117], forming an intricate web [90] that necessitates a unifying framework.
A first probe to quantum chaos can be via the study of frame potentials, namely the adherence of an ensemble of unitary operators to the full unitary group. This quantity also serves as measure of randomness. We discuss the Isospectral twirling of the frame potential in Sec. 3.1.
A second possibility is to characterize quantum chaos through the growth of out-of-time-ordered correlation functions (OTOCs), thus analogous to the definition of diverging trajectories with the initial states [90]. OTOCs have various formulations, of which we discuss few in this paper [29,31,32,37,42,67,118].
Loschmidt Echos have also been proposed to probe the high sensitivity of a quantum dynamics. Recently [1,86], it has been understood that there is a direct connection between these quantities and the OTOCs. We study their Isospectral twirling in Sec. 3.2.
Entanglement and quantum correlations have been understood as important quantities for the description of complex quantum dynamics. Quantum chaotic dynamics is supposed to generate universal patterns of entanglement, on the one hand, leading to thermalization [49,119], while, on the other, to high entanglement complexity and emergence of irreversibility [85]. Correlations, both quantum and classical, in a quantum many-body system are described by the mutual information. Moreover, a more complex form of correlations defined by the tripartite mutual information [64,88] has been advocated as a measure of scrambling and thus of quantum chaos. These quantities are studied in Sec. 3.3 .
Recently, quantum coherence [120] has received a renewed interest as a resource theory [121,122] for quantum algorithms and quantum thermodynamics. In this paper, we show that also quantum coherence can probe chaotic dynamics. Quantum coherence has been studied with a similar goal in the recent paper [123]. The Isospectral twirling of quantum coherence and skew information is shown in Sec. 3.4.
One of the most important motivations for the study of quantum chaos is quantum thermodynamics. We have already mentioned the aspects of thermalization and entropy production. Thermalization in closed quantum systems has been associated to the typicality of the entanglement generated by typical Hamiltonians [17,65,67,68,70,71,74,75,78,79,119]. More generally [66], even in presence of conserved quantities, equilibration in a closed quantum system can ensue in the form of typical expectation values for local observables. We show how this kind of approach to equilibrium tells apart chaotic from integrable Hamiltonians in Sec. 3.5.1. Moreover, an important topic in quantum thermodynamics is the understanding of quantum engines [119] and their capability of storing energy or performing work, e.g., quantum batteries [96]. These systems are discussed under the lens of the Isospectral twirling in Sec. 3.5.
Finally, one can consider more general quantum evolutions, like those described by generic CPmaps and quantum channels [124]. We show in Sec. 3.6 how one can take the Isospectral twirling of such channels.
For all these quantities, we show how to cast them in the form of Eq. (1). Then, applying the results from random matrix theory of the previous Sec. 2.2, we show how these quantities discriminate quantum chaotic from non chaotic dynamics.

Frame potential
Random unitary operators have been used in the literature to approximate quantum chaotic dynamics, for instance, in the context of black holes [125] to study the "scrambling" of information. An important question is how a given set of unitaries samples well the unitary group, in its capability of representing the group averages up to the t−th moment of the distribution. This forms the notion of unitary t−design [23,118,[126][127][128][129][130][131][132][133][134][135][136][137][138]. A t−design is therefore 'random enough' up to the t−th moment. Quantum chaotic dynamics is supposed to be able to realize high t−designs unlike a less ergodic, non-chaotic dynamics.
Given E = {q j , U j } an ensemble of unitary operators U j ∈ U(H) weighted by q j , the k−fold channel of A ∈ B(H ⊗k ) is defined as: it is the weighted average of A under the adjoint action of the U j s. By definition, the ensemble of unitaries E is a k−design iff the k−fold channel Eq. (19) equals the Haar k−fold channel: where dU is the Haar measure over the unitary group. Recall that the Haar measure is the uniform measures over groups [139], see App. 4.1 for more details. In other words, a k−design reproduces k moments of the uniform distribution on U(H). A convenient characterization for k−designs is provided by the Frame potential F E of the ensemble E, introduced in [126] and defined: This quantity measures the 2−norm distance between the k−fold channel of E Eq. (19) and the Haar k−fold channel Eq. (20). Scott [126] proved that the Frame potential is minimized, namely F (k) This makes the frame potential a natural measure of the randomness for ensembles E.
Given one unitary U , we can construct the ensemble of all the unitaries with the same spectrum, namely {G † U G, G ∈ U(d)} and ask whether this ensemble realizes a unitary t−design. Since the unitaries are generated by a Hamiltonian, we can write G † U G = exp(−iG † HGt) and we can write the ensemble as We now compute the k−th frame potential of the ensemble E H to quantify its randomness, i.e. to see how well E H replicates the moments of the Haar distribution. The k-th frame potential of the ensemble E H reads: As it turns out, the frame potential of E H is the square Schatten 2−norm of the Isospectral twirling of Eq. (2) (see Proposition 1 in [1]): F (k) for completeness, a proof is reported in App. D.1. Recall that the computation of frame potential for the uniform Haar distribution returns F (k) Haar = k! and it provides a lower bound for Eq. (24), namely k! ≤ F (k) E H ; this consideration, together with Eq. (24), implies that a lower value for the 2−norm of the Isospectral twirling Eq. (2) can be associated to higher randomness in the dynamics for E H (t). One of the goal of this paper is to discriminate between chaotic and integrable Hamiltonians from spectral properties only; one of the key feature for the distinction comes from quantifying the randomness generated during the dynamics, i.e. as a function of t. A lower bound for the frame potential of E H which turns useful as a measure of the deviation from the (Unitary) Haar average is the following (see This bound holds for any t ≥ 0, e.g. if one takes the infinite time average: of the both sides of Eq. (25), it is possible to compute the lower bound for the asymptotic value of the frame potential, in the hypothesis of generic spectrum (see Sec. 2.2); from Proposition 3 in [1]: which is far from the Haar value k!. The non generic spectrum is defined in 2.2 and in practice means a particularly strong form of non resonance condition. The asymptotic value computed in Eq. 27 is the same for any Schwartzian spectral probability distribution, e.g. for GUE, GDE (cfr. E H is non trivial in its time evolution. For k = 1, we have : where the coefficientsc a (t) do depend on the particular choice of the spectrum of H; they are defined in Eq. (7) and computed in Sec. 4. From [1] we report the plot of the spectral average of the frame potential Eq. (28) for GUE, GDE and Poisson, see Fig. 3. While for Poisson and GDE the frame potential has a similar behavior, these are quite distinct from the non-trivial dynamics in the GUE ensemble. The frame potential for Poisson and GDE never crosses the asymptotic value 3 + O(d −1 ), and it always stays away from the Haar value 1. This is consistent with Proposition 4 in [1], which states that i.e. the randomness in time of the GDE ensemble is always far from being equal to the Haar value; a proof of this result can be found in [1], here we present an alternative proof, see App. D. In the case of GUE instead, the frame potential reaches the Haar value 1, and remains close to it during the temporal , only after reaching the asymptotic value 3. We conclude that, the chaotic spectra (GUE) present more randomness during the dynamics than the integrable ones (GDE, Poisson) which reach their maximum value of randomness asymptotically.

OTOCs and Loschmidt-Echos
In this section, we discuss two important classes of probes to chaos taking the form of correlation functions. In Sec. 3.2.1 we compute the Isospectral twirling of Loschmidt-Echos, while in Sec. 3.2.2 we compute the Isospectral twirling for the Correlators (OTOCs).

Loschmidt-Echo of the first and second kind
The Loschmidt Echo quantifies the non-trivial dynamics of a quantum system [91,92,140,141], it can be used to study revivals occurring when an imperfect time-reversal procedure is applied to a complex quantum system [142]. The Loschmidt-Echo of the second kind is defined as the squared Hilbert-Schmidt scalar product between the unitary time evolution e iHt and the perturbed backward evolution e −i(H+δH)t : It is interesting to use Hamiltonians H and δH such that the time evolution is non-trivial (for instance, chaotic), and thus the system is sensitive to the changes in initial conditions. In [86,90], it was found that, under suitable conditions, the OTOC (which we discuss next) and LE are quantitatively equivalent. Via the Isospectral twirling it was shown in [1] that the OTOC and the LE are very similar measures. We comment on this issue in the next subsection when we discuss OTOCs. For the Loschmidt Echo of the second kind, we now use the following connection to the Isospectral twirling (see Proposition 6 in [1]). Provided that Sp(H) = Sp(H + δH) the Isospectral twirling of the LE is given by: where A = A † ⊗ A. We thus see that the 4−point Isospectral twirling enters directly into the Loschmidt echo as well. If one sets A to be a Pauli operator on qubits one gets in the large d limit( see App. E): A plot of the Loschmidt-Echo of the second kind is in Fig. 4

(top).
A similar formula also applies to Loschmidt Echo of the first kind [92] L 1 (t), which is defined as the modulus square the overlap between |ψ and its time evolution |ψ(t) ≡ U |ψ , with U ≡ exp{−iHt}. Defining ψ = |ψ ψ|, it is straightforward to write it as: We can now obtain the Isospectral twirling; with the usual trick |tr (U ψ)| 2 = tr (U ⊗ U † ψ ⊗2 ), we obtain: We thus see that a 2-point Isospectral twirling enters into the Loschmidt Echo of the first kind. If we insert in this formula the average, obtained in Eq. (104), we obtain We immediately note that L 1 (t) depends onc 2 (t), while L(t) depends onc 4 (t) and therefore the Loschmidt-Echo of the second kind is more efficient in the distinction among chaotic and non-chaotic dynamics, cfr. Table 2. A plot of the quantity in Eq. ( 35) is shown in Fig. 5 (left); it is clear that we can distinguish chaotic from integrable time evolution.

Out-of-time-order correlators (OTOCs)
Among the many measures of information propagation in quantum many body physics, Out-of-Timeorder Correlators (OTOC) have recently found many new applications, ranging from black holes to disordered spin systems [33]. In particular, the "Information Scrambling" can be measured either via the OTOC [15], which we discuss here, or the tripartite mutual information (TMI) [88] which we discuss later (however, the decay of OTOC with time implies the decay of the TMI [64,88], and thus the two measures are related). In this section we show how the OTOCs are described by the Isospectral twirling. The definition of the 2k−OTOC is the following. Let us consider 2k local, non-overlapping operators A l , B l , l ∈ [1, k]. Then, the 4k−point OTOC is defined via the following expression: where operators evolve in the Heisenberg representation as A l (t) = e iHt A l e −iHt . Define A l := A † l ⊗ A l and similarly for B. The Isospectral twirling enters the 4k−point OTOC as (see Proposition whereT (4k) is a normalized permutation operator of S 4k defined asT (4k) π = ST 1 4k···2 S † , with S a permutation operator. For completeness we report the proof of Eq. (37) for k = 1 in App. E.1. For k = 1 we obtain the 4-point OTOC: If one sets A and B to be non-overlapping Pauli operators on qubits, one finds in the large d limit( see in App. E for details): As the expression above shows, the 4−point OTOCs distinguish chaotic from integrable behavior, through the timescales ofc 4 (t) (see Table 2(Right)); we observe that the integrable distributions, GDE and Poisson, reach the asymptotic value 1/d 2 in a time O( √ log d) and O(d 1/2 ) respectively; GUE instead reaches a negative value in a time O(d 1/2 ) and then the asymptotic value 1/d 2 in a time of O(1/d) ( see Fig. 4).
We now note that in this setting the LE and OTOC are equal to the spectral functionc 4 (t) up to a constant O(d −2 ). As a result, this implies that both the OTOC the LE are probing how scrambling the quantum system; in this sense, such an approach supports the findings of [94] on the quantitative agreement between OTOC and LE.
Another quantity of interest is the OTOC 2 (t). Let A be a local operators on B(H), the 2−OTOC is defined as: where Taking the Isospectral twirling is straightforward to see: where we have defined A = A † ⊗ A. A straightforward application of Eq. (104) shows that: The operator A is just a general Hermitian operator. In particular we can consider the case that A is a state. If we take A ≡ ψ to be a pure state, we find that the 2−OTOC and the Loschmidt-Echo of the first kind are the same quantity. The plot of the 2−OTOC for Poisson, GDE and GUE are presented in Fig. 5 (Right) for d = 2 16 . The convergence time can be observed to be different between the case of Poisson and GUE, and GDE.

Entanglement
The Isospectral twirling also enters the evolution of entanglement [106,[143][144][145][146][147][148][149][150][151][152][153][154][155] assuming that the evolution can be well described by a random Hamiltonian with a given spectrum . The setup we consider is the following. The unitary time evolution of a state ψ ∈ B(H A ⊗ H B ) is given by ψ → ψ t ≡ U ψU † . Given a bipartition of the total Hilbert space H = H A ⊗ H B of the Hilbert state, the entanglement of ψ t is computed by the 2−Rényi entropy S 2 = − log tr (ψ A (t) 2 ). As proven in [1] (see Proposition 7, and App. F.1), the 2−Renyi entropy has the following structure in terms of the Isospectral twirling: In the equation above, one gets, in the large d limit, ( see the details in App. F.1): The result we obtain using this technique is similar to the one obtained in [145]. It is thus clear that the time evolution of the entropy S 2 G is dictated byc 4 (t). An analysis of the behavior of the lower bound can be obtained via the asymptotic analysis of the coefficientc 4 (t) (see Table 2), for clarity a plot of the lower bound is shown in Fig. 6

Mutual information
The mutual information (MI) of two random variables is a measure of the dependence of one variable on the other, and viceversa, and it quantifies how much information can be obtained from measuring a certain subsystem on its complementary [89,124], e.g. a measure of correlations. Given a certain joint distribution p(x, y), the mutual information is defined as log p(x,y) p(x)p(y) , and quantifies the level of factorization of a certain joint probability [156] with respect to two random variables X and Y . Similarly, in quantum information theory, the quantum mutual information is a measure of correlation between subsystems and is defined via the von Neumann entropy S(A) = log ψ ψ , where ψ is the density matrix of a certain system A. It is the quantum mechanical analog of Shannon's mutual information, and is thus easily defined via the reduced density matrices of two complementary subsystems of the total Hilbert space. Specifically, let us thus consider a certain bipartition of the total Hilbert space H = H A ⊗ H B and a state ψ ∈ B(H). The mutual information between the two partitions of the system A and B is quantified by the mutual information, defined as: where S(A) is a measure of quantum entropy of the reduced density matrix ψ A ≡ tr B (ψ), the same for S(B). Instead of using von Neumann entropy, here we consider the 2-Rényi entropy S 2 (ψ) as a lower bound, for which we can perform the analysis analytically. We thus have: If one sets ψ = |ψ ψ| to be a pure state, the mutual information reduces to the entropy I(A : Therefore, taking a state ψ with a given purity tr (ψ 2 ) = p, let it evolve unitarily with a random unitary U G (t) with a given spectrum, from Eq. (43) we obtain the isospectral average: The bound on the mutual information is similar to the lower bound on the entanglement introduced in Sec. 3.3.1, as for the entanglement the behavior of the bound on the mutual information is dictated byc 4 (t). If we consider a bipartite system

Tripartite mutual information.
Among the quantum measures of scrambling and chaotic behavior for many-body systems, the tripartite mutual information (TMI) is an information-theoretic measure that has been recently proposed in the literature [64,88]. Consider a quantum channel that maps an input state in a bipartite system A good scrambling channel must be such that measurement on a local part of the output, say C, cannot detect any operation performed on a local part of the input, say A. For this reason, a good measure of scrambling is given by the tripartite mutual information [64]defined as where A, B and C, D are the fixed bipartitions of past and future time slices respectively. To compute this quantity, one first maps the quantum channel into a state by the Choi isomorphism. Because of the unitarity of the evolution, the reduced density matrices ρ AB and ρ CD of the Choi state are maximally mixed, and thus I(A; B) = I(C; D) = 0. Then, since the tripartite mutual information is defined as a conditional mutual information, one necessarily must have I 3 ≤ 0, which is guaranteed thanks to the sub-additivity property of the entropy, whichever one decides to use. The unitary is thus minimally scrambling when I 3 = 0. While one can work with many measures of entropy, here we focus on the 2-Rényi entropy measure I 3 (2) . The tripartite mutual information has several properties, see App. F.2 for more details. More specifically, ρ AC(AD) = tr BD(BC) (ρ U ), where ρ U is the Choi state of the unitary evolution U ≡ exp{−iHt} and H a random Hamiltonian with a given spectrum. Set A = C and B = D, then, by defining can be written as [1]: The second term of Eq. (49) is reminiscent of the entanglement of quantum evolutions, which has been introduced in [145]. It should come at this point as no surprise that also TMI can be written in the form of an Isospectral twirling. In fact (following from Proposition 8 in [1]), one has that the Isospectral twirling of I 3 (2) is upper bounded by Given the fact that, as mentioned above, the TMI is a negative-definite quantity, the decay of the r.h.s. of Eq. (50) implies scrambling, that is, its value gets close to the minimum. In the large d limit, Eq. 50 becomes (for details see App. F.2.3): It follows that we can, also in this case, evaluate I 3 (2) G E over the spectra GUE, GDE and Poisson. We In general, the time evolution of I 3 (2) G E depends on the spectral functionsc a (t) which we have defined above, and calculated explicitly in Sec. 4. The results are however shown in Fig. 7, where we can quantitatively see how the chaotic and integrable behaviors are clearly different. Clearly, the time-behavior of I 3 (2) (t) depend on the timescales ofc 4 (t) Table. 2(Right). The equilibrium value of Eq. (51) can be calculated explicitly in our approach. We find that, for large d: An interesting feature which can be observed in Fig. 7 is that the fluctuations on the mean value of the tripartite mutual information seem to decay with time. In order to quantify such decay, we fit the time scale of fluctuations (due to the oscillations around the mean) decay with a function of the form obtaining the following values of a, b in the case d C = d D = √ d for the GUE ensemble a = −3.9, b = 0.8, while for the Poisson distribution a = −16.3, b = 3.2. For GDE instead, there are no oscillations around the mean. It is possible to found similar behaviors in bipartitions of the system in which both d C and d D scales with the system size d, not necessarily being equals. When the system size of one of the partition is much greater than the other, the oscillations go to zero.
In order to show the connection between the Isospectral twirling and the Hilbert-Schmidt coherence, we define the dephasing superoperator D B (X) := i Π i XΠ i , where we introduced the spectral basis B = {Π i }, with Π i 's are rank-one projectors. This definition of coherence is natural, as it measures how the non-diagonal is an operator with respect to a certain basis defined by Π's. We define ψ U = U ψU † to be the state given a certain unitary transformation. Then, the l 2 measure of coherence in the basis B is given by In the case in which the state ψ is pure, then the first term is just one, otherwise, the quantity above is the purity of the initial state. Let us set tr (ψ 2 U ) = 1. The connection between the coherence and the Isospectral twirling of C B has been derived in detail in App. G.1, and we have the following expression where D B ≡ i Π ⊗2 i . It is interesting to see how the coherence depends on the dimension of the Hilbert space. The expression for the coherence for large d takes the following asymptotic value: As a function of time, a plot for d = 2 16 of the coherence is shown in Fig. 8 for the case of GUE, Poisson and GDE. We can see that while for t → ∞ the asymptotic value approaches the maximum value 1 in all the three cases, the convergence to equilibrium is characterized by the behavior ofc 4 (t); looking at Eq. (56), we note that the spectral coefficient c 4 (t) become negligible compared to 1 in a time O(1) (cfr. Table 2), for this reason we do not observe any scaling difference between the three ensembles, see Fig. 8.

Wigner-Yanase-Dyson skew information
In this section, we consider also upper bounds to the Wigner-Yanase-Dyson (WYD) skew information, which was introduced in [174,175]. The WYD has been historically studied in order to test how "hard" it is for an observable to be measured, but it has been used in a variety of contexts, including various generalizations of Heisenberg's uncertainty relations to mixed states [176]. The WYD is defined as It is easy to see that we can write the expression as Taking the Isospectral twirling of such a quantity we obtain: The Isospectral twirling is evaluated in App. G.2, where we see that the Isospectral twirlings of order 2 and 4 enter.
Recently, it has been found [163] that the Wigner-Yanase-Dyson skew information can be set to be a good measure of quantum coherence. Let B = {Π i } be a basis of the Hilbert space H, the coherence of ρ can be measured as: where D B (·) is the dephasing super-operator on B. It is worth noting that for pure states C(ρ) = C B (ρ) reduces to the previously defined l 2 measure of coherence. The Isospectral twirling of C(ρ) follows immediately from the Isospectral twirling of I η (ρ, X), see Eq. (58).

Quantum Thermodynamics
In this section we discuss some topics in quantum thermodynamics that result in a tell-tale of quantum chaotic behavior by means of the Isospectral twirling. First, we study the approach to equilibrium in a closed quantum system, following the setup of [95]. Then we turn our attention to quantum batteries. Among the quantum devices capable of performing work, quantum batteries are of central importance, also in view of the fact of the possible technological applications of nanoscale devices, and of possible quantum advantages in storing and extracting energy from the systems. Quantum batteries thus play a relevant role in the study of thermal machines and are an area of intense study [96,[177][178][179][180][181][182][183][184]. We apply the formalism of the Isospectral twirling to quantum batteries following the setup of [96] calculating average work and fluctuations. In the case of a closed quantum system, which we discuss in Sec. 3.5.2, the work is the amount of energy stored in the population levels defined of a certain Hamiltonian H 0 . The population levels can be changed via the interactions, defined via a certain operator V . In Sec. 3.5.3 we discuss the framework of open systems, explicitly calculating the free energy in a given bipartition of the Hilbert space H = H A ⊗ H B (and its convergence to equilibrium, following the setup of [95]). The application of the random matrix theory results of Sec. 2.2 to these quantities shows that the approach to equilibrium differs in chaotic and non-chaotic systems, and moreover that also the capability of storing and releasing energy in quantum batteries is affected by, and capable of discriminating, the quantum dynamics at play.

Convergence to Equilibrium
We are interested in how the Isospectral twirling enters the convergence to the equilibrium of a quantum many-body system. For this purpose, we follow the setup of [95], in which an equilibrium state ω in a bipartite quantum system As it is well known, in a closed quantum system there is no equilibration at the level of the wave-function. What is possible is equilibration in probability for subsystems. In other words, local observables attain typical expectation values [68] or the reduced density matrix to a local subsystem has almost always a very small trace distance from some typical state [119]. Consider the reduced state ψ A (t) = tr B ψ(t) of a global state evolving as ψ(t) = U ψU † and U ≡ exp{−iHt}. We can compute the distance of the state of the subsystem A, ψ A (t) from the equilibrium reduced state ω A = tr B ω, through the Schatten 2−norm [185]: It is preferable to study f (t) the square of the 2-norm instead of f (t) to keep the calculations simple. The square of the norm can be written as: Since the Isospectral twirling randomizes over the eigenstates of H, leaving the spectrum invariant, we need to calculate the Isospectral twirling of the dephasing superoperator D H ; in Sec. 3.6 we introduce the Isospectral twirling for CP maps: In order to take the correct Isospectral twirling for f (t), since we cannot consider the Isospectral twirling of U ⊗1,1 and the Isospectral twirling of D H to be independent, we need to define the following hybrid Isospectral twirling, which will be used ad hoc for the computation of f (t):R then we can take the Isospectral twirling of f (t): where we have defined: From Eq. (163) and Eq. (64) we obtain the isospectral twirling for f (t): The expression for large d are presented below:  Table 2).
The behavior of the convergence to equilibrium, for , as found in [95], GUE and Poisson equilibrate in a time O(d 1/12 ) and O(d 1/8 ) respectively. A plot is presented in Fig. 9. Equilibration in closed systems can also happen in integrable systems [68]. We see from our result that indeed the approach to equilibrium does not separate the behavior in GUE from that of Poisson. A similar analysis can be done for Eq. (67).

Quantum Batteries in Closed Systems
The model of quantum batteries we discuss in this section is that of a Hamiltonian H 0 which sets the scale of the energy via its quantum levels and of a quantum state ρ evolving in time as C t (ρ) = ρ t by means of a quantum evolution. The work on the battery, since the evolution is unitary and the entropy constant, results from populating the levels of H 0 in a different way from the initial state [186].
A quantum battery can be modeled by a Hamiltonian of the form H = H 0 + V (t). Let ψ t = U(ρ) = U ψU † the state at the time t obtained by unitary evolution U induced by the Hamiltonian H(t) [96]. The work extracted by (or stored in) the quantum battery is defined as the difference between expectation value of H 0 at time 0 and t: Following [96], one has that in the interaction picture the work can be written as where We can compute an average work by averaging over isospectral unitary evolutions. We obtain (see the App. H.1 for details) the expression where we have defined δW 0 ≡ (E 0 − E HT ), where E HT = tr (H 0 )/d and E 0 = tr (ψH 0 ). The fluctuations of work are defined by ∆ 2 W G := W 2 (t) G − W (t) 2 G and are also related toR (4) Define the normalized work fluctuations as The normalized work fluctuations in the large d limit are given by: where 1. An analysis of the work can be done, similarly to what has been done previously, on the basis of the dependence on the coefficient c 2 (t). In Fig. 10 we observe that the normalized work, defined as W (t) δW 0 (and normalized between 0 and 1), reaches the asymptotic value 1 in a time O(1); indeed, since the work depends on c 2 (t) we can see in Table 2 (left) that c 2 become negligible with respect to d 2 in a time O(1). We conclude that the average work does not distinguish chaotic from integrable dynamics. Conversely, the fluctuactions around the mean, see Eq. (75), are able to discriminate between chaotic and integrable behaviors.

Quantum Batteries in Open systems
An application for the Isospectral twirling is the study of quantum batteries in open systems. In the case of open systems storable and extractable energy does not correspond to extractable work [187], and thus in addition to the stored energy one needs also to consider the entropy [188], leading to the free energy. The free energy operator can be written as:  It is possible to apply the Isospectral twirling to the extractable work. In order to average with the Isospectral twirling we should take the average of the Von-Neumann entropy term, for this we choose instead to use the hierarchy of Rényi entropies to bound the Von-Neumann entropy consequently follows a bound for the extractable work, then we use the Jensen inequality to average and we obtain: where δW 0 ≡ (E 0 − tr (H 0 )/d) has been defined in Sec. 3.5.2. Defining ≡ δW 0 β, a dimensionless parameter depending on the temperature T of the bat, we get in the large d limit( see App. H.2): We see that the lower bound of the extractable work depends on thec 2 coefficient, while the upper bound onc 2 andc 4 . It is possible to observe that the two bounds differ for the term containing c 4 (t); this term become negligible after a time O(1) for all the three ensembles (cfr. Table 2). It is interesting to point out that when this term equals zero, the lower and upper bound equals each other and the extractable work is exactly determined and it reads: The normalized extractable work after the equilbration time O(1) ofc 2 E andc 4 E , will depend only from the temperature of the bath −1 = β −1 /δW 0 . A plot of the normalized free energy ∆ F G for the three different ensemble is reported in Fig. 11.

Isospectral twirling and CP maps
In this subsection, we briefly discuss the application of the formalism of the Isospectral twirling to completely positive (CP) maps. Given Q(·) = α K α (·)K † α a quantum map, i.e. a completely positive trace preserving map, with K α Kraus operators, the 2k-Isospectral twirling is defined as: where K ≡ α K α ⊗ K † α is a non hermitian bounded operator on H ⊗2 ; this operator has a peculiar property: taking T ∈ B(H ⊗2 ) the swap operator, then tr (T K) = d. We compute the Isospectral twirling for quantum maps with the usual techniques of the Haar average, see Sec. 4.1, and get: i.e. it is a linear combination of permutation operators weighted by coefficients depending on the spectral properties of the Kraus operators. The universality of the Isospectral twirling follows almost identically when it is applied to quantum maps. The 2−Isospectral twirling for quantum maps reads: see App. I for details. We see that it depends on the trace of K only; we will discuss the dephasing channel D B and explicitly compute the 2−Isospectral twirling for D B , see Eq. (92). Let us first give two examples of applications of the above formalism: the Loschmidt-Echo of the first kind L 1 , defined in Sec. 3.2.1, and purity Pur(ψ) = tr (ψ 2 ).
We define the Loschmidt-Echo of the first kind for quantum maps as: 2 where ψ ∈ B(H) is a generic quantum state. Taking the Isospectral twirling we get: see App. I; plugging the expression for the 2−Isospectral twirling Eq. (86) we have: Now consider the purity a state ρ = Q(ψ) with ψ ∈ pure states: and taking the Isospectral twirling we get: we see that it is proportional to the projector Π + = (1l ⊗2 +T )/2 onto the symmetric subspace of H ⊗2 , defined in Eq. (106); one important remark is worth making here: the asymptotic limit (t → ∞) for R (2) the above expression makes sense: the Loschmidt-Echo is quantifying the probability of finding the state D B (ψ) in ψ; after the action of the dephasing D B on a random basis, the state becomes highly mixed; indeed, for the properties of the projectors of rank 1, we see that the purity in Eq. (90) for is equivalent to the Loschmidt-Echo; D B requires the 2−Isospectral twirling instead of the 4−Isospectral twirling as the other quantum maps, see Eq. (90). We obtain the average coherence on a random basis B, see Sec. 3.4.1 for the definition: i.e. given a pure state ψ the average coherence on a random basis is nearly maximal in the large d limit. We remark that the formula Eq. (95) can be regarded as the average coherence of a random pure state ψ on a fixed basis B, for this reason, an analog result can be found in [164,169].

Isospectral twirling: techniques
This section is devoted to the description of all the techniques we used to derive the results in Sec

Haar average and Weingarten functions
The goal of this section is to calculate Eq. (2). The Haar average over the unitary channel can be performed via the following well known procedure [189][190][191][192], which we briefly review here for completeness. Let A ∈ B(H) a bounded operator on the Hilbert space H and G ∈ U(H) a unitary operator on H chosen uniformly at random. The Haar average of A ⊗k on the unitary group reads: where dG is the Haar measure over the unitary group; it has the following properties dG = 1 and dG = d(GG ) = d(G G) for any G ∈ U(H) unitary on H; the latter is called left(right)-invariance of the Haar measure. The general formula to compute the Haar average is [23]: where π, σ are permutations of the permutation group S k while W g (πσ) are the Weingarten functions, and πσ is the product of the permutation π and σ, see Sec. 4.2 for details on permutation operators. The Weingarten functions [190] are defined as where the sum runs up to the number of conjugacy class of S k . χ λ is the irreducible character of S k associated with the conjucacy class λ of S k , e is the identity element of the group. s λ (1) ≡ s λ (1, . . . , 1) is the Schur function evaluated on d = dim(H) elements. While the definition of the Weingarten functions is given above, we present an easy way to compute them due to Roberts and Yoshida [23]. Any permutation operator T κ ∈ S k satisfies [A ⊗k , T κ ] = 0, where [·, ·] is the commutator between operators of the algebra B(H ⊗k ) of the bounded operators on H ⊗k ; therefore taking the average of the permutation operator T κ , from Eq. (96) one has: inserting the above condition in (97), we get: therefore if we define Ω a k! × k! real symmetric matrix with components: we can simply express the Weingarten functions, treated as a matrix with components W g (πσ), as the inverse of Ω, i.e: the above is a straightforward way to compute the Weingarten functions: we calculate traces of all the products T π T σ and build the matrix Ω, then we calculate the components of the inverse matrix (Ω −1 ) πσ .

Computation and properties of the Isospectral twirling
In this subsection, provided with tools and techniques of the Haar measure, we compute the Isospectral twirling in Eq. (2). Here we explicitly present the simplest case of Isospectral twirling, i.e k = 1, for k = 2 see App. A. For k = 1 we have only two permutation operators because the permutation group S 2 contains only two elements: the identity permutation 1l ⊗2 and the swap operator T . The matrix Ω, introduced in Eq. (101), and its inverse read: π (U ) = C d (tr U α ) γ (tr U †β ) δ , the corresponding powers α, β, γ, δ and the constant C d . See Eq. (7) for the short notation c 2 (t). Table 3:

From Eq. (3) we obtain, see
we see thatR (2) (t) depends only on c 2 (t), which we introduced in Eq. ( 7).R (2) (t) is a timedependent hermitian operator; for t = 0 it equals the identity, while the asymptotic value t → ∞ for the ensemble averages we consider reads (cfr. Table 2): i.e. it is proportional to the projector Π + = (1l ⊗2 + T )/2 onto the symmetric subspace S + of H ⊗2 , defined as: from Eq. (92) we can also write the following relation: whereR (2) (D B ) is the Isospectral twirling applied to the dephasing superoperator, see Sec. 3.6 for the definition; the reason is twofold: first, it is easy to see that the infinite time average of E T U ⊗1,1 ∝ D B and we prove in Sec. 4.5 that the asymptotic value for the ensemble average ofR (2k) coincides with its infinite time average: more precisely we prove that the infinite time average and the spectral average with Schwartzian probability distributions do commute. If one, instead of considering the Isospectral twirling of a unitary evolution U (t) generated by an Hamiltonian H, compute the Haar average of the Isospectral twirling in Eq. (2) with respect to U , one obtains: where the last equality follows from the left-invariance of the Haar measure. The general result for Eq. (108) can be found in Corollary 3.21 of Ref. [193]: where π, σ ∈ S k and T π −1 ,σ |i 1 , . . . , i k , j 1 , . . . , j k = j σ(1) , . . . , j σ(k) , i π −1 (1) , . . . , i π −1 (k) . For k = 1, this expression reduces to: where T is the swap operators. Looking at Eq. (104) we see that the closer c 2 (t) gets to 1, the closer the Isospectral twirling for k = 1 is to its Haar average value.

Algebra of the permutation operators
For the sake of clarity, here we would like to briefly introduce the cycle representation of the permutation operators and their algebra used ubiquitously in the text. First, the permutation operator T π corresponding to the permutation π ∈ S k reads: These are generalization of the swap operator. Since the permutation operator acts on the right to left rather than left to right. In this paper, we are mostly concerned with permutation elements in the groups S 2 and S 4 . Since we are using the cycle representation of the permutation group, if the element is (23) in the group S 4 , clearly we are using the standard shorthand notation for (1)(23)(4), meaning that only the channels 2 and 3 are being swapped. The permutation operators thus follow the multiplication rules of the permutation group. For in the permutation group S 4 , (12)(23) = (123), and thus T (12) T (23) = T (123) . An intuitive and graphical representation of how the permutation operators multiply is shown in Fig. 12 (left). The permutation operator can be represented as a series of interchanging lines. For instance T (12) can be represented as the interchange of the first and second channel. Thus if we combine T (12) and T (23) , their product results is an interchange between the first and third channel, and the second and the first and the third and the second. This approach can also be used to calculate traces. Since we also have traces of permutation operators in S k combined with operators in H ⊗k , the traces can also be performed graphically. The graphical permutations of the channels are then combined with boxes, representing the action of operators on the channel. A trace is then represented as a closed sequence of lines that connects the final and initial channel in the same position. For instance tr (T (324) (A ⊗ B ⊗ C ⊗ D) = tr (A)tr (BCD), is shown in Fig. 12 (right). The global trace connects the nth line of the input to the nth line of the output. In the example above, one then has a continuous line in which one can read the line as a trace of A multiplying the trace of B, C and D. Two properties which we use a lot throughout the paper are the following; let T π , T σ ∈ S k two permutation operators with π be a generic permutation and σ = (1 k k − 1 . . . 2); let A i ∈ B(H) with i = 1, . . . , k be bounded operators on H. Then: the proofs are given in App. C.

Spectral average over eigenvalues distributions: GDE and GUE
In this section, we perform the spectral average of the spectral functions c a (t), a = 2, 3, 4, in the two ensemble of random matrices, GUE and GDE. Let us report Eq. (11) and Eq. (12) for convenience: the calculations for GUE were previously performed in [94,194] and we discuss the technique here, see Sec. 4.3.1, for completeness. We compute the spectral average over the GDE distribution in Sec. 4.3.2; similar results can be found in [95].

GUE
In this section, we present the calculation of the spectral average of c 2 (t) for GUE Hamiltonians, while the calculation for c a (t) for a = 3, 4 are reserved to App. B. By definition Eq. (9) the spectral average Figure 12: Graphical representation of the swap operator algebra and how the traces work. As a notational note, since the swap operator is defined as T σ = ijkl |σ(i)σ(j)σ(k)σ(l) ijkl|, the permutation group element in the cycle representation is intended to work from the right to the left.
of the coefficient c 2 (t) reads: First of all, note that since E i , E j are dummy variables, the sum reduces to the counting i = j, while for i = j the integral sum to 1 thanks to the normalization of the probability distribution. In the third equality in Eq. (115) we have the 2-point marginal probability distribution. The n−point marginal probability distribution is defined as: ρ it can be expressed in terms of the determinant of the kernel K [8-10]: K is a n × n matrix with the following matrix elements in the large d limit: the matrix elements K ij for i = j are the well known Wigner semicircle law [3]. The calculation of c 2 (t) requires the 2−point correlation function: where the θ−function is defined as θ(x) = 1 for x ≥ 0, θ(x) = 0 for x < 0 and J 1 (x) is the Bessel function of the first kind. A plot of c 2 (t) GUE is reported in Fig. 13(Left), while the scalings can be found in Table 2.

GDE
In this section we calculate the spectral average of c 2 (t) for GDE Hamiltonians. From Eq. (12) we immediately see that the probability distribution P GDE factorizes, i.e therefore the calculations for the spectral averages reduce to combination of one dimensional Fourier transforms of Gaussians; the marginal probability for the GDE distribution is trivial since the eigenvalues are uncorrelated from each other. Let us compute the spectral average for c 2 (t): The calculations for c a (t), a = 3, 4, can be easily made following the same techniques, see App. B.
A plot of c 2 (t) GDE can be found in Fig. 13(Right).

Nearest level spacing distribution: technique
In Sec. 2.2 we claimed that the spectral functions c (2k) π (U ) are functions of the nearest level spacing distribution P (s); in this section we illustrate the technique used to average over families E of Hamiltonians with a given nearest level spacings distribution P E (s). The nearest level spacings distribution we consider are listed in Eq. (14) and reported here for convenience:

Wigner-Dyson from GUE
We use the same techniques of [99], i.e. we define the characteristic functions of the probability distributions P E (s); a calculation for the spectral average of c 2 (t) for Poisson recently appeared in [195]. Here we compute c 2 (t), while the details of c a (t) for a = 3, 4 can be found App. B. The gaps between the energy levels are there are d(d − 1)/2 gaps, but only d − 1 are independent from each other. Indeed, defining the nearest level spacings as: we can express the gaps ∆ ij Eq. (124) for any i, j as a linear combination of s α , namely: To compute the spectral average of c 2 (t) with respect to an family of Hamiltonians E with a given nearest level spacings distribution, we first need to express c 2 (t) in terms of s α : Let us use the expansion of the cosine in terms of exponentials: The next step is to take the spectral average of c 2 (t). In order to make the average of c 2 (t) Eq. (129) tractable, we use the hypothesis of molecular chaos, that is, the stosszahlansatz [196], for which we consider s α 's to be independent identically distributed stochastic variables; let us define the operation formally: let f ({s α }) be a function of the nearest level spacings, we compute its spectral average as: this hypothesis works well for Integrable Hamiltonians and therefore for E ≡ P [8,104]. Given a probability distribution on the nearest level spacings P E (s) we define the characteristic function of the level spacing s α [99]: which, thanks to the stosszahlansatz, does not depends on α. Taking the spectral average of Eq. (129), we get: Making the sum with the usual techniques, one has the final result: The characteristic functions for Poisson, Wigner-Dyson from GOE and GUE read:

Normalization of the energy scale.
Throughout the paper we compare the spectral average of c a (t) for various ensembles of Hamiltonians.
To this aim we need to compare Hamiltonians whose largest eigenvalues scale in the same way with the Hilbert space dimension d. In Sec. 2.2, we have defined GUE and GDE ensembles such that the largest eigenvalue is O(1). We need to normalize the energy scale of the nearest level spacings distributions; indeed, we claim that, in our settings, the largest eigenvalues for the nearest level spacings distributions is O(d); we give the following argument: compute the average maximum spacing ∆ d1 : where we used the fact that s = 1; indeed the distributions in Eq. (14) are normalized such that dsP E (s) = 1 and s = ds P E (s)s = 1. This implies the maximum eigenvalue to be O(d): without loss of generality choose E 1 = 0: i.e. the nearest level spacings distribution for H and H − E 0 1l is the same. The normalization of the largest eigenvalue can be done by the replacement t −→ t/d in the expression of c a (t).

Envelope curves
In this section present the envelope curves of the spectral average of c 2 (t) and c 4 (t), for Poisson and GUE ensembles, from which the scaling behavior in Table 2 can be computed. We evaluate the envelope curves using the method of [94], i.e. we calculate the asymptotic behavior of c 2 and c 4 for large d and then we suppress the oscillating part. Using this technique, the envelope curves forc 2 (t) E , for E = P, GUE are given by: where we define: the envelope curves forc 4 (t) E are theñ In order to study the timescales of GUE up to the dip time, as claimed in [94], it is sufficient to consider simpler expressions: c 2 GUE ≈ 1/πt 3 and c 2 GUE ≈ 1/π 2 t 6 . Let us give some example of the calculations for the timescales presented in Table 2; let us compute the equilibrium time for Poisson, see Eq. (137). We define the equilibrium time t p as the time for which the function 2/t 2 become comparable with the factor 1/d and therefore we need to solve the following equation for t: The equilibrium time for GDE can be obtained with a similar argument from the exact expression of the spectral average, see Eq. (122):

Asymptotic value collapse
In this section we give the mathematical motivation for which the asymptotic value of GUE and GDE are the same. The reason lies on the properties of the spectral distributions. Indeed here we prove that if the spectral probability distribution P ({E i }) is Schwartz, namely it decays to zero faster than any polynomial of generic degree, then the asymptotic value of the spectral average does not depend on the particular choice of P ({E i }). Proposition. Let E to be an ensemble of isospectral Hamiltonians characterized by a Schwartz probability distribution P (E) with E ≡ (E 1 , . . . , E d ), then the spectral average of the coefficients c (2k) π (U ) reads: where f π (d) depends on the permutation element π and on the dimension of the Hilbert space H, but on the particular choice E. Proof. In general, the spectral function c (2k) π (U ) can be written as: where α, β, γ, δ and C d are assigned to π according to Tables 3 and 4 and the overall sum contains d γ+δ elements. From decompositions in Eq. (169) and Eq. (209) is clear that the these coefficients split in two parts: one independent from t (which is the asymptotic value f π (d)) and one dependent from t. This decomposition does not depend on the particular choice of the spectral probability distribution, but on the powers α, β, γ, δ. Indeed, because of the linearity of the exponent in Eq. (145) with respect to the energies E i , we can recast it in the form of a scalar product in R d : where we have defined vectors C j 1 ,...,j δ i 1 ,...,iγ ∈ R d , one for each element of the sum in Eq. (145); their components depend on α, β and take into account the coefficients of the linear combination of the E i s. Defining the multi-index a = (i 1 , . . . , i γ , j 1 , . . . , j δ ) ∈ [1, d γ+δ ] to clean the notation, we can rewrite Eq. (145) as: let us split the sum in two parts: C a · E = 0 and C a · E = 0: where we used δ(x) = 1 iff x = 0 and we defined: we re-labeled as C b the first D = d γ+δ − f π (d) vectors satisfying C b · E = 0. Taking the spectral average (see Eq. (9)) of the r.h.s. of Eq. (148) we finally get: whereP (x) is the Fourier transform of P (E); taking the limit t → ∞ of both sides of Eq. (150) one get the desired result. Indeed the Fourier transform of a Schwartz function is a Schwartz function and dies asymptotically. We remark that the asymptotic value for the ensemble average, Eq. (144), coincides with the infinite time average, defined in Eq. (26), of c (2k) π (U ); indeed, from Eq. (147), we have: it is well known that the above integral is always null, but for C a · E = 0; therefore: i.e. we proved that the ensemble average for Schwarzian probability distributions and the infinite time average do commute.

Lévy lemma and typicality
In this section we apply the Lévy lemma on P O (G) = tr (T (2k) OG †⊗2k U ⊗k,k G ⊗2k ), where G ∈ U(H) and prove the bound Eq. (14). The Lèvi lemma states that if f : G −→ C is η−Lipschitz, then [129]: where the average f (G) G is taken with the Haar measure over the unitary group. This statement is extremely useful in the context of many body physics; the Hilbert space dimension d, appearing at the exponent in Eq. (153), grows exponentially with the system size, therefore, if the Lipschitz constant η does not scale with d, the probability of finding a value f (G) different from its average f (G) G is double exponential small in the system size. This concentration bound is way stronger than the Chebyshev inequality. Let us prove the bound in Eq. (14).
Consider the following series of inequalities: The first equality follows from the linearity of the trace. The first inequality follows from |tr (AB)| ≤ A p B q , where p −1 + q −1 = 1, · p is the Schatten p-norm [185] and from AB p ≤ A p B p ; the second inequality has been proven in [128], we include here the proof for completeness: consider A, B, C, D to be unitary operators, then indeed for any unitary operator U , U ∞ = 1. By iterative application of the latter inequalities one get: The third inequality in Eq. (156) is obtained by similar techniques: The last inequality in Eq. (156) follows easily from the hierarchy of the Schatten p−norm A p ≤ A q if p > q. This concludes the proof.

Conclusions
Quantum chaos is a property of quantum dynamics that should, on the one hand, feature a very diverse list of properties and high sensitivity to the details of interactions -the butterfly effect -complex growth of entanglement and coherence, information scrambling, and rapid decay of the out-of-timeorder correlation functions, and, on the other hand, tell apart the dynamics generated by an integrable Hamiltonian, which should not be chaotic from that of non integrable Hamiltonians, that are supposed to be chaotic. These properties are usually understood in terms of the statistics of the gaps in the energy spectra, resulting in Poisson or universal distributions. In this paper we have unified in a single framework the description of the numerous probes to quantum chaos in a way that allows for a clear classification in terms also of the spectral properties. The (unitary) evolution operator belongs to a family of isospectral operators with different eigenvectors. By averaging over these eigenvectors, one obtains quantities that only depend on the spectrum. The isospectral twirling indeed consists in practice in averaging over the eigenvectors of a given quantum dynamics. By performing calculations in random matrix theory for the appropriate ensembles of integrable or non-integrable Hamiltonians, we showed how the unified probes to quantum chaos neatly separate chaotic from non chaotic behavior in terms of how their values and corresponding time scales scale with the dimension of the Hilbert space d. We also showed why the infinite time, large d values do not depend on the type of Hamiltonian chosen and why they detach from the prediction from the Haar measure.
In perspective, a number of open questions can be explored using the several techniques used in this paper. We would like to understand how to incorporate the locality of interactions and what effects this has on quantum chaos. We showed how to extend the isospectral twirling to general completely positive maps. In perspective it would be interesting to understand in what sense these channels or master equations described by Lindblad operators can be chaotic. Finally, a very interesting subject is the transition from an integrable dynamics to a chaotic one. In the context of quantum circuits, it has been understood that one can dope a random Clifford circuit and obtain higher order unitary t-designs or transitions to universal entanglement [130,197]. In the context of quantum many-body systems, one could study integrable spin chains perturbed by an integrability-breaking term [142] and study the typical behavior in the sense of the isospectral twirling of the probes to quantum chaos. The insights gained in this way could open the way to formulate a quantum KAM theorem, that is, understand to what extent a chaotic system obtained by perturbing an integrable one [17] shows remnants of its lost integrability. Perhaps gazing enough into quantum chaos will have it gaze back at us.
In order to write an explicitly the 4−Isospectral twirlingR (4) (U ), we define the following elements: ≡ T (1324) + T (1423) Y (2)(2) ≡ T (14)(23) + T (13) (24) it is now possible to insert the c U π in Eq. (160) and then the final formula forR (4) (U ) reads: is a time-dependent operator asR (2) (U ); for t = 0 it equals 1l ⊗4 , while the asymptotic limit t → ∞ for the ensemble average we consider is (cfr. Table 2): In Sec. 4.5 we proved that the infinite time average and the ensemble averages commute; for this reason we can express Eq. (163) in terms of the isospectral twirling of the dephasing superoperator. First, consider the infinite time average of U ⊗2,2 : now, taking the Isospectral twirling in the sense of quantum maps, defined in Sec. 3.6, we get: where: where Π + 1 (2) is the projector onto the symmetric subspace of the first(second) two copies of H, while Π + (4) = π T (4) π /24 is the projector onto the symmetric subspace of H ⊗4 . It is possible to compute the Haar average of the Isospectral twirling with respect to U for as shown in Eq. (108) for k = 2, from Eq. (109) reads: Looking at Eq. (163), it is possible to see that if c 4 (t) = c 2 (2t) = 2,c 2 (t) = 1 and Re c 3 (t) = 0 we obtain the Haar value.

B Higher-order spectral functions
In this section, we compute spectral averages for the higher-order spectral functions c 3 (t) and

B.1 Spectral function c 4 (t)
In this section, we develop the calculation of the spectral function c 4 (t) introduced in Eq. (7), before to move on the spectral average for each distribution, we give a decomposition of c 4 (t), that will help us in the calculation of the spectral average. First, we start recalling the definition of c 4 (t): it is possible to decompose the sum in Eq. (168) as: that can also be written in terms of cosines as: B.1.1 Calculation of c 4 (t)

GDE
The calculations for the spectral function c 4 (t) GDE repeat similarly to what done for the spectral function c 2 (t) GDE , we decompose first as Eq. (170) and after we take the average over the GDE ensemble, we obtain the following terms: The final expression for c 4 (t) GDE is: the time for which it achieves the equilibrium is the same as for c 2 (t) can be found in Fig. 15.

B.1.3 Calculation of c 4 (t) E for nearest level spacing
In this section, we calculate the coefficient c 4 (t) E for the nearest level spacing distribution, and we present the plot of the averaged spectral function c 4 (t) for the Poisson distribution, Wigner-Dyson GUE distribution and Wigner-Dyson GOE distribution. We use the same technique presented in Sec. 4.4. The calculation for c 4 (t) E start decomposing c 4 (t) as in Eq. (170) and subsequently we take the average over the nearest level spacing distributions Taking the average over the ensembles E, we have four dynamics terms, two of them are already evaluated in 133: We first deal with the average of i =j =k =l cos( It is possible to point out that there are 4! ways to order four items, but with the conditions, i > j and k > l they can be reduced to six, the following ones: noting that there is another symmetry, that is the transformation i → k and j → l. Due to this symmetry, we have only three remaining different orders: Therefore, we can split the sum into 3 pieces: The first term S 1 (t) given by the combination (i > j > k > l) is: cos(s l +· · ·+s k−1 +2s k +· · ·+2s j−1 +s j +· · ·+s i−1 ) (181) where ∆ il are the gaps between the energy levels. An equal result could have been obtained using cos(∆ ik + ∆ jl ). The cosine reads: cos(s l + · · · + s k−1 + 2s k + · · · + 2s j−1 + s j + · · · + s i−1 ) = cos( Recalling the definition of g E (t) in Eq. (131) and making use of the hypothesis of molecular chaos discussed in Sec. 4.4 for which s i can be treated independently from each other let us take the ensemble average: The result is: The term S 1 (t) becomes: The second term S 2 (t) corresponding to the case(i > k > l > j) read as: Therefore the cosine can be rewritten as: Following the same procedure as before we obtain: The third term S 3 (t) corresponding to the case(i > k > j > l): then the cosine can be written as: that becomes: This results allows us to write S 3 (t) The final result is: The last term that we have to study is , we can divide this sum into three parts: i > j > k, j > k > i and j > i > k: It can be exploited and we obtain that v 1 (t) is: The cosine reads: cos (s k + . . . + s j−1 + 2s j + 2s i−1 ) = cos Since, all the terms in the cosine are independent from each other we have The result is: For the case j > k > i: Making use of the cosine property we have that We can proceed as before and we obtain For the case j > i > k: The cosine can be rewritten as Just as before we obtain: The final result reads: It is possible to combine all the obtained result and obtain the spectral function c 4 (t) for a nearest level distribution. We give a plot of the results in Fig. 16.

B.2 Real part spectral function c 3 (t)
In this section, we are going to calculate the real part of the spectral function c 3 (t), since in the calculations of our probes of chaos only the real part enters. The calculations are similar to the ones made for c 4 (t), we start recalling the definition of c 3 (t) introduced in Eq. (7).
that can be decomposed as: , for this reason in the literature [94,194] the 3−point coefficient can appear labeled by 4, 1, where the 1 labels the contraction. The average of the real part of the coefficient c 3 (t) can be written as: If we decompose c 3 (t) as done in Eq. (209) and remember the definition of ρ (n) (E 1 , . . . , E n ) given in Eq. (116) Re the integrals are the same of the c 4 (t) GUE . The final result is: A plot of the Re c 3 (t) GUE for the GUE ensemble can be found in Fig. 17(Left)

GDE
In order to calculate c 3 (t) GDE , we first decompose as in App. 209, after we take the average over the GDE distribution Eq. (12) and making use of the calculations made in Eq. (174), we obtain: The asymptotic value is: lim A plot of the spectral function c 3 (t) GDE can be found in Fig. 17(Right).

B.2.3 Calculation of c 3 (t) E
In this section we calculate c 3 (t) E . As done before we decompose with Eq.209 and the average of Re c 3 (t) is given by: It is possible to recognize that the first two dynamic terms are the ones obtained in and are estimated in Eq. (177), while the last term is calculated in Eq. (207). The results are plotted in Fig. 18: The equilibrium value of Re c 3 (t)

C Auxiliary results for algebra of permutations
The aim of this section is to derive the property introduced in Eq. (112) and (113)

C.1 Permutations of operators
In order to prove the property of Eq. (112) we have two first to prove the result for a swap operator and then generalize it. For this simple case, we aim to prove that T (12) The action of T (12) is This means that we can rewrite Eq. (218) as: This concludes our proof. This proof can be generalized for a generic permutation T π ∈ B(H ⊗k ), regarding T π as a product of swaps, between H i ⊗ H j , with i, j ∈ 1, . . . , k.

C.2 Traces of order k
In this section we want to give a proof of Eq. (113). We have to proof the following statement: where A 1 , . . . , A k ∈ B(H). The first term can be written as Analyzing the second term, we have (223) The action of T (1 k (k−1)...2) on the product state is It is possible to rewrite Eq. (223) as This concludes our proof. The definition of the frame potential is [23,126]: The frame potential can be rewritten as: It is possible to employ the equivalence (tr A) k = tr A ⊗k , the frame potential becomes: where we have defined U ⊗k,k ≡ U ⊗k ⊗ U †⊗k . The trace operator and the integral commute, and employing the definition Eq.

D.2 Proof of the bound Eq. (25)
A proof for the bound Eq. (25) we make use the following bound [185]: The 2-norm of 2k-fold Haar channel of a unitary operator can be bounded as: The property rank(A) = rankA † A allows us to write rank(R (2k) (U )) = d 2k , this means: In order to proceed we employ the following bound [185]: where p −1 + q −1 = 1, if we set B ≡ 1l, p = 1 and q = ∞, we obtain |tr A| ≤ A 1 , from which follows: The obtained bound, for the frame potential becomes: This concludes the proof.

D.3 Proof of Eq. (29)
From the bound Eq. (25) we have that: using the decomposition Eq. (147) proven for the Proposition in Sec. 4.5, we can write the above expression as: for some C a s. Taking the spectral average of the above expression with respect to the GDE ensemble P (E), one has, see Sec. 4.5: the spectral distribution for the GDE ensemble factorizes in Gaussians (cfr. Eq. (12)); the Fourier transform of a Gaussian is a Gaussian in the conjugate space and therefore we haveP (Ct) ≥ 0 ∀t ≥ 0. This concludes the proof.

E OTOC and Loschmidt-Echo
The definition 4−point OTOC is: with A, B ∈ B(H) two local operators. The introduced definition can be written as: The product can be linearized thanks to the replica trick C.2 as: This can be rewritten as: (23) (242) where we used that T 23 T 23 = 1l. Thanks to the property Eq. (112) it is possible to write: The adjoint action of T (23) on T (1432) is equal to substituting and making the Isospectral twirling we obtain: where A = A † ⊗ A, and B = B † ⊗ B.

E.2 Loschmidt-Echo of second kind
The definition of Loschmidt-Echo of the second kind is: where we introduced a small perturbation δH of the Hamiltonian H. If Sp (H) = Sp (H + δH), we can write: with A a local unitary operator close to the identity. Eq. (247) and the unitarity of A help us to rewrite the LE as: For the Loschmidt-Echo writing U G := G † U G: In order to get A ⊗A = A † ⊗A⊗A † ⊗A, we act adjointly with T (12) , noting that [T (12) , (U ⊗2,2 G )] = 0: Taking the average over G we obtained the desired result.

E.3 Estimations for 4−point OTOC and second kind of Loschmidt-Echo
Let us now consider A, B to be non-overlapping Pauli operators on qubits and let us make this calculation explicit. For Pauli operators on qubits, the OTOC and the LE read: In this settings, we have tr (A) = tr (B) = tr (AB) = 0, tr (ABAB) = d and tr (A 2 ) = tr (B 2 ) = d. Computing the traces for both the quantities the results read: The entanglement of ψ t in the given bipartition computed by the 2−Rényi entropy is given by where we used the identity tr A 2 = tr (T A ⊗ A). In order to simplify the notation we introduced the the operator T (A) defined as T (A) ≡ T A ⊗ 1l ⊗2 B . This means that the entanglement can be written as As usual, we want to linearize the product in the trace, in order to do this we set a = U ⊗2 ψ ⊗2 , while b = U †⊗2 T (A) and we use tr (ab) = tr T (13)(24) a ⊗ b and we obtain We can map U → G † U G and average over G. The result is where The lower bound in Eq. (260) is obtained thanks to the Jensen inequality and we used the equivalence U ⊗2,2 G =R 4 (U ). Using the explicit form ofR (4) , see Eq.( 162), we find If our state ψ is pure we obtain A plot of the r.h.s of Eq. (262) can be found in Fig. 6. The large d expressions are: The equilibrium value for the case d A d B is It is possible to map isomorphically an operator O ∈ B(H) into a state |O ∈ H ⊗2 , the isomorphism between an operator and a state is called Choi isomorphism. If we apply the map to the unitary operator U introduced above, we obtain the normalized state: where we labeled the states with 1(2) the two copies H ⊗∈ . The state |U can be rewritten as: where |I is the Bell-state on H ⊗2 defined as In what follows we use the density matrix associated with this state: It is interesting to point out a property of the Choi state if we trace the input or the An important property is that if one traces out the input or the output, the resulting state be maximally mixed.
The above property can be generalized for any unital CP −map f acting on the output ρ f = 1l ⊗ f (|I I|), one can prove that: where the labels in and out labels the two copies of H. The proof is first developed for input. The r.h.s of Eq. (274), can be written explicitly written as: since f is unital the proof for the input is concluded. We have now to prove the second statement, the one involving the trace of the output. The action of a unital map can be written as where U i ∈ U(H) and i c i = 1, c i ≥ 0, this conditions assure that f is trace preserving. From this follows: The last equation concludes the proof.

F.2.2 Tripartite Mutual Information
This section aims to provide a brief introduction on the tripartite mutual information, as pointed out in [64] and then in [125], the tripartite information can be considered a measure of scrambling. The simplest model available is to consider a unitary time evolution U AB→CD defined on two fixed bipartitions A, B the input bipartition and C, D the output bipartition, with the constraints that In this set-up, the definition of tripartite mutual information is: in terms of Von Neumann entropy the TMI reads, [64,88]: It is possible to use the property for the sum of entropies to obtain It is possible to prove from the above definition and properties that S(C ⊗ D) is equal to log d for the Choi state and therefore: This quantity given is a constant minus a sum of two entropies. All these entropies are calculated, from the Choi state ρ U associated to U . Now, we condensate all the results proven in [125] that we need in the following. First, the tripartite mutual information has the following bounds: so it's a negative-definite quantity. The more negative is, the more scrambled the information will be. The upper bound is achieved, by all and only the unitaries made as: where A L , A R , . . . , D L , D R are sub-partition of A, . . . , D respectively, see [125].

F.2.3 2-Renyi TMI
Since it is not easy to work with the Von Neumann entropy, we focus our attention on the 2-Renyi entropy, as previously done in [64]: The choice of the two Renyi entropy as a substitute of the Von-Neumann entropy gives us an upper bound for the mutual information. Let's avoid the subscript ρ U and write I 3 (2) as where ρ AC = tr BD ρ, while ρ AD = tr BC ρ, we have: From one hand, indeed, the quantity we are going to deal with is not exactly the TMI, but from the other hand, if one is interested in unitary evolutions that scramble the information, the 2-Renyi TMI (2RTMI) is sufficient. Moreover, it's easy to convince ourselves that the TMI and the 2RTMI share the same bounds: Proposition The unitaries of the type U C ⊗ U D , where U C ∈ U(H C ), U D ∈ U(H D ) achieve the upper bound of the 2RTMI and therefore, according to this measure of scrambling, do not scramble the information.
Proof Firstly, rewrite Eq. (283) as where T (A) ≡ T A ⊗ 1l BCD and T A a swap operator with support on H A similarly for T (C) and T (D) . Now, we calculate the scrambling give by: and insert it in Eq. (286) and we obtain Now making use of the cyclic property of the trace It is possible to observe that for the second term U ⊗2 D does not act on T (C) , so it cancels with U † ⊗2 D , similarly for the third term but with U ⊗2 C , T (C) = 0, and U ⊗2 D , T (D) = 0, we obtain at this point a straightforward calculation shows that I 3 (2) (U C ⊗ U D ) = 0. This is not true if we employ the Von Neumann entropy as proved in [125], as explained in the first proposition.
Proof of Eq. (49) Starting from Eq. (283) we can rewrite it as Eq. (286): where ρ = (1l AB ⊗ U CD ) |I I| (1l AB ⊗ U † CD ). Let's focus on the second term alone, the third one is similar: In order to proceed with the proof we need two facts. First of all note that |I I| = ψ AC ⊗ψ BD where ψ AC , ψ BD are Bell states. Then, note that the following identity holds: let T be the swap operator and let ψ ∈ H a pure state: Now turn again to Eq. (293): first we use the above identity so we obtain: then we can trace out H A ⊗ H B and Eq. (293) becomes: where the factor d −2 comes from the trace on the Bell state: the proof repeats identically for the second term Eq. (292). This concludes the proof.
Isospectral twirling of 2RTMI If we write the complete form ofR (4)  Define U G = G † U G. The coherence in B of a state evolved by U G can be written as where we used the trick tr ab = tr (T a ⊗ b) calling a = Π i ⊗ Π i U ⊗2 G , b = ψ ⊗2 U †⊗2 G . Now, taking the average over G by twirling we obtain Using the explicit form ofR (4) , see Eq.( 162), we find [ d 2 (3 + d)(d − 1) 2 + 2 (c 4 (t) + 2 Re c 3 (t) + c 2 (2t) − 8c 2 (t)) + (d + 1)tr (D B ψ) 2 (4c 2 (t) − c 2 (2t) − 2 Re c 3 (t) + c 4 (t)) ] The expression of C B (ψ U ) G for large d becomes: It is interesting to observe that tr (D B ψ) 2 ) is the purity of the dephased initial state and for t = 0 we have the coherence in the basis B of the initial state. The equilibrium value is in the large d limit: Since we are interested to see how the evolution affects the coherence, we choose ψ as one of the rank-one projector of B, Eq. (302) becomes The result is plotted in Fig. 8 in the main text.

H.1.2 Work fluctuations
In this section, we aim to develop the calculations for the work fluctuations. We start from the definition of fluctuations If we analyze the term W (t) 2 and make use of the property Eq. (113) and we can write: To see the correct order, let us insert T 2 (23) = 1l and make use of the property shown in Eq. (112): Now we take the Isospectral twirling of (K ⊗2 ⊗ K † ⊗2 ): The fluctuations of the work are given by Making use of Eq. (317) we obtain Making the hypothesis that ψ is a pure state and that it is an eigenstate of H 0 , the work fluctuations becomes: For large d the fluctuations are given by It is interesting to point out that for t −→ ∞, the ensemble averages of the Isospectral twirling of the work fluctuations for the studied distributions converge to the same value that is Let us define the dimensionless control parameter h = 4E HT E 0 /(tr (H 2 0 )/d − E 2 HT ) and study: where ∆W 2 (t) G = ∆W 2 (t) G /(tr (H 2 0 )/d − E 2 HT ). In large d limit Eq. (337): the results are plotted in Fig. 10 in the main text.

H.2 Quantum batteries in Open systems
This section will expand the calculations made in Sec. 3.5.3, we begin recalling the definition of free energy operator is Eq. (76)F defined on a Hilbert space H A ⊗ H B . The A subsystem represents the system from which one extracts work and the B is the bath of the system with temperature T = 1/(k B β). The reduced evolution on A is not unitary and reads ψ U,A = tr B (U G ψU † G ) and is induced by a Hamiltonian of the form H = H A + H B + H AB . The extractable work is defined as where we have defined ψ U,A = tr B U ψU † = tr B ψ U . Let us manipulate this quantity: If we use this identity tr (G(A ⊗ 1)) = tr A (tr B (G)A) we obtain: It's possible to recognize the Von Neumann entropy in the second term of the r.h.s. We can regard it as a Renyi entropy S α , defined as: with α = 1 is: and use the hierarchy of Renyi entropies: S α (ψ) ≥ S α (ψ), α < α ∈ N 0 (346) to bound it: S 2 (ψ) ≤ S 1 (ψ) ≤ log rank(ψ) (347) therefore the extractable work is: Now we can take the Isospectral twirling of this quantity U → U G = G † U G, defined in Eq. (2): Now taking a look at the result for the entanglement Eq. (43) and following a procedure identical to the one used for the work H.1.2, we find: Recall Eq. (262), Eq. (263) and Eq. (72) and consider the difference between the extractable work for t = 0 and in t: where we defined δW 0 ≡ E 0 − tr (H 0 )/d, that is the difference between the expectations value of the Hamiltonian evaluated on the state ψ and on the infinite temperature state 1l/d. The bounds for the extractable work in the large d limit can be written as: where we introduced a dimensionless parameter ≡ ∆W 0 β, and the quantity∆ F G ≡ ( F G (t) − F G (0))/δW 0 . In the large d limit we have:

I Isospectral twirling and CP maps
Recall the definition of Isospectral twirling for CP-maps Q(·) = α K α (·)K † α : Let us first give the expression for k = 1. From Sec. 4.1, the computation follows straightforwardly: Now let us prove the expression Eq. (88). Starting from the definition: L 1 (Q) = tr (ψQ(ψ)) = α tr (ψK α ψK † α ) and using the trick of the swap operator T ∈ B(H ⊗2 ) we can write it as: then, twirling the Kraus operator as K α → G † K α G, taking the Haar average, recalling the definition K = α K α ⊗ K † α : Now, let us prove the expression for the purity, Eq. (91). Starting from the definition: setting ψ ∈ pure states, use the property T ψ ⊗2 = ψ ⊗2 and T 2 = 1l to write the above expression as: recall that tr ((A ⊗ B)(C ⊗ D)) = tr (T (13)(24) A ⊗ B ⊗ C ⊗ D) and write: Pur(Qψ) = αβ tr (T (13) (24) in order to have the correct expression in the definition Eq. (84), we need to manipulate it: act adjointly with T (234) inserting T (234) T (243) = 1l twice and get: taking the Isospectral twirling we finally have: