Quantum Chaos and Level Dynamics

We review the application of level dynamics to spectra of quantally chaotic systems. We show that the statistical mechanics approach gives us predictions about level statistics intermediate between integrable and chaotic dynamics. Then we discuss in detail different statistical measures involving level dynamics, such as level avoided-crossing distributions, level slope distributions, or level curvature distributions. We show both the aspects of universality in these distributions and their limitations. We concentrate in some detail on measures imported from the quantum information approach such as the fidelity susceptibility, and more generally, geometric tensor matrix elements. The possible open problems are suggested.


Introduction
It is a great pleasure to be able to contribute to the volume of Giulio Casati's 80th birthday. Since the very beginning of my encounter with quantum chaos, Giulio Casati was one of those whose works inspired younger people. As an example, let me mention a contribution of late Prot Pakoński with whom I had a pleasure to consider the Kepler map (one of the toy models of Giulio), extending it to arbitrary polarization of the microwaves [1]. In this review I will discuss, however, a different topic-the statistical measures related to level dynamics in quantally chaotic systems. This is where we met scientifically, contributing to a single paper I had an honour to coauthor with Giulio [2]. The subject of level dynamics developed mainly in the 1980s and 1990s of the former millennium, yet it has recently found interesting extensions and applications in the modern many-body physics.
Level dynamics, described in a pedagogical way in the books of late Fritz Haake [3] and Hans-Jürgen Stöckmann [4], considers the motion of levels as a function of some arbitrary scalar parameter λ which characterizes the Hamiltonian H(λ) of the system. It may be viewed as the motion of interacting fictitious particles (represented by levels) with λ being the effective time, as described originally by Pechukas [5] and followed by Yukawa [6] who built a corresponding statistical mechanics picture. This was applied in different ways to either justify the random matrix theory application to quantum chaotic spectra (this aim has been, in fact, accomplished by a beautiful periodic orbit theory application by the late Fritz Haake and collaborators-for a review see the Ref. [3]) or to define new statistical measures of energy spectra and find their distributions. Without any claim for completeness, some of these applications will be reviewed below.
We begin with defining the notation for level dynamics and the corresponding statistical mechanics picture in Section 2 showing how standard results from this approach provide a prediction for level statistics in a chaotic integrable transition. The resulting distribution seems to work surprisingly well for the data in many-body localization crossovers [7]. To somehow complete the picture, we review in the next Section several other models for the statistics in the transition, notably the banded matrix model, developed by Casati and coworkers [8][9][10][11][12]. In Section 4 we review the universality conjecture in level dynamics [13] showing how it is reflected in the so-called curvature distributions [14] in Section 5. We mention the velocity correlations in Section 6 and the avoided crossings statistics [15,16] in Section 7. Applications of these measures are discussed, stressing their limitations in real systems. Section 8 describes, on the other hand, very recent findings on the distribution of fidelity susceptibility [17], while the extensions to many parameter dynamics with geometric tensor matrix elements distributions are reviewed in Section 9. We mention briefly the parametric measures in the transition to localized regimes in Section 10. We finish with conclusions discussing future perspectives.
While the presentation that follows is only theoretical, we should mention that the predictions concerning level dynamics were tested, to some extent, in experiments. As typical for quantum chaos, those experiments were not carried out on eigenvalues of the Schrödinger equation, but rather on related models of quasi-2D microwave cavities or propagation of acoustic waves. We provide an incomplete list of references to beautiful experiments [18][19][20][21][22][23][24] stressing the work of the Stöckmann group [18] where a rather complete comparison of different measures with experimental microwave resonance data was carried out.

Level Dynamics
Let us recall some basic details on level dynamics, to fix the notation. Let the Hamiltonian, H(λ) = H 0 + λV, depend on some parameter λ for arbitrary H 0 and V. The eigenvalue equation (where E(λ) a is the eigenvalue corresponding to eigenvector |a(λ) ) upon differentiation with respect to λ gives Let us call p a ≡Ė a a velocity of level E a where λ becomes a fictitious time. Such a "motion" for a set of levels {E a (λ)} is visualized in Figure 1. A λ derivative (denoted by a dot over the variable) of p a yieldsṗ where we introduced f ab = V ab (E a − E b ). Following the procedure one step further, we find equations forḟ ab :ḟ One notices that no new quantities appear, and the set of Equations (2)-(4) is closed. It is sometimes called the Pechukas-Yukawa equations following original works [5,6]. This set of nonlinear equations is integrable-this is not surprising as a problem is equivalent to a diagonalization of the Hamiltonian matrix as noted in the Ref. [3]. Each level may be visualized as a particle with the position given by the energy, E a (λ), the velocity given by the level slope dE a /dλ, and acceleration (curvature of the level, d 2 E a /dλ 2 . The eigenvalues of H = H 0 + λV are V-dominated for large λ and the dynamics become trivial. Haake [3] introduces a different λ dependence (equivalent for small λ): we shall follow the "trigonometric choice" [4,14,16] H = H 0 cos(λ) + V sin(λ).
This results in an additional harmonic binding of eigenvalues which prevents them from escaping to infinity. In effect, the equations of motion become: with Since the dynamics are integrable, the appropriate statistical description should involve all possible constants of the motion. Such an approach would be a formidable task. The Yukawa simplified way is just to consider the simplest integrals of motion, the total energy H and the trace of the square of f ab matrix, Q = 1 2 Tr(F 2 ) [4] with The phase space density, according to Gibbs, is: The reader may be surprised that we use α as an effective inverse temperature. The simple reason is that, in order to follow the sacred quantum chaos notation, we reserve β for a level repulsion parameter with β = 1, 2, 4 characterizing different universality classes of Dyson and corresponding, for Gaussian ensembles to Gaussian Orthogonal, Unitary, and Symplectic Ensembles (denoted as GOE, GUE, and GSE, respectively). The density ρ may be explicitly written out as (10) By integrating out the variables p n and f nm , we can compute the joint probability distribution (JPD) of eigenvalues [4] P(E 1 , with β = 1, 2, 4 corresponding to three universality classes. The β appears in (11) as the structure of the integrated variables f nm depends on the universality class with F being orthogonal, unitary or symplectic. The similar ensemble was considered by Gaud in the Ref. [25], as well as Forrester [26] and Hasegawa and Ma [27]. They considered mainly two point correlation functions for the unitary ensemble. We rather concentrate on the time-reversal invariant case, as most commonly met in current many-body localization studies. Equation (11) is simplified in different limiting cases. The Poissonian distribution was reached in the γ/α >> 1 limit. The distribution becomes On the other hand, for γ/α 1 the distribution yields the Gaussian ensemble result. We have Finally, to reach (GOE) in this limit, we fix β = 1, and we also fix, for simplicity, α = 1 (the latter choice affects the global energy scale only). We also denote γ/α = 10 p . The distribution (11) takes the final form where p = log 10 γ β is the single parameter interpolating between GOE (p → −∞) and Poisson (p → ∞) limits. The first term in (14) represents the pairwise interaction between the particles and the exponential term provides the harmonic binding of the eigenvalues. The resulting distribution, obtained using Monte-Carlo sampling for different p, was shown to faithfully reproduce statistics of eigenvalues on the transition between ergodic and many body localized situations [7].

Other Interpolating Ensembles
It is interesting to review several interpolating statistics models proposed in the past. An early work of Rosenzweig and Porter [28] is certainly worth mentioning. In this model the variance of off-diagonal elements in random matrices is controlled by a matrix dimension-dependent parameter. Its value interpolates between the Gaussian orthogonal ensemble (GOE) value to vanishing values for the Poissonian case. The other approach was proposed on the basis of a Wigner-inspired 2 × 2 matrix approach by Lenz and Haake [29]. Yet another was an ad hoc expression known as Brody distribution [30] which fit to low-resolution experimental data surprisingly well. Berry and Robnik [31] proposed the distribution based on sound physical assumption of the separation between "chaotic" wave functions faithful to GOE and those localized in the regular part of the phase space. The corresponding distribution was shown to work well in the so-called deep semiclassical limit [32]. Another proposition due to Seligman and coworkers [33] assumed that the variance of off-diagonal elements a ij should scale as exp −(i − j) 2 /σ 2 . For σ → 0 one recovers the Poisson case while σ → ∞ becomes GOE. Yet another well-known approach is that of Guhr [34] who used supersymmetric techniques to express the two-level correlation function in the Poisson-GOE ensemble in terms of a double integral. It is worth stressing that another popular proposition was advocated by Giulio Casati and coworkers [8,9]. They considered banded Gaussian random matrices as a useful tool in describing the transition, where the corresponding parameter was y = b 2 /N with b being the matrix bandwidth and N its rank. While the (unfolded) level spacing statistics was the main object of quantum chaos studies, in a many-body localization (MBL) context Huse and Oganesyan [35] introduced a new measure, the gap ratio, defined as r n = min[δ n , δ n−1 ]/max[δ n , δ n−1 ], where δ n = E n+1 − E n is the energy gap between the consecutive energy levels. The dimensionless gap ratio does not require unfolding. The MBL transition description was addressed by Serbyn and Moore [36] who proposed two stages of the GOE-Poisson transition: (1) A Short Range Plasma Model (SRPM) and (2) a semi-Poissonian level statistics [37,38]. Recent efforts worth mentioning are a β-Gaussian (β − G) model [39]. A comparison of the performance of different models was given in the Ref. [40], while the Ref. [41] proposes a more complicated, two-parameter β − h model, where the pairwise interaction between the levels is limited to a number h.
Comparison of some of these models with numerics for interacting disordered spin systems modelling ergodic to the MBL transition is given in the Ref. [7]. The interested reader should consult the Ref. [7] for details, but it suffices to say here that the single-parameter Yukawa-like model described above compared favorably with other single-parameter models and quite faithfully reproduced the disordered spin data for the MBL-ergodic crossover.

Universality of Parametric Dynamics
A simple inspection of Equation (10) shows that the velocities, p n , have, in this approach, a Gaussian distribution with the variance determined by the "inverse temperature" α. This is the essence of level dynamics universality, as was determined and thoroughly studied by Simons and Altschuler [13,42]. The level spacings have a single scale-the mean level spacing, ∆. The unfolding then corresponds to rescaling the energy levels i = E i /∆. Level dynamics introduces a novel scale determining how fast the eigenvalues change with the parameter λ. The original definition [13] involves the velocity-velocity correlation function for unfolded levels averaged over eigenstates n. Then C(0) yields the second, apart from the mean level spacing, important scale. When the levels are unfolded using mean spacing, ∆, and the parametric dependence is "unfolded" using C(0) as [13] x = C(0)λ the spectral properties of different systems should be similar. Clearly, C(0) in our notation is directly related, modulo-unfolding to the "inverse temperature" α in the Gibbs ensemble.
There is a hidden assumption in the level dynamics universality as formulated by Simons and Altschuler [13,42] and apparent in the Pechukas-Yukava statistical approach-the parameter change is global in a sense that both H 0 and V are of similar "strength" i.e., belong to the same ensemble. This does not have to be so, as nicely exemplified in the Ref. [43] where it is shown that small local perturbation modifies strongly the velocity distribution in the billiard example studied. Assuming random wavefunctions it is derived that in such a case P(v) ∝ K 0 (A|v|) where A is a constant and K 0 (.) denotes a modified Bessel function. The distinction between local and global perturbations in the context of unviersality were further studied in the Ref. [44].
It seems natural to review now the properties of velocity correlation function. For the reasons that become obvious later, it is more convenient to consider first the second derivatives of energy levels with respect to the parameter, the so-called curvatures.

Curvature Distributions
The curvatures of levels K n =ṗ n as derivatives of velocities should be called in the dynamics language "level accelerations". We stick to curvatures as this is a commonly used name. Large curvatures appear in the vicinities of avoided crossings in the system. Then, essentially only two levels are involved. Following this strategy [45] showed that the large curvature tail behaves as |K| −(β+2) for all three universality classes.
The full analysis of curvature distributions, not limited to large curvature tail, was carried out in the Ref. [14]. Large numerical data collected for all three ensembles suggested the following simple and analytic form: (with β = 1, 2, 4 for GOE, GUE and GSE, respectively) and where, recall, ∆ is the mean level spacing (i.e., an inverse of the mean density of states). Defining the dimensionless curvature, k, as we have explicitly These expressions, which could be claimed as being determined via Monte-Carlo integration and inspired guess, were soon proven analytically for all three ensembles of Gaussian random matrices [46,47] and by an alternative technique in the Refs. [48,49]. Let us remark that the above definition of k differs from the form suggested by the universality rule, (16), d 2 /dx 2 by a multiplicative factor πβ which simplifies (20)- (22).
The distributions (20)- (22) appear to work well for circular ensembles as well as some quantally chaotic systems such as kicked tops [14] or periodic band random matrices in the metalic regime, as shown by Casati and coworkers [11]. The question remains, however, to what extent these RMT predictions are universal and to what extend the particular quantally chaotic systems are faithful to them. The first aspect was clarified by Li and Robnik [50] who pointed out that a nonlinear transformation from λ to some other parameter µ(λ) leads to a different curvature distribution as the transformation between curvatures is nonlinear. It reads [50]: In the expressions above k µ and k λ are normalized curvatures calculated with respect to the corresponding parameters, p λ -the slope and prime denotes derivative with respect to λ. As Li and Robnik [50] point out since velocities are Gaussian distributed (fast decaying) the universality of curvatures may be restored for large curvatures but the nonlinearity of the transformation (23) prevents universality at all scales, see also [51]. The same argument shows, however, that for any "local" linear transformation the universality may hold. As long as changes of H(λ) are linear in λ, as assumed in the derivation above, one might expect the universality to hold. There is, however, another origin for the lack of universality which gives us insight into the physics involved. Already Takami and Hasegawa [52] suggested that the presence of scarring, i.e., strong localization of eigenstates in the space where unstable periodic orbits exist in the classical limit [53,54] may affect curvatures. Numerical studies of several examples such as the hydrogen atom in a magnetic field or quantum billiard proved that this is indeed the case [14]. While referring the reader to an original paper for numerical details it suffices to say that strong scarring leads to a peculiar level dynamics with some levels (scarred eigenstates) have quite different slope than the rest and interact with other levels only locally in narrow avoided crossings. Those levels behave like solitons and may be described as such [55,56]. Their behavior leads to an excess of small curvatures (outside of these avoided crossings) as well as very large curvatures (at the centers of avoided crossings).
It is worth stressing (which we just do with a single sentence) that the curvatures are strongly linked with transport and conductance [11,57]. Particularly interesting for this case are situations where the parameter breaks time reversal invariance as it happens for twisted boundary conditions.

Velocity Correlations
Let us come back to the level slopes, i.e., velocity correlations. Simons and Altschuler [13,42] in their analysis considered the autocorrelator of velocities at a some energy difference, ω, c(x, ω) (note-this is a different object than C(x) which involves correlations for the same level n) that involves all level velocities in a given interval studied. We refer the reader to original papers for details [13,42]. The C(λ) or rather C(x) was studied numerically [58] for all three unitary classes. A simple analytic approximation for C(x) was proposed in terms of the plasma error function, see [58].
Application of Taylor series expansion of C(x) allows one to link the velocity correlator to the variance of the rescales curvatures. Explicitly, one obtains [58] Defining normalized correlation c(x) = C(x)/C(0) one reproduces the results [42] for GUE: c GUE (x) = 1 − 2π 2 x 2 and gets c GSE (x) = 1 − 8 3 π 2 x 2 + . . . for the symplectic ensemble [58]. Interestingly, for the most common orthogonal universality class one encounters the problem as the variance of curvatures, following (20), does not exist. This indicates that the small x behavior may be singular and the Taylor expansion is questionable. This issue has been studied further in the Ref. [2] where it was shown that in fact c(x) reveals singularities around x = 0. Taking the parametric dependence (5) one may show that Fourier components of the Fourier transform of c(x) have algebraic tails which directly indicates singularities at x = 0 of the velocity correlator. Again we just quote the the final result which shows that with b i being coefficients of the order of unity. One may observe that the singularity at x = 0 becomes weaker with growing level repulsion β, being most severe for GOE. An even more in depth analysis of singularities appears in the Ref. [59] where explicit values for the parameters, b i are found. As mentioned in Section 4 local (instead of global) perturbation affect strongly velocity distributions [43,44]. This has also a pronounced effect on velocity correlations [43,44] as further analysed in detail for the unitary ensemble [60].

Avoided Crossings Distributions
Another statistical property with interesting links to level spacings is the distribution of avoided crossings, i.e., minima of distances between neighboring levels. The problem of finding the corresponding distribution was formulated by Wilkinson [61] who has shown that integrated distributions for small minimal distances D for GOE (GUE) show similar repulsion as present in spacing distributions. Avoided crossings for billiard models were numerically studied by Goldberg and Schweizer [62]. While in a general case the exact distributions are not known in some analytic form, a well working approximations based on two-levels approximation may be easily derived following the Wigner approach for level spacings themselves [15]. For GOE case it is written down immediately as the two-level Hamiltonian H = H 0 + λV may be written in the eigenbasis of V as The minimal distance between levels is simply 2|d|. Since H 0 is assumed to correspond to GOE, d is Gaussian distributed, so we get the distribution (for D = 2|d|) normalized to unit mean avoided crossing: Situation is only slightly more complicated for other ensembles. For GUE d in (26) should be complex, d = d 1 + id 2 , with independently Gaussian distributed (with the same variance) d i . A simple integral leads to a normalized distribution which is identical to the so-called Wigner surmise for spacings for GOE. We observe a simple rule that the avoided crossings in two-level approximation share the same distribution as the nearest neighbor spacings but the the repulsion parameter β reduced by unity. So for GUE with β = 2 we get the Wigner formula corresponding to spacings for β = 1. This is in full agreement with small D perturbative prediction of the Ref. [61]. Similarly, an explicit calculus shows that for β = 4 GSE the avoided crossing distribution behaves as D 3 for small D.
Numerical tests (which have to be carefully done to correctly estimate and avoided crossing values [16]) show excellent agreement between two-level approximate formulae and numerical data for random matrices of larger sizes. The agreement is in fact better than the spacings and the Wigner surmise. The reason is simple-the two-level approximation works better for minimal distances between levels.

Fidelity Susceptibility
Rapidly developing in the past 20 years, the quantum information field brought yet another measure which may be related to level dynamics, the fidelity, F [63]. While generally defined for mixed states, a pure-state definition [64] suffices for our purposes Here, ψ(λ) is an eigenstate at the value of the parameter equal to λ. Taylor expansion for small λ leads to the definition of fidelity susceptibility, χ Fidelity susceptibility became an indicator of quantum phase transitions. At the transition point, the ground state properties change, leading to the enhancement of χ [64][65][66]. Apart from ground states, thermal states were also considered [67][68][69]). The fidelity statistics, taking into account a set of eigenstates, was discussed for the first time in the Ref. [17]; note that an attempt to identify many-body localization transition was due to the Ref. [70]. We shall briefly review the results of the Ref. [17] that provide one of the rare situations when exact analytic results are available for the arbitrary size of random matrices. Consider H = H 0 + λV with both H 0 , V belonging to GOE or GUE. Fidelity susceptibility of state |n with energy E n of H 0 is easily derived as showing some similarity to curvaturesṗ n (3)-the difference is just a power in the denominator.
The probability distribution of the fidelity susceptibility at energy E reads: which we consider at the center of the spectrum E = 0. Following the technique developed in the Refs. [48,71,72], one arrives [17] at, for the GOE case, where C O N is a normalization constant and Equation (33) is an exact result for an arbitrary rank N of the random matrix from GOE. This is one of the rare situations when analytic formulae for arbitrary N and not only for the N = 2 or N → ∞ limit are known. The N → ∞ limit is interesting. As I O N is asymptotically proportional to N, one can define a scaled fidelity susceptibility y = χ/N. Its distribution, in the N → ∞ limit takes a rather simple form As tested numerically, this expression works well for N ∼ 200 already.
Similarly, see the Ref. [17] for the derivation, where one obtains an analytic, valid for arbitrary N results, for GUE. We quote here just the N → ∞ limit for the scaled fidelity susceptibility while for a full expression for χ valid for arbitrary N as well as for a comparison with numerical data, we refer to the Ref. [17]. The distributions (33) and (34) for GOE and GUE, respectively, are presented in Figure 2. For large scaled fidelities y a power law decay is observed 1/y 2 for GOE and 1/y 5/2 for GUE.

Generalization to More Parameters
A natural extension of parametric level dynamics occurs in the presence of more than one parameter. One may define where H 0 , V i are statistically independent and drawn from the same (as in this review) or different universality classes. Obviously, novel problems appear in that case which we describe briefly only. Probably it was Michael Wilkinson and his student, Elisabeth Austin, who addressed first such a situation in their study of density of degeneracies, for example, diabolical points [73]. A three-parameter family was considered for Chern integer fluctuations [59,74]. Steuwer and Simons [75] found the distribution of adiabatic curvature (related to Berry phase) for GUE. The multiparameter dynamics was recently revisited by Berry and Shukla who discussed Berry curvature deriving two-level and three-level distributions [76][77][78]. They found the large curvature scaling, P(c) ∼ c −2 for the orthogonal and P(c) ∼ c −5/2 for the unitary class, the result already present in the Ref. [75]. Importantly, however, Berry and Shukla [78] linked the problem with the quantum information concept of geometric tensor and the distance between quantum states [79][80][81]. The Fubini-Study distance in the Hilbert space between states differing by a small change of parameters form λ ≡ (λ 1 , λ 2 , . . . , λ n ) to λ + d λ is where g (n) ij ( λ) is the so-called geometric tensor [79][80][81] which governs the quenches (in λ) of the system. For (36) Note that the distance between states is determined by the real part of the geometric tensor only- (37). The imaginary part, Img ij , gives the Berry curvature [78,82] related to changes of two parameters λ i , λ j . For a single-parameter problem the geometric tensor reduces to a scalar proportional to the fidelity susceptibility discussed in the previous Section. Additionally, if V i belongs to the same ensembles, then the distributions of g ii are equal. One may consider, therefore, the trace G = Trg (n) ij as an equivalent of the fidelity susceptibility while the distribution of the imaginary part of the geometric tensor reduces to the Berry curvature distribution.
The distribution of trace, G, is, therefore, given by (33) discussed above, valid for arbitrary N [17] after taking into account the fact that, for equal diagonal elements the trace is simply the matrix rank multiplying the diagonal element. The alternative derivation using supersymmetric techniques is provided in the Ref. [82] for GUE in the N → ∞ limit. The same authors obtained for the Berry curvature the result derived earlier by the Ref. [75].

Towards the Localization Limit
While in the introductory part we have considered the transition between an integrable (localized) regime in the context of level spacings and the dynamics of Pechukas-Yukawa gas, later we mainly described the results for the ergodic regime well-described by Gaussian random ensembles. Here we shall briefly mention some of the results for level dynamics measures that involve the transition.
Here the seminal contributions were provided by the Como group centered around Giulio Casati [8][9][10][11]. The banded random matrix ensemble provided a natural tool to study the transition from ergodic (metallic) to localized transition by varying the width of the band. The team addressed also curvature distributions [2,11]. The analytic approach to the problem was pursued by Yan Fyodorov, who, starting around 1994, considered comprehensively level dynamics features close to the localization transition studying velocity correlations [83,84] or curvature distributions [85]. In particular, the velocity distribution for one-dimensional disordered wire is derived using the supersymmetric approach to be for the scale velocity v s . The curvature distributions were more recently addressed in the context of MBL studies [86,87]. The level dynamics across the many-body localization transition for the paradigmatic XXZ spin model was considered in the Ref. [88]. Velocity, curvature and fidelity susceptibility distributions were considered. Interestingly, while velocities depended on the choice of the parameter (being the interaction strength or the kinetic tunneling), curvature exhibited universal behavior in the delocalized regime. In the localized regime, curvature distributions reveal system specific characteristics that exemplify the presence of local integrals of motion in the localized phase. Large curvature or large fidelity susceptibility tails change their slope when entering the localized regime. Such behavior is well-understood qualitatively and linked to weaker level repulsion. One may also mention here a proposition to use adiabatic eigenstate deformation as a probe of entering the localized regime [89].
Recent studies address avoided crossing distributions and level dynamics on the localized side of the MBL transition [90][91][92]. Such an analysis, however, has to face the problem of assumed integrability which leads necessarily to exact crossing of levels associated with different quantum numbers.

Conclusions-Where Do We Stand
With this traipse through the last 40 years, starting with the Pechukas model [5], we hope to have shown that many fascinating results have been obtained in the studies of level dynamics of complex systems, but that there are still many open questions and unsolved problems. There are at least two areas that await a more decisive attack and, hopefully, solutions. One is the transition between different ensembles [93]. We have not, on purpose, reviewed only a few studies in this domain for the reader to formulate their own problems. A second related area with several white spots lies in multiparameter level dynamics. The latter has been mostly limited to studies within unitary ensembles, notoriously easiest to treat. We do not know the Berry curvature distribution for the orthogonal ensemble. We do not know about the geometric tensor properties when different parameters induce different transitions. Even the simplest questions remain unanswered. For example, simple analysis shows that fidelity susceptibility of Berry curvature decays with reverse quadratic power for GOE while the corresponding power is −5/2 for GUE. Can we say something when we couple GOE and GUE ensembles by some parameter, such as weakly breaking the time reversal invariance? Can we generalize the findings to the symplectic ensemble? What about the ten-fold method [94]? It is my belief that we may expect in the future some very interesting results coming from new people entering the subject. I already anticipate the excitement.