Frame potential of Brownian SYK model of Majorana and Dirac fermions

We consider the Brownian SYK, i.e. a system of $N$ Majorana (Dirac) fermions with a white-noise $q$-body interaction term. We focus on the dynamics of the Frame potentials, a measure of the scrambling and chaos, given by the moments of the overlap between two independent realisations of the model. By means of a Keldysh path-integral formalism, we compute its early and late-time value. We show that, for $q>2$, the late time path integral saddle point correctly reproduces the saturation to the value of the Haar frame potential. On the contrary, for $q=2$, the model is quadratic and consistently we observe saturation to the Haar value in the restricted space of Gaussian states (Gaussian Haar). The latter is characterised by larger system size corrections that we correctly capture by counting the Goldstone modes of the Keldysh saddle point. Finally, in the case of Dirac fermions, we highlight and resolve the role of the global $U(1)$ symmetry.


Introduction and main results
The phenomenon of information scrambling in quantum many-body systems has recently gained significant attention due to its deep connections with holography, black hole physics, and quantum information [1][2][3][4][5][6][7][8].While some quantum scrambler quantifiers are inspired by the classical theory of chaos, as is the case for the OTOCs (out-of-time-order correlations), a slightly different way to look at information scrambling is via k-design realisations.Namely, one asks when and how the unitary time evolution from a pure initial state can generate an ensemble of random quantum states (under a given definition of such an ensemble) which, up to the k-th moment, is uniformly (Haar) distributed [9][10][11][12].Such a question is directly related to the emergence of eigenstate thermalization (ETH) in closed quantum systems [13][14][15][16][17].More recently, a related problem has been posed in the context of deterministic Hamiltonian systems undergoing final read-out measurements: the randomness induced by quantum measurements can be approximated with different accuracy by a Haar random ensembles, a feature recently dubbed deep thermalization [11,[18][19][20].Generic chaotic quantum systems are indeed expected to thermalise at a late time, i.e. the reduced density matrices obtained by tracing out most of the degrees of freedom are expected to converge to the standard Gibbs statistical density matrix, which reduces to the identity matrix in the absence of conserved quantities.However, the presence of different sources of stochasticity (e.g.noise or measurements), raises the question of whether, beyond the convergence of the stochastic average, individual realisation (or trajectories) of the system are exploring uniformly the accessible Hilbert space of pure states.However, the study of strongly interacting systems, in particular with the additional ingredient of stochasticity, is generally analytically (and also numerically) very challenging, except for a few specific systems made of random matrices or dual unitary quantum circuits [21][22][23].For this reason, it is important to develop an intuition on some sufficiently generic but tractable cases.Among these, a prominent example is the Brownian SYK model [24,25], describing L Majorana (or Dirac) fermions with all-to-all interactions (as for the standard SYK), with a time-dependent white-noise interaction coupling.The model is well known to serve as a paradigmatic model for quantum scrambling and chaos [7,[25][26][27][28][29].
In this paper, we focus on investigating the degree of mixing in the time evolution of the Brownian SYK model.To address this question, we employ a measure of mixing called the Frame potential.Specifically, we consider two identical copies of an initial density matrix ρ 0 and compute the overlap between the time evolutions of these copies under two independent realisations of the noise.The Frame potential, denoted as F (k) (T ), is given by the expression: Here, U 1 and U 2 represent the (noisy) evolution operators of two independent Brownian SYK systems with randomly chosen couplings and F (k) equals the k-th moment (averaged over both noise realizations) of their overlap.In the case of sufficiently mixing time evolutions and starting from a reference pure state ρ 0 = |Ψ⟩ ⟨Ψ|, the Frame potential exhibits an exponential decay from its initial value one.This decay manifests the decreasing probability that the two density matrices remain close as time progresses.Eventually, at infinite time T → ∞, it is expected that the density matrix ρ(T ) fully explores the entire Hilbert space in a uniform manner, modulo the constraints imposed by global symmetries as U (1), implying that its statistics can be effectively obtained by replacing the time-evolution operator with a Haar distributed unitary operator in the Hilbert space H as ρ(T ) → U † ρ 0 U, U ∼ Haar(H). (1.2) Consequently, the Frame potentials for all k's are expected to converge at large times to their Haar-averaged value, i.e.where both U 1 (T ) and U 2 (T ) are replaced by two independently drawn Haar-distributed unitary matrices.We denote this limit case as F

(k) Haar
and it is easy to prove (see for instance [23]) that it corresponds to the minimum overall The rescaled Frame potential with q = 2 in the logarithmic scale, exact numerical simulations compared with the theoretical predictions.The number of fermions in both cases is L = 120.Brown colour corresponds to the number of replicas k = 3, green colour corresponds to k = 2.The dotted line is the theoretical prediction of the Keldysh calculation at early times, and the dashed line is the prediction for the late time behaviour.From the numerical data, we can see that the time of the saturation can be estimated directly from the theoretical prediction, modulo finite size corrections (which are larger for larger k as expected).The intersection between the short-time saddle point (dotted line) and the long-time saddle point (dashed line) gives the estimation for the saturation time.In the case q = 2: as compared to the case q > 2, where we have possible distributions of unitary operators on H.As detailed in this paper, this convergence indeed occurs for any q > 2 (with interactions coupling more than 2 fermions), thus generalising the results of [8].An exception arises in the Gaussian integrable case q = 2, with exactly two-fermions interaction [30][31][32].In this case, the evolution only exhibits mixing behaviour within the space of Gaussian states.There, we can introduce a Gaussian Haar (gHaar) [33] measure and demonstrate that the Frame potentials converge for all k to the latter at late times.Intriguingly, the value of the gHaar-integrated Frame potential is equivalent to the generic case in the limit of a large fermion number L, with logarithmic corrections in L. As demonstrated in this work, these corrections arise from fluctuations of Goldstone modes around the late-time saddle point, which are only present in the Gaussian case where both the Keldysh action and the saddle point's symmetries are generated by continuous groups.Moreover, we show that in the case of Dirac SYK the global U (1) symmetry leads to a late-time convergence (1.2) in each sector of a given fermionic charge, giving a radically different behaviour between the non-Gaussian and the Gaussian case.In both cases, the Frame potential converges at late times to a larger value compared to its Majorana equivalent, as the result of the global constraints given by the U (1) symmetry.
Summarizing our results here, we find the following: • the Frame potential decays at short times as Haar = e Dirac.
• For q > 2, given the dimension of the Hilbert space N , the Frame potential saturates at large times to the values Dirac.
• For q = 2, including the powers in L due to massless fluctuation around the saddle point, the Frame potential saturates at large times to the values, whose leading terms in L are given by The structure of this paper is as follows.In Section 2, we derive the Frame potential for the Brownian Majorana SYK system using the framework of the thermofield double state.We outline the Keldysh path integral approach, which facilitates averaging over the Brownian noise, and obtain an effective action for the problem.Moving to Section 3, we explore the Swinger-Dyson equations by performing a variation of the action.Initially, we consider a naive saddle point solution that does not mix different Keldysh contours and replicas.The resulting Frame potential exhibits exponential decay, as depicted by the dotted lines in Fig. (1).Subsequently, we investigate a non-trivial solution in the large time limit, wherein different replicas and contours are mixed.Notably, this solution exhibits distinct group symmetries for cases with q > 2 and q = 2.The continuous symmetry group characterizes the q = 2 case, while discrete symmetries are present for q > 2. We then proceed to analyse the fluctuations around the saddle points.For the q > 2 case, these fluctuations trivially cancel with the normalization.However, in the q = 2 case, fluctuations governed by the Goldstone theorem lead to massless modes in the action.Consequently, the integration over these fluctuations deviates from the normalization, resulting in polynomial corrections (in system size L) to the Frame potential.During the evolution, the latter reaches a saturation value exponentially small with respect to the system size.This behaviour is illustrated by the dashed lines in Fig. (1).In Section 4, we extend our analysis to the complex, Dirac, Brownian SYK model.Here, we show how to properly take into account the global U (1) symmetry in the late-time Keldysh saddle point solution.Finally, in Section 5, we conclude the paper by summarizing our findings.We also discuss potential future directions and prospects for further exploration within this research framework.

Brownian SYK model and Keldysh path integral representation of the Frame potential
We consider the time evolution generated by the time-dependent Hamiltonian where χi are L Majorana fermions and anticommutation relations are applied { χi , χj } = δ ij , here q is an even integer constant.This model has close analogies with the celebrated SYK [25,34,35], but with the important difference that h i 1 •••iq are not constant in time here.On the contrary, we take them to follow a white noise distribution, where we denote the collective set of indices i = (i 1 , . . ., i q ) (similarly for j) and set δ i,j = δ i 1 ,j 1 . . .δ iq,jq .For this reason, the current model is named Brownian SYK.We shall here assume L to be an even integer, such that the operators χi can be embedded into a N = 2 L/2 dimensional Hilbert space.We are interested in studying the scrambling dynamics induced by the time evolution (2.1).In order to do so, we focus on a specific initial condition, known as thermofield double (|TFD⟩) [6,36,37], which has already been employed in the context of Brownian SYK in [26].In practice, we consider two copies of the system, prepared in a maximally entangled state.In the following, we shall address the two copies as left (L) and right (R) halves.The selection of the |TFD⟩ state is not uniquely determined and relies on the choice of a basis.However, different definitions yield states that are connected through local unitary transformations.To provide a concrete illustration, we offer an explicit characterization of |TFD⟩ in the context of the SYK model.Within this framework, the Majorana fermion operators are denoted as χj,L and χj,R , where j takes values from 1 to L. (It should be noted that L must be an even number to ensure a welldefined Hilbert space).A convenient specification for the state can be formulated using Dirac fermions ĉj,L − ĉ † j,R |TFD⟩ = 0 and ĉ † j,L − ĉj,R |TFD⟩ = 0, ( which can be expressed in terms of Majorana fermions as Then we have the standard definition in terms of Majorana fermions: χj,L(R) + i χj,L(R) |TFD⟩ = 0, ∀j = 1, ..., L. (2.5) From the eq.( 2.3) we see that the TFD state in the eigenbasis of n j,L(R) = ĉ † j,L(R) ĉj,L(R) basis can be expressed as a product of |0⟩ L |0⟩ R + |1⟩ L |1⟩ R on each site.More explicitly, we consider the initial state: where each |j⟩ represents a possible string of zeros and ones.Here N is the dimension of the Hilbert space.This state has the following property, where the trace on the left is on the doubled Hilbert space, while the one on the right is on the single copy.This identity holds for arbitrary operators A, B/C, D acting respectively on the first/second copy of the Hilbert space.Except for the initial entanglement, the two halves evolve according to two uncoupled unitary operators, where U h L and U h R are the unitary time evolution operators acting on the two halves.We assume that the time evolution in each half is generated by an independent realisation of the Brownian SYK where the subscript σ = L, R labels the corresponding half and H L (t), H R (t) have the form of eq.(2.1) in terms of two sets of Majoranas χi,σ and independently generated white noises h σ i (t).Starting with this initial state, we calculate the k-th moments of the Frame potential averaged over Hamiltonian realisations with measure dη(h): here, (T ), (and the same for h ), where the couplings h 1(2) are independent random variables with the distribution dη (h) = e − h 2 2σ 2 dh and variance defined in eq.(2.2).The initial condition ρ 0 is the double thermofield state that we discussed earlier.Then, the Frame potential takes the form, for generic k, where the evolution operator U h (T ) is now only acting on a single copy R (or L) and where N is the dimension of its Hilbert space, coming from the definition of the initial state eq.(2.6).In order to get to this expression, we used the property eq.(2.7) together with The way of calculating this integral is the following.First, we move to the path integral representation using Keldysh technique.Then, we shall do an average with respect to Gaussian random variables, and finally, calculate the path integral using the saddle-point approximation.
In order to use Keldysh path integral formalism, we need to construct proper coherent states [38].To do so, we introduce an additional set of Majorana operators η i and create fermionic operators in the form: (2.12) With this notation, the coherent state can be defined as where ψ is a Grassmann variable that parameterizes our state.The set of these states is not orthonormal, and their overlap is given by: Concerning the path integral formalism, we need to get two more useful identities.First, using the eq.( 2.14) we can express an identity operator as an integral over the Grasmann fields: and secondly we can express the trace of an operator Ô in terms of Grassmann variables as: here −1 comes from commuting left and right coherent states.Using eq. ( 2.15) and eq.( 2.16) we can rewrite eq.(2.11) as following where s = ± indicates one of the two Keldysh contours, l counts replicas, and i labels time.
Each matrix element here has the form, Let us notice that for any matrix element of this form, we can simply replace ĉ † with ψ and ĉ with ψ.We do not need normal ordering here, as H in eq. ( 2.1) has one term in which the order is opposite.We can reverse the order at the cost of a minus sign (we stress that i k ̸ = i l everywhere and all fermion operators anti-commute), obtain the action, and once again reverse the order of the corresponding Grassmann variables ψ, ψ at the cost of an extra minus sign.Therefore, we shall write the usual Keldysh action: Here the continuous limit was taken ϵ → 0, N → ∞, and N ϵ → T .We can now move back to Majorana fermions: note that the variables we consider now are not operators but Grassmann variables and i is the index of the Fermionic degrees of freedom.In this limit, actions for χ and ξ fields separate [38], and the partition function factorizes.Therefore, the part that depends on η is merely a spectator and can be integrated out eq. ( 2.11) giving a contribution to the overall normalisation.This can be fixed using the fact that the frame potential must be equal to 1 at time zero.Therefore, with correct normalisation, we write where we have and F (k) (T ) is un-normalised value of frame potential.This action is a part of the standard path integral formalism for the SYK model [25,35,39].Here, we are interested in averaging the Frame potential over the random noisy couplings h.To do so, we calculate the Gaussian integral over each of the variables and use the property given by equation (2.2), which defines the variance σ 2 : where Let us introduce the bi-local fields Ĝ(t 1 , t 2 ), Σ(t 1 , t 2 ) and perform Hubbard-Stratonovich transformation1 : ) here we suppose the summation over the repeating indexes.Inserting this delta function in the expression for the frame potential we obtain the following expression for an effective action: where Notice that in order to form bi-local fields in the expression eq.(2.23), we performed q(q−1) 2 transpositions of the Grassman variables, which gave us the prefactor in the interaction term (−1) together with (−1) q 2 from the Gaussian integration, resulting in (−1) 2 , which is unity due to considering q as an even integer.In order to proceed further, we integrate out the fields χ i , and the action takes the form: where the Pfaffian is taken with respect to indices in time, in Keldysh contours and replica spaces, and the same is applied to the trace.Here, the indexes s and s ′ take two different values, + and −, which indicate one of the Keldysh contours.These contours are two loops, as can be seen from equation (2.17).Notice that after averaging over the noise, the fields from the different contours become mixed, and in the action, we have not only G ++ and G −− , but also G +− and G −+ .Therefore, we can think of these operators as matrices in the Keldysh contour s, s ′ space.Also, let us take into account that the derivative term 1 2 χ∂ t χ has a positive sign for the + contour and a negative sign for the − contour.Therefore, this term can also be expressed as a 2 × 2 matrix with Keldysh contour indexes in the following way: (2.28) Here, we introduced the Pauli matrix σz in the ss ′ space, and we assume a Kronecker delta in the replica space for this term.In the next section, we are going to calculate this functional integral over the fields Ĝ and Σ using the saddle-point approximation at large L.

Majorana SYK model and saddle points
We are in a position to obtain the equations of motion at the saddle points.To do so, we consider the variation of the action with respect to the bi-local fields: with the tensor ĉs,s ′ = c ss ′ .In this expression, we suppressed Keldysh and replica indexes, for simplicity.Setting the variation to zero and summing over all possible solutions, we can find the time dependence of the Frame potential.

Short-time behaviour
Let us now proceed naively by looking for a saddle point that is diagonal in the +− space.
Then the obvious solution is: With this solution, the Frame potential can be expressed as: This functional determinant can be calculated using diagonalisation of the operator σz ∂ t .

4T
, which prevents us from treating matrices as different matrices.Notice that this determinant should also be taken in the replica space, where we have k copies of the system.Taking the limit N → ∞ and following the regularization technique from the book [40] (notice that we have ω = 0), we get where 2 L 2 is the dimension of the Hilbert space for the initial operator U .Notice that the field Ĝ is also diagonal in the replica space.Therefore, G ++ ll ′ = 1 2 δ ll ′ = −G −− ll ′ , and we have e The Frame potential on this solution takes the form We can observe that in this solution, the Frame potential exhibits an exponential decay with respect to time, and all the replicas are factorised.This behaviour implies that the contribution from this part of the Frame potential is relevant during the early times of the system's evolution but becomes increasingly suppressed as time progresses.

Long-time behaviour
In this section, we consider two fundamentally different cases: q > 2 and q = 2.The difference arises from the fact that in the case of q = 2, the action is invariant under the continuous symmetry group SO(k) × SO(k) (i.e. the replica indices in each Keldysh contour can be rotated separately), therefore the fields Σ and Ĝ can be rotated in the replica space by the generators Âi is a skew-symmetric matrix of the size k and by Âσ z we suppose the tensor product.
Whereas for q > 2, the replica indices are invariant only under the discrete permutation group P k .Therefore, in the first case, the saddle point physics will be characterised by the Goldstone theorem, and the integral over the fluctuations around the saddle point solution will provide us with non-trivial corrections in L, which we also aim to study.Another difference between these two examples is that the space of saddle point solutions for q = 2 forms a continuous manifold, while for q > 2, they are simply a discrete set.

3.2.1
The non-Gaussian case q > 2 Here, we begin by considering the non-diagonal solution in the s, s ′ Keldysh contours space.
To do so, let us examine the following ansatz Here, f (|t 12 |) is an arbitrary function that depends on the module of difference between time t 1 and t 2 , which will be found later and Θ(t 12 ) is a Heaviside theta function.The matrix τ is a permutation matrix in the replica space.let us discuss the function sgn(t 12 ) in the ansatz.In the discretised time one can see that the action consists of the terms G i,i−1 , where i labels time.As the solution depends on the difference between consequent times sgn(t i − t i−1 ) = sgn(ϵ × (i − i + 1)) = sgn(ϵ).Therefore, the sgn(t 12 ) is never equal to zero in the action, which as we will see later helps us to regularise the action.The choice of this ansatz is motivated by the fact that we can permute replicas in the action without altering its form.Notice that in the case q = 2 we have continuous SO(k) × SO(k) symmetry of the action instead of the permutation symmetry for general q.This case will be considered later in this work.By using equation (3.1), we obtain: Now, our objective is to determine the function f (t).To achieve this, we perform the Fourier transformation Σ(t) = dω Σ(ω)e −iωt and substitute this ansatz into the first equation of motion (3.1).This procedure gives: Ĝ(ω) = 1 By performing the inverse Fourier transform and comparing it with the ansatz, we arrive at the following result: (3.12) Our objective is to understand the behaviour of the action on this solution.Firstly, let us note that the second term in the action, eq.(2.27), cancels due to the c ss ′ prefactor On the other hand, the third term will give us a non-trivial contribution Finally, let us proceed with the calculation of the first term in the action.In fact, the Pfaffian can be understood as a fermionic (quadratic) thermal partition function with the real inverse temperature equal to T : where A similar expression was calculated in [41], and here we use the analogous method.Let us observe that any arbitrary permutation can be expressed as a direct sum of disjoint cycles, then Tr(e , where the sum goes over the cycles.Also, we notice that each cycle can be represented in the canonical form, therefore we have where n c represents the length of a cycle, and δ nc+1,β = δ 1,β .Here, sgn(β − α) merely affects the permutation on the boundaries.The boundary conditions have been chosen in a way to keep the parity P equal to 1.It can be observed that for any τ , this solution does not respect the parity of the initial state (the initial state is a thermofield double state between, where for each site, the Fermi parity is P = 1).Therefore, the boundary conditions eq. (3.20) arise from this fact.Then, for each cycle, we have: Introducing complex fermions and inserting this into the Hamiltonian we obtain this is the Kitaev chain [42] with periodic/antiperiodic boundary conditions, eq.(3.20) and odd/even length of chain respectively.The diagonalisation of the chain is readerly done, giving Tr(e −T Ĥ(τc) ) ∼ e Consequently we find the following expression for the Pfaffian: and the Frame potential of this solution can be expressed as: We shall now perform the summation of all the possible saddle points, which is the space of permutation τ of size k.Moreover, we notice that we are selecting only permutations of a given parity, fixed by the initial state.Therefore, we find, for large L, namely the same value of the Frame potential for haar-distributed density matrices (with fixed parity).We have then shown here that the late-time solution converges to the Haar distribution.

The Gaussian case q = 2
In the case where q = 2, the action has the continuous symmetry, SO(k) × SO(k) as discussed at the beginning of the section.Therefore, the ansatz that we examine here consists of orthogonal matrices instead of permutations.Using the same procedure as in the previous subsection, one can see that for q = 2, the solution is: where θ is now an orthogonal matrix in the replica space and not a permutation matrix as in the previous case.Again, we are interested in the calculation of the action on this solution.Notice that the second and third terms in the action are the same as in for permutation matrices.Therefore, our objective now is simply to extend the calculation to the first term, namely the Pfaffian.Let us then apply the same procedure to the orthogonal matrices.Any arbitrary orthogonal matrix can be expressed in block diagonal form: , and θ2m+1 = where, m is a positive integer, and Ri is an orthogonal matrix in two-dimensional space.Then, Tr(e Here, we should also take into account that not all orthogonal matrices will preserve the parity of the initial state.Therefore, Ri should be rotation matrices.Also, in the case of an odd size, the last element on the diagonal is positive unity as we want to preserve parity.Now, let us consider the Hamiltonian for the two-dimensional orthogonal sub-matrix: Here ϕ is an angle that parametrizes the rotational matrix, and we can again move to the complex fermions, see eq. (3.22), which gives a simple two-fermions Hamiltonian: after the diagonalisation the trace can be easily calculated.Using the well-known formula for the trace of quadratic density matrices Tr(e c † i Γ ij c j ) = det(1 + e Γ ij ), we obtain: Tr(e −T ĤR ) = det 1 + exp T q Therefore the Frame potential on the particular solution which depends on an orthogonal matrix θ is As in the previous subsection, the saddle point solution with a given rotational matrix does not depend on the latter.Therefore, we should first consider the fluctuations around the saddle point, and sum over all possible solutions given by the space of orthogonal matrices.As we shall see in the coming section, the zero-mass fluctuations are responsible for different finite-size effects between the q > 2 and the q = 2 case.In the q = 2 case, they carry log(L) L corrections, while with q > 2 there are no zero-mass modes that can fluctuate around the saddle point solution.

Symmetries and Goldstone modes
In the next subsection, we shall calculate fluctuations around the saddle point, which we discussed earlier.For the case q = 2, we need to understand how many modes will be massive and carry non-trivial corrections to the Frame potential.To do so, we can use the Goldstone theorem.Let B be a continuous group of a global symmetry of the action and H is a subgroup of B which leaves the solution of the equations of motion unchanged.Then the number of massless modes is equal to R ml = dim(B) − dim(H).Therefore, the number of massive modes is R m = d − R ml where d is the number of degrees of freedom of our system.Let us start with the identification of the degrees of freedom.By the definition: and , where L is the number of fermionic modes, first we notice that the matrix G ll ′ is a skew-symmetric matrix due to the anticommutation relations of the field χ, which contains 2k(2k − 1)/2 degrees of freedom.Let us check the commutation relations for this matrix expanding the commutator and using anticommutation relations which gives the commutation relations of rotation generators in 2k dimensions.Therefore, we see that by the definition our matrix G ll ′ is a generator of so(2k) algebra.In [43], it was shown that S ll ′ = iG ll ′ satisfy orthogonality relation in the large L limit: which for the matrix Ĝ this implies Ĝ ĜT = I. (3.39) Therefore this matrix should be not only skew-symmetric but also orthogonal.The space of these matrices isomorphic to the SO(2k)/U (k) [44], which gives us d = 2k(2k−1) 2 − k 2 = k(k − 1) degrees of freedom.Now let us consider the group of the global symmetry of the action in more detail.Remind that at finite times, the action looks like eq. (2.27).Therefore, after the integration over the fields Ĝ we have One can see that the problematic term here is L log(Pf(σ z ∂ t − Σ)), which as we already noticed has a rotation symmetry in replica space O 1 = e i Â1 σz , O 2 = e i Â2 1 with Âi skewsymmetric matrix of the size k and by Âσ z we suppose the tensor product.Notice that matrices O 0 = e i Â form a SO(k) symmetry group, therefore the space of the symmetries is SO(k) × SO(k).Therefore, since matrices O i commute with σz , we have where θ is an orthogonal matrix k × k.From the form of the solution, we see that the group that allows us to transform one solution into another one is B/H = SO(k), which means that H = SO(k) too.Therefore, giving the number of massive modes, which will allow us to calculate the integral over the relevant fluctuations.The saddle point approximation gives us the following expression for the Frame potential where the on-shell action S 0 was found in previous chapters and ⟨...⟩ θ means the summation other all the solutions.Now let us consider the quadratic fluctuation term S (2) (δG, δΣ): doing variation we can find The first term here belongs to S 0 , the second term vanishes due to the equations of motion, and the third term is the quadratic fluctuation that we are looking for.Therefore, where the field Ĝ is a solution of the saddle point equation, and we suppose the summation over the repeating indexes.Integrating out the fluctuations, we obtain: Here R m is a number of massive modes.We also need to calculate the normalisation term, which is a standard Gaussian integral and d is a number of degrees of freedom.Putting it all together, we obtain for q = 2: where 2 is a number of massless modes, taking the logarithm of the Frame potential we can get: which coincides with the result for Gaussian Haar calculation, see Appendix (A.8).Notice that for the generic case of q interaction, the symmetry of the action and the solution is discrete, therefore there are no Goldstone modes in this case.Due to this fact, the integration over the fluctuation must be carried over all the degrees of freedom which are dim(SO(2k)/U (k)), also, we have the same degrees of freedom at zero and at finite time.Therefore, the normalization will cancel the fluctuation part from the finite time action.Hence, for the generic q > 2 case we have

Dirac SYK model and Kelsdysh saddle points
We shall now consider the Brownian evolution given by the complex SYK Hamiltonian Here, ĉ † and ĉ are Dirac Fermions, and standard anti-commutation relations are applied {ĉ k , ĉ † l } = δ kl .As in the previous chapter, h i 1 ...iq are normally distributed random variables, but this time, they are complex and have the variance: We shall then repeat the analysis of the previous chapter.The model now possesses an extra symmetry compared to the Majorana case, which is the global U (1) charge conservation We shall then show how this impacts the saddle point solutions at late time, both for q > 2 and q = 2.

Keldysh path integral and saddle points
We again choose the TFD state as the initial state, and the object of our interest has the form given in the eq.(2.11), which can be expressed as two Keldysh contours eq.(2.17).The evolution operator is given by U i (ϵ) = e −iϵH(t i ) , where the Hamiltonian is from the eq.(4.1).Now we want to derive an effective action of the theory.We proceed similarly as in the Majorana case.Therefore, the fermion density can be written as n = ⟨Q⟩ L .The Frame potential can be expressed as Since each of the Hamiltonian conserves the total number of particles Q = i ĉ † i ĉi , we can split each trace into smaller traces over each charge sector, namely with Namely, for each Keldysh contour and each replica we can introduce a chemical potential µ s ℓ , and we denote by µ the whole set of them.Writing down the path integral representation for each trace and by averaging over Hamiltonian realisations, we obtain with the action S(ψ, ψ) already introduced in eq.((2.19)) and with the projectors now written as Antiperiodic boundary conditions are enforced to give the expression of traces.By performing the Hubbard-Stratonovich transformation as usual eq.( 2.24) and integrating over Grassmann fields, we obtain the action: Here μ = µ s l δ ss ′ δ ll ′ is a matrix diagonal in the Keldysh contours space s, s ′ and in the replica spaces l, l ′ .And the fermionic two point function is given by i,l ′ (0) .We are interested in the equations of motion.Repeating the procedure from the section with Majorana fermions, one can find: In the first expression, we suppressed Keldysh and replica indexes, and in the second expression, we suppressed only replica indexes.Notice that for Majorana fermions, by definition, matrices Ĝ and Σ were real-valued.Now, they are allowed to be complex.In the case of single replica k = 1, in the limit t → 0 ± the Green function G ss ′ was found already in [45,46] but here we extend to generic replicas.The early-time behaviour can also be calculated similarly 3.1, such that the corresponding Frame potential is given by Now let us consider the late-time solution for q > 2. Notice that the last two equations of motions, eq.(4.10) define a boundary conditions for Ĝ, fixed by the charge content.Also, notice that we have the sum over all possible values of n +(−) in the eq.(4.7).Therefore, let us first consider the terms with n + l = n − l = n l .In this case, we have a replica diagonal solution, and it is given by charge charge-dependent matrix, see Appendix and with This matrix consists of the fermion densities from different replicas and Γ = (n(1 is also a diagonal matrix in replica space.As in the case of Majorana Fermions, other possible solutions of eq.(4.10) consists of permutation matrices, as these equations have permutation symmetry in the replica space.However, these solutions are valid only for the terms in the eq.( 4.10) which are given by permutation of the charge content n− = τ n+ τ T , due to the last two equations in eq.(4.10).This gives the solution: Let us calculate the on-shell action using this solution.We shall do our calculation for a replica diagonal solution, as the first three terms in the action are invariant under permutations.Therefore, we can always cast our solution to the replica diagonal one and get: where where Γ l is a matrix element of the diagonal matrix Γ.Which gives the answer for the determinant calculation: Now let us move to the second term in the action: 2 )δ(t) = 0, (4.22) notice that both side limit for the sum of Heaviside functions is not well-defined, therefore we assume: Therefore, as in the previous case, this term is equal to zero on the solution due to the forward and backward evolution, which comes with different signs.The third term in the action gives, We find therefore the final expression for the Frame potential where we used Q + = L k .We have then recovered the expression obtained by integration over the Haar group in each charge sector, see Appendix A.2.

4.2
The Gaussian case q = 2 Now let us consider the case q = 2: Similar to the case with q > 2 we start with a replica diagonal solution with n where f (t 12 ) = e − |t 12 | 1 2 −iμt 12 , since Γ = 1 in the case q = 2 and Σdiag = Since our equations of motion have rotation symmetry in replica space we can consider the solution with an arbitrary unitary matrix û: However notice that this solution implies boundary conditions n− = ûn + û † , but in the sum over all possible charge content, there are no terms that would satisfy these boundary conditions.Therefore, in the q = 2 case, this form of the solution does not respect boundary conditions from 4.26.Another possible solution is given by even charge content for each replica n = nδ ij = n 1.Notice that for this form of the solution, we have a generic rotated solution, with an arbitrary unitary matrix û.It is clear that boundary conditions are simplified in this case, as n − 1 = û(n + 1)û † = n + 1, and terms of this form are presented in the sum over the charge content.Therefore, the solution takes the form: where f (t 12 ) = e − |t 12 | 1 2 −iµ 1t 12 , with µ being just a number here.Here û is a unitary matrix in k × k space.Notice also that not all possible matrices û respect the parity of the initial state.Here we need to assume that det(û) = 1, in order to take into account the parity.Therefore, û matrix belongs to SU (k).As in the case of Majorana fermions, the difference between the two cases q > 2 and q = 2 will come from the consideration of the fluctuations around the saddle point solution and the presence of the Goldstone modes.The calculation of the action gives the same results as in the previous section and can be done the same way replacing n → n 1.

Symmetries, and fluctuations
Let us start with the consideration of the Goldstone theorem for the case q = 2.The vector fields ψ and ψ have 2k components, therefore the group acting in this space is U (2k), with d = 4k 2 generators.This is the number of degrees of freedom for our system.Now let us consider the group of the global symmetry of the action.The symmetries of the action are restricted by the term L log(det(σ z ∂ − Σ + iσ z μ)) and other terms with chemical potential.Notice that in the case n+ = n− = n 1 and μ+ = μ− = µ 1 the action again has rotation symmetries in the replica space U 1 = e i Â1 σz , U 2 = e i Â2 1, where Âi is a Hermitian matrix of size k.By Âσ z , we again mean the tensor product.Notice that matrices U 0 = e i Â form the U (k) group.Therefore, the group of symmetries of the action is B = U (k) × U (k).
Remember that the saddle point solution has the form eq. (4.32)where û is a unitary matrix.From the form of the solution, we see that the group that allows us to transform one solution into another is B/H = SU (k).Therefore, the number of massless modes is R ml = dim(B) − dim(H) = k 2 − 1.Then we have the number of massive modes: Let us consider fluctuations in our system.To do so, we again expand the effective action around the saddle point solution S = S 0 (G 0 , Σ 0 ) + S (2) (δ Ĝ, δ Σ) and integrate out the fluctuations, which are expected to give us the next order correction in the fermionic degrees of freedom.The integration over fluctuation in the highest order in L gives us The normalisation term can be calculated the same way as eq.(3.51).Which gives the scaling with the system size: Therefore in the final expression for the Frame potential, we get: and the logarithm is which coincides with the result for Gaussian Haar calculation, eq.(A.15).As in the case of Majorana fermions, for the generic case of q interactions, there are no Goldstone modes.Due to this fact, the integration over the fluctuations extends over all the degrees of freedom, which are dim(U (2k)).This coincides with the number of degrees of freedom at both zero and finite time.Therefore, the normalization will cancel the fluctuation part from the finite-time action.Hence, for the generic case of q > 2, we have: Notice that this result also coincides with the gHaar integrated result in eq.(A.11).

Conclusion and discussion
In conclusion, this paper explores the degree of quantum mixing given by the unitary time evolution of the Brownian SYK model, which serves as a paradigmatic model for quantum scrambling and chaos.Our investigation focuses on the Frame potential, a direct quantifier of k-th design realisation.By considering two identical copies of an initial density matrix and calculating the overlap between their independent time evolutions, the Frame potential, denoted as F (k) (T ), is evaluated.From a study of the short and late-time saddle points of the associated Keldysh action, we find that the Frame potential exhibits exponential decay at early times and saturates to an exponentially small, with the size of the system value at a timescale only dependent on the number q of interacting fermions.Our analysis reveals that in the case of non-Gaussian SYK models with q > 2, the Frame potential converges to its Haar-averaged value, denoted as F Haar , at late times.Such late-time behaviour manifests that the systems fully explore the available Hilbert space in a uniform manner.However, a difference arises in the Gaussian integrable case with two-fermion interactions.In this scenario, the evolution exhibits mixing behaviour solely within the space of Gaussian states.Here, a Gaussian Haar (gHaar) measure is introduced, and it is demonstrated that the Frame potential converges to the gHaar measure at late times, associated with the value of the Frame potential denoted as F Haar .Remarkably, the value of the gHaar-integrated Frame potential is equivalent to the generic case in the limit of a large fermion number L, albeit with logarithmic corrections in L. The presence of these corrections arises from Goldstone fluctuations around the late-time saddle point.Notably, these fluctuations are only present in the Gaussian case, where both the Keldysh action and the saddle point have a continuous group symmetry.In the complex (Dirac) SYK case, the total charge conservation forces the time evolution to reach Haar only within each charge sector.We have shown that by introducing projectors over the different charge sectors, one can introduce saddle point equations with different charge content between different replicas.While in the interacting q > 2 case the solution with different charges dominates, in the q = 2 case instead the one with all the replicas in the same sector dominates, as the latter possesses the largest number of massless modes.Overall, as expected, the constraint of the global charge conservation leads to values of the Frame potential at large times which are larger than the corresponding ones in the Majorana case.
The study of scrambling dynamics in the context of continuous monitoring has received significant attention in recent years, with a particular focus on measurement-induced phase transitions in fermionic systems [47][48][49][50].Exploring for example the effects of continuous monitoring on the Frame potential in the SYK model represents a promising avenue for extending the current work.Additionally, investigating the behaviour of the Frame potential in non-Brownian dynamics for the SYK model presents another compelling extension.orthogonal group for Majorana fermions and the unitary group for Dirac fermions.These Haar values correspond to the evolution with an arbitrary number of interacting fermions (q > 2).Second, we consider the case where the evolution matrix is quadratic in fermions (q = 2), i.e., U = e c † hc .In this case, the matrix values are again chosen with respect to the Haar measure, but the space of integration is smaller compared to the full Haar average.

A.1 Gaussian-Haar averages for Majorana fermions
If we consider Gaussian evolution, then its Haar distribution is the unitary evolution is a Gaussian U , i.e.U = e 1 2 h ij χ i χ j , where χ j are Majorana fields, then we can express Tr(U ) = Pf(1 + e h ij ) (as for example shown in [51]) where e h ij = û is an orthogonal matrix with the size L × L (L is even).We are interested in the case when matrices û are random with respect to the Haar measure.Notice that we are averaging over matrices with det(û) = 1, as û = e 1 2 h where h is an antisymmetric matrix (all the eigenvalues of h have conjugated pairs and lie on the imaginary axis, it means that eigenvalues of û are λ 1 = e ib 1 , λ * 1 = e −ib 1 , ... and their product is 1).This gives us the following averaging In the third line, we used the fact that all eigenvalues have their conjugated pair.The measure in this expression comes from the fact that we integrate other orthogonal matrices with the determinant equal to 1 and even size(see chapter 2.6 in [52]).Doing change of variables cos(θ) = t and calculating the Jacobian of this transformation we can get Now we need to calculate these integrals.For simplicity, we can introduce the function to name these integrals : Here we can notice that it is just a special case of Selberg integral (see chapters 3.6, 3.7, 4.1, 4.7 in [52]) where f ll ′ (t) is a time-dependent matrix in replica space with initial condition f ll ′ (0) = f (0)δ ll ′ .First, we consider the equation Σ ss ′ ll ′ (t) = c ss ′ q (2G ss ′ ll ′ ) q 2 −1 (t)(−2G s ′ s l ′ l ) q 2 (−t)δ(t), (B.2) which gives Σd = f q−1 (0) q where Γ = f q−1 (0) q Γ.Here we can use the formula for the block matrix and assume µ + l = µ − l = µ l (since Q + l = Q − l ), therefore we find

Figure 1 .
Figure1.The rescaled Frame potential with q = 2 in the logarithmic scale, exact numerical simulations compared with the theoretical predictions.The number of fermions in both cases is L = 120.Brown colour corresponds to the number of replicas k = 3, green colour corresponds to k = 2.The dotted line is the theoretical prediction of the Keldysh calculation at early times, and the dashed line is the prediction for the late time behaviour.From the numerical data, we can see that the time of the saturation can be estimated directly from the theoretical prediction, modulo finite size corrections (which are larger for larger k as expected).The intersection between the short-time saddle point (dotted line) and the long-time saddle point (dashed line) gives the estimation for the saturation time.In the case q = 2:T q=2 Dirac = 4 2 log(2) − k L log(L) + O 1 L , T q=2 Majorana = 4 2 log(2) − k−1 L log(L) + O1L ; as compared to the case q > 2, where we have T q>2 Dirac = T q>2 Majorana = q 2 • 2 log(2) + O 1 L .

. 17 )
Let us explain the form of the Hamiltonian.The derivative in the −− sector eq.(3.15) comes with an extra minus, therefore the Majorana fields in the partition function should come with an extra i.Finally, we obtain (Pf[σ z ∂ t − Σ]) L = Tr(e

24 )
Now we can finally sum over all the cycles, noticing that the size of the permutation matrix k should be equal to the sum of the length of cycles of this permutation k = c n cTr(e −T Ĥ ) =

. 42 )
Other possible transformations are not symmetries of the action, therefore the group is B = SO(k) × SO(k) and dim(B) = k(k − 1).Remind that the saddle point solution has the form Ĝ(t 1 , t 2 ) = e − |t 12 |