<figure><img src="image_1.png" /></figure>

HAL
open science

# Frame potential of Brownian SYK model of Majorana and Dirac fermions

Anastasiia Tiutiakina, Andrea de Luca, Jacopo de Nardis

► To cite this version:

Anastasiia Tiutiakina, Andrea de Luca, Jacopo de Nardis. Frame potential of Brownian SYK model of Majorana and Dirac fermions. Journal of High Energy Physics, 2024, 01, pp.115. 10.1007/JHEP01(2024)115. hal-04151996

HAL Id: hal-04151996

https://hal.science/hal-04151996v1

Submitted on 23 Oct 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_3.png" /></figure>

# Frame potential of Brownian SYK model of Majorana and Dirac fermions

Anastasiia Tiutiakina, Andrea De Luca and Jacopo De Nardis

Laboratoire de Physique Théorique et Modélisation, CNRS UMR 8089, CY Cergy Paris Université, 95302 Cergy-Pontoise Cedex, France

E-mail: anastasiia.tiutiakina@cyu.fr, andrea.de-luca@cyu.fr, jacopo.de-nardis@cyu.fr

**ABSTRACT:** 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.

**KEYWORDS:** Non-Equilibrium Field Theory, Quantum Dissipative Systems, Random Systems, Stochastic Processes

ARXIV ePRINT: 2306.11160
# Contents

<table><tr><td>1</td><td>Introduction and main results</td><td>1</td></tr><tr><td>2</td><td>Brownian SYK model and Keldysh path integral representation of the Frame potential</td><td>4</td></tr><tr><td>3</td><td>Majorana SYK model and saddle points</td><td>9</td></tr><tr><td>3.1</td><td>Short-time behaviour</td><td>9</td></tr><tr><td>3.2</td><td>Long-time behaviour</td><td>10</td></tr><tr><td>3.2.1</td><td>The non-Gaussian case $q > 2$</td><td>10</td></tr><tr><td>3.2.2</td><td>The Gaussian case $q = 2$</td><td>13</td></tr><tr><td>3.2.3</td><td>Symmetries and Goldstone modes</td><td>14</td></tr><tr><td>4</td><td>Dirac SYK model and Kelsdysh saddle points</td><td>17</td></tr><tr><td>4.1</td><td>Keldysh path integral and saddle points</td><td>18</td></tr><tr><td>4.2</td><td>The Gaussian case $q = 2$</td><td>21</td></tr><tr><td>4.2.1</td><td>Symmetries, and fluctuations</td><td>22</td></tr><tr><td>5</td><td>Conclusion and discussion</td><td>23</td></tr><tr><td>A</td><td>Haar-averaged Frame potential</td><td>24</td></tr><tr><td>A.1</td><td>Gaussian-Haar averages for Majorana fermions</td><td>24</td></tr><tr><td>A.2</td><td>Haar average of frame potential for Dirac fermions</td><td>25</td></tr><tr><td>A.3</td><td>Gaussian-Haar averages for Dirac fermions</td><td>26</td></tr><tr><td>B</td><td>Dirac SYK saddle point with fixed charge density</td><td>26</td></tr></table>

# 1 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–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–12]. Such a question is directly related to the emergence of eigenstate thermalization (ETH) in closed quantum systems [13–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–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–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–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 $\rho_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:

$$F^{(k)}(T) = \langle \text{Tr} (U_1(T) \rho_0 U_1^\dagger(T) U_2(T) \rho_0 U_2^\dagger(T))^k \rangle. \quad (1.1)$$

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 $\rho_0 = |\Psi\rangle\langle\Psi|$, 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 \to \infty$, it is expected that the density matrix $\rho(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 $\mathcal{H}$ as

$$\rho(T) \to U^\dagger \rho_0 U, \quad U \sim \text{Haar}(\mathcal{H}). \quad (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_{\text{Haar}}^{(k)}$ and it is easy to prove (see for instance [23]) that it corresponds to the minimum overall possible distributions of unitary operators on $\mathcal{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–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
<figure><img src="image_4.png" /><figcaption>**Figure 1.** 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_{\text{Dirac}}^{q=2} = 4(2\log(2) - \frac{k}{L}\log(L)) + O(\frac{1}{L})$, $T_{\text{Majorana}}^{q=2} = 4(2\log(2) - \frac{k-1}{L}\log(L)) + O(\frac{1}{L})$; as compared to the case $q > 2$, where we have $T_{\text{Dirac}}^{q>2} = T_{\text{Majorana}}^{q>2} = q^2 \cdot 2\log(2) + O(\frac{1}{L})$.</figcaption></figure>

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

$$F_{\text{Haar}}^{(k)} = e^{-\frac{LTk}{2q^2}} \quad \text{Majorana,}$$

$$F_{\text{Haar}}^{(k)} = e^{-\frac{LTk}{q^2}} \quad \text{Dirac.}$$

* For $q > 2$, given the dimension of the Hilbert space $\mathcal{N}$, the Frame potential saturates at large times to the values

$$F_{\text{Haar}}^{(k)} = \frac{k!}{\mathcal{N}^{2k}} \quad \text{Majorana,}$$

$$F_{\text{Haar}}^{(k)} = \frac{L^k k!}{\mathcal{N}^{2k}} \quad \text{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

$$F_{\text{gHaar}}^{(k)} = \frac{c_k L^{k(k-1)/2}}{\mathcal{N}^{2k}} + \dots \quad \text{Majorana,}$$

$$F_{\text{gHaar}}^{(k)} = \frac{\tilde{c}_k L^{k^2}}{\mathcal{N}^{2k}} + \dots \quad \text{Dirac,}$$

where $c_k = 2^{1-\frac{(k-1)(k-2)}{2}} \prod_{i=0}^{k-1} 1/(2i-1)!!$ and $\tilde{c}_k = \text{sf}^2(k-1)/\text{sf}(2k-1)$, with $\text{sf}(n) = 1!2!\dots n!$.

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 Swingar-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 figure 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 figure 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.

## 2 Brownian SYK model and Keldysh path integral representation of the Frame potential

We consider the time evolution generated by the time-dependent Hamiltonian

$$H(h,t) = i^{\frac{q}{2}} \sum_{1 \le i_1 < i_2 < \dots < i_q \le L} h_{i_1 \dots i_q}(t) \hat{\chi}_{i_1} \dots \hat{\chi}_{i_q}, \quad (2.1)$$

where $\hat{\chi}_i$ are $L$ Majorana fermions and anticommutation relations are applied $\{\hat{\chi}_i, \hat{\chi}_j\} = \delta_{ij}$, here $q$ is an even integer constant. This model has close analogies with the celebrated
SYK [24, 25, 34], but with the important difference that $h_{i_1...i_q}$ are not constant in time here. On the contrary, we take them to follow a white noise distribution,

$$ \langle h_i(t_1)h_j(t_2) \rangle = \frac{2^{q-1}q!}{q^2L^{q-1}} \delta_{i,j}\delta(t_1-t_2) = \sigma^2\delta_{i,j}\delta(t_1-t_2), \qquad (2.2) $$

where we denote the collective set of indices $i = (i_1, ..., i_q)$ (similarly for $j$) and set $\delta_{i,j} = \delta_{i_1,j_1} ... \delta_{i_q,j_q}$. For this reason, the current model is named Brownian SYK. We shall here assume $L$ to be an even integer, such that the operators $\hat{\chi}_i$ can be embedded into a $\mathcal{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, 35, 36], 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 $\hat{\chi}_{j,L}$ and $\hat{\chi}_{j,R}$, where $j$ takes values from 1 to $L$. (It should be noted that $L$ must be an even number to ensure a well-defined Hilbert space). A convenient specification for the state can be formulated using Dirac fermions

$$ (\hat{c}_{j,L} - \hat{c}_{j,R}^\dagger) |\text{TFD}\rangle = 0 \quad \text{and} \quad (\hat{c}_{j,L}^\dagger - \hat{c}_{j,R}) |\text{TFD}\rangle = 0, \qquad (2.3) $$

which can be expressed in terms of Majorana fermions as

$$ \hat{c}_{j,L} = \frac{1}{\sqrt{2}} (\hat{\chi}_{2j-1,L} + i\hat{\chi}_{2j,L}) \quad \text{and} \quad \hat{c}_{j,R} = \frac{1}{\sqrt{2}} (\hat{\chi}_{2j,R} - i\hat{\chi}_{2j-1,R}). \qquad (2.4) $$

Then we have the standard definition in terms of Majorana fermions:

$$ (\hat{\chi}_{j,L(R)} + i\hat{\chi}_{j,L(R)}) |\text{TFD}\rangle = 0, \quad \forall j = 1, \dots, L. \qquad (2.5) $$

From the eq. (2.3) we see that the TFD state in the eigenbasis of $n_{j,L(R)} = \hat{c}_{j,L(R)}^\dagger \hat{c}_{j,L(R)}$ basis can be expressed as a product of $|0\rangle_L |0\rangle_R + |1\rangle_L |1\rangle_R$ on each site. More explicitly, we consider the initial state:

$$ |\text{TFD}\rangle = \frac{1}{\sqrt{N}} \sum_{j=1}^{N} |j\rangle |j\rangle, \quad \rho_0 = |\text{TFD}\rangle \langle \text{TFD}|, \qquad (2.6) $$

where each $|j\rangle$ represents a possible string of zeros and ones. Here $N$ is the dimension of the Hilbert space. This state has the following property,

$$ \mathrm{Tr}^{(2)}[A \otimes B |\text{TFD}\rangle \langle \text{TFD}| C \otimes D] = \frac{1}{N} \mathrm{Tr}[AB^t D^t C], \qquad (2.7) $$

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,

$$|\text{TFD}\rangle(t) = U_{h^L}(t) \otimes U_{h^R}(t) |\text{TFD}\rangle, \quad \rho(t) = |\text{TFD}(t)\rangle \langle \text{TFD}(t)|, \quad (2.8)$$

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

$$U_{h^\sigma}(t+dt) = e^{-iH(h^\sigma,t)dt}U_{h^\sigma}(t), \quad (2.9)$$

where the subscript $\sigma = 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 $\hat{\chi}_{i,\sigma}$ and independently generated white noises $h_i^\sigma(t)$. Starting with this initial state, we calculate the $k$-th moments of the Frame potential averaged over Hamiltonian realisations with measure $d\eta(\mathbf{h})$:

$$F^{(k)}(T) = \int d\eta(\mathbf{h}_1) d\eta(\mathbf{h}_2) \text{Tr} \left( U_{h_1}(T) \rho_0 U_{h_1}^{\dagger}(T) U_{h_2}(T) \rho_0 U_{h_2}^{\dagger}(T) \right)^k, \quad (2.10)$$

here, $U_{h_1}(T) = U_{h_1^L}(T) \otimes U_{h_1^R}(T)$, (and the same for $h_2^{L(R)}$), where the couplings $h_{1(2)}^{L(R)}$ are independent random variables with the distribution $d\eta(\mathbf{h}) = e^{-\frac{\mathbf{h}^2}{2\sigma^2}} d\mathbf{h}$ and variance defined in eq. (2.2). The initial condition $\rho_0$ is the double thermofield state that we discussed earlier. Then, the Frame potential takes the form, for generic $k$,

$$F^{(k)}(T) = \frac{1}{\mathcal{N}^{2k}} \int d\eta(\mathbf{h}) |\text{Tr}[U_{\mathbf{h}}(T)]|^{2k}, \quad (2.11)$$

where the evolution operator $U_{\mathbf{h}}(T)$ is now only acting on a single copy $R$ (or $L$) and where $\mathcal{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 $U_{\mathbf{h}} = U_{h_1^L}(t)U_{h_1^R}(t)U_{h_2^R}^\dagger(t)U_{h_2^L}^\dagger(t)$. 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 [37]. To do so, we introduce an additional set of Majorana operators $\eta_i$ and create fermionic operators in the form:

$$\hat{c}_i = \frac{\hat{\chi}_i + i\hat{\xi}_i}{\sqrt{2}}, \quad \hat{c}_i^\dagger = \frac{\hat{\chi}_i - i\hat{\xi}_i}{\sqrt{2}}. \quad (2.12)$$

With this notation, the coherent state can be defined as

$$|\psi\rangle = |0\rangle - \psi |1\rangle = e^{-\psi \hat{c}^\dagger |0\rangle}, \quad (2.13)$$

where $\psi$ is a Grassmann variable that parameterizes our state. The set of these states is not orthonormal, and their overlap is given by:

$$\langle\psi|\psi'\rangle = (\langle 0| - \langle 1|\bar{\psi}\rangle)(\langle 0| - \psi'\langle 1|) = 1 + \bar{\psi}\psi' = e^{\bar{\psi}\psi'}. \quad (2.14)$$
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:

$$ \int d\bar{\psi} \int d\psi e^{-\bar{\psi}\psi} |\psi\rangle\langle\psi| = 1, \qquad (2.15) $$

and secondly we can express the trace of an operator $\hat{O}$ in terms of Grassmann variables as:

$$ \mathrm{Tr}(\hat{O}) = \sum_{n_1, n_2, \dots = 0, 1} \langle n_1 \dots n_N | \hat{O} | n_1 \dots n_N \rangle = \int d\bar{\psi} d\psi e^{-\bar{\psi}\psi} \langle \psi | \hat{O} | -\psi \rangle, \qquad (2.16) $$

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

$$ F^{(k)} = \frac{1}{\mathcal{N}^{2k}} \int d\eta(h) d\bar{\psi} d\psi \prod_{l=1}^{k} \langle \psi_1^{-,l} | U_1(\epsilon) | \psi_2^{-,l} \rangle \langle \psi_2^{-,l} | U_2(\epsilon) | \psi_3^{-,l} \rangle \dots \langle \psi_{2N}^{-,l} | - \psi_1^{-,l} \rangle \cdot \langle \psi_{2N}^{+,l} | U_{2N}^{\dagger}(\epsilon) | \psi_{2N-1}^{+,l} \rangle \langle \psi_{2N-1}^{+,l} | U_{2N-1}^{\dagger}(\epsilon) | \psi_{2N-2}^{+,l} \rangle \dots \langle \psi_1^{+,l} | - \psi_{2N}^{+,l} \rangle e^{-\bar{\psi}_i^{s,l} \psi_i^{s,l}}, \qquad (2.17) $$

where $s = \pm$ indicates one of the two Keldysh contours, $l$ counts replicas, and $i$ labels time. Each matrix element here has the form,

$$ \langle \psi_1^{-,l} | U_1(\epsilon) | \psi_2^{-,l} \rangle = \langle \psi_1^{-,l} | e^{-i\frac{q}{2}+1} \epsilon h_{i_1 i_2 \dots i_q} (\hat{c}_{i_1}^{\dagger} + \hat{c}_{i_1}) \dots (\hat{c}_{i_q}^{\dagger} + \hat{c}_{i_q}) | \psi_2^{-,l} \rangle. \qquad (2.18) $$

Let us notice that for any matrix element of this form, we can simply replace $\hat{c}^{\dagger}$ with $\bar{\psi}$ and $\hat{c}$ with $\psi$. 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 \neq i_l$ everywhere and all fermion operators anti-commute), obtain the action, and once again reverse the order of the corresponding Grassmann variables $\psi, \bar{\psi}$ at the cost of an extra minus sign. Therefore, we shall write the usual Keldysh action:

$$ S(\bar{\psi}, \psi) = \sum_{l,s=\pm} s \int_0^T dt \left( \sum_j \bar{\psi}_j^{s,l} \partial_t \psi_j^{s,l} + H(\bar{\psi}, \psi) \right). \qquad (2.19) $$

Here the continuous limit was taken $\epsilon \to 0$, $N \to \infty$, and $N\epsilon \to T$. We can now move back to Majorana fermions:

$$ \chi_i^{s,l} = \frac{\bar{\psi}_i^{s,l} + \psi_i^{s,l}}{\sqrt{2}}, \qquad \xi_i^{s,l} = i \frac{\bar{\psi}_i^{s,l} - \psi_i^{s,l}}{\sqrt{2}}, \qquad (2.20) $$

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 $\chi$ and $\xi$ fields separate [37], and the partition function factorizes. Therefore, the part that depends on $\eta$ 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

$$ F^{(k)} \equiv \frac{\tilde{F}^{(k)}(T)}{\tilde{F}^{(k)}(T=0)} = \frac{Z(\xi)}{\tilde{F}^{(k)}(T=0) \mathcal{N}^{2k}} \int d\eta(h) d\chi e^{-S(\chi)}, \qquad (2.21) $$
where we have

$$S(\chi) = \sum_{l,s=\pm} s \int_0^T dt \left( \frac{1}{2} \sum_j \chi_j^{s,l} \partial_t \chi_j^{s,l} + i^{\frac{q}{2}+1} \sum_{i_1<i_2<\dots<i_q} h_{i_1\dots i_q}(t) \chi_{i_1}^{s,l} \dots \chi_{i_q}^{s,l} \right), \quad (2.22)$$

and $\tilde{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, 34, 38]. 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 $\sigma^2$:

$$\int dh_{i_1 \dots i_q} e^{-\frac{1}{2\sigma^2} h_{i_1 \dots i_q}^2 + i^{\frac{q}{2}+1} h_{i_1 \dots i_q} M_{i_1 \dots i_q}} = e^{-\frac{(-1)^{\frac{q}{2}} \sigma^2}{2} M_{i_1 \dots i_q}^2}, \quad (2.23)$$

where $M_{i_1 \dots i_q} = \sum_l (\chi_{i_1}^{-,l} \dots \chi_{i_q}^{-,l} - \chi_{i_1}^{+,l} \dots \chi_{i_q}^{+,l})$. Let us introduce the bi-local fields $\hat{G}(t_1, t_2)$, $\hat{\Sigma}(t_1, t_2)$ and perform Hubbard-Stratonovich transformation:<fnref n="1" />

$$\int d\hat{G} \prod_{s,s'} \delta \left( G_{ll'}^{ss'}(t_1, t_2) - \frac{1}{L} \sum_i \chi_i^{s,l}(t_1) \chi_i^{s',l'}(t_2) \right) = \int d\hat{\Sigma} d\hat{G} e^{-\frac{L}{2} \int dt_1 dt_2 \sum_{ll'}^{ss'} \left( G_{ll'}^{ss'} - \frac{1}{L} \chi_i^{s,l} \chi_i^{s',l'} \right)} \quad (2.24)$$

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:

$$S = - \sum_{s,l} s \int_0^T \frac{1}{2} \chi_j^{s,l} \partial_t \chi_j^{s,l} + \frac{L}{4q^2} \sum_{l,l',s,s'} \int_0^T dt_1 dt_2 c_{ss'} \left( 2G_{ll'}^{ss'} \right)^q \delta(t_{12}) - \frac{L}{2} \sum_{l,l',s,s'} \int dt_1 dt_2 \sum_{ll'}^{ss'} (t_1, t_2) \left( G_{ll'}^{ss'}(t_1, t_2) - \frac{1}{L} \chi_i^{s,l}(t_1) \chi_i^{s',l'}(t_2) \right), \quad (2.25)$$

where

$$c_{ss'} = -ss'. \quad (2.26)$$

Notice that in order to form bi-local fields in the expression eq. (2.23), we performed $\frac{q(q-1)}{2}$ transpositions of the Grassman variables, which gave us the prefactor in the interaction term $(-1)^{\frac{q(q-1)}{2}}$ together with $(-1)^{\frac{q}{2}}$ from the Gaussian integration, resulting in $(-1)^{\frac{q^2}{2}}$, which is unity due to considering $q$ as an even integer. In order to proceed further, we integrate out the fields $\chi_i$, and the action takes the form:

$$S = L \log(\mathrm{Pf}(\hat{\sigma}_z \partial_t - \hat{\Sigma})) + \frac{L}{4q^2} \sum_{l,l',s,s'} \int_0^T dt_1 dt_2 c_{ss'} \left(2G_{ll'}^{ss'}\right)^q \delta(t_{12}) - \frac{L}{2} \mathrm{Tr}(\hat{\Sigma}^T \hat{G}), \quad (2.27)$$

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

<footnote n="1">¹Notice that here in the case of Brownian SYK one could also more simply introduce fields $\hat{G}(t)$ of only one time instead of two, as the interactions are local in time.</footnote>
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 $\frac{1}{2}\chi\partial_t\chi$ has a positive sign for the $+$ contour and a negative sign for the $-$ contour. Therefore, this term can also be expressed as a $2 \times 2$ matrix with Keldysh contour indexes in the following way:

$$ -\sum_{s,l} s \int_0^T \frac{1}{2} \sum_j \chi_j^{s,l} \partial_t \chi_j^{s,l} = -\frac{1}{2} \sum_l \sum_j \chi_j^l (\hat{\sigma}_z \partial_t) \chi_j^l. \qquad (2.28) $$

Here, we introduced the Pauli matrix $\hat{\sigma}_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 $\hat{G}$ and $\hat{\Sigma}$ using the saddle-point approximation at large $L$.

### 3 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:

$$ \frac{\delta S}{\delta \hat{\Sigma}} : \hat{G}^{-1}(t_1, t_2) = \hat{\sigma}_z \partial_t - \hat{\Sigma}(t_1, t_2), \qquad (3.1) $$

$$ \frac{\delta S}{\delta \hat{G}} : \hat{\Sigma}(t_1, t_2) = \frac{\hat{c}}{q} (2\hat{G})^{q-1}(t_1, t_2)\delta(t_{12}), \qquad (3.2) $$

with the tensor $\hat{c}_{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.

#### 3.1 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:

$$ \hat{\Sigma} = 0, \qquad \hat{G}^{-1} = \hat{\sigma}_z \partial_t. \qquad (3.3) $$

With this solution, the Frame potential can be expressed as:

$$ F^{(k)}(T=0) = \frac{1}{N^{2k}} \det(\hat{\sigma}^z \partial_t)^{\frac{L}{2}} e^{-\frac{L}{4q^2} \sum_{ll'} \int_0^T dt (2G_{ll'}^{++})^q + (2G_{ll'}^{--})^q} . \qquad (3.4) $$

This functional determinant can be calculated using diagonalisation of the operator $\hat{\sigma}^z \partial_t$. By establishing antiperiodic boundary conditions, we obtain eigenvalues $\lambda_n = \pm \frac{\pi i (2n+1)}{T}$, where $n = 0, \pm 1, \pm 2, \dots$. In the $++$ space, we have two eigenvalues $\lambda_n$ and $-\lambda_n$. However, coherent states are over complete, therefore we consider $\lambda_n = \frac{\pi i}{T}$, $\frac{3\pi i}{T}$, ..., $\frac{(4N+1)\pi i}{4T}$, which prevents us from treating matrices

$$ \begin{pmatrix} -\lambda_n & 0 \\ 0 & \lambda_n \end{pmatrix} \quad \text{and} \quad \begin{pmatrix} \lambda_n & 0 \\ 0 & -\lambda_n \end{pmatrix}, \qquad (3.5) $$
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 \to \infty$ and following the regularization technique from the book [39] (notice that we have $\omega = 0$), we get

$$ \det(\hat{\sigma}_z \partial_t)^{\frac{L}{2}} = 2^{\frac{2kL}{2}}, \qquad (3.6) $$

where $2^{\frac{L}{2}}$ is the dimension of the Hilbert space for the initial operator $U$. Notice that the field $\hat{G}$ is also diagonal in the replica space. Therefore, $G_{ll'}^{++} = \frac{1}{2}\delta_{ll'} = -G_{ll'}^{--}$, and we have

$$ e^{-\frac{L}{4q^2} \sum_{ll'} \int_0^T dt (2G_{ll'}^{++})^q + (2G_{ll'}^{--})^q} = e^{-\frac{LTk}{2q^2}}. \qquad (3.7) $$

The Frame potential on this solution takes the form

$$ F^{(k)} = e^{-\frac{LTk}{2q^2}}. \qquad (3.8) $$

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.

## 3.2 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) \times SO(k)$ (i.e. the replica indices in each Keldysh contour can be rotated separately), therefore the fields $\hat{\Sigma}$ and $\hat{G}$ can be rotated in the replica space by the generators $O_1 = e^{i\hat{A}_1\hat{\sigma}_z}$, $O_2 = e^{i\hat{A}_2\hat{1}}$ (which form $SO(k) \times SO(k)$ group), here $\hat{A}_i$ is a skew-symmetric matrix of the size $k$ and by $\hat{A}\hat{\sigma}_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

$$ \hat{G}(t_1, t_2) = \frac{f(|t_{12}|)}{2} \begin{pmatrix} \Theta(t_{12}) - \Theta(-t_{12}) & -\hat{\tau}^T \\ \hat{\tau} & -\Theta(t_{12}) + \Theta(-t_{12}) \end{pmatrix}. \qquad (3.9) $$

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 $\Theta(t_{12})$ is a Heaviside theta function. The matrix $\hat{\tau}$ is a permutation matrix in the replica space. Let us discuss the function $\text{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($\epsilon \times (i - i + 1)$) = sgn($\epsilon$). 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:

$$ \hat{\Sigma}(t_1, t_2) = \frac{1}{q} \delta(t_{12}) f(0)^{q-1} \begin{pmatrix} 0 & -\hat{\tau}^T \\ \hat{\tau} & 0 \end{pmatrix}. \qquad (3.10) $$

Now, our objective is to determine the function $f(t)$. To achieve this, we perform the Fourier transformation $\hat{\Sigma}(t) = \int d\omega \hat{\Sigma}(\omega) e^{-i\omega t}$ and substitute this ansatz into the first equation of motion (3.1). This procedure gives:

$$ \hat{G}(\omega) = \frac{1}{\omega^2 + \left(\frac{f(0)^{q-1}}{q}\right)^2} \begin{pmatrix} i\omega & -\frac{f(0)^{q-1}}{q} \hat{\tau}^T \\ \frac{f(0)^{q-1}}{q} \hat{\tau} & -i\omega \end{pmatrix}. \qquad (3.11) $$

By performing the inverse Fourier transform and comparing it with the ansatz, we arrive at the following result:

$$ \hat{G}(t_1, t_2) = \frac{e^{-\frac{|t_{12}|}{q}}}{2} \begin{pmatrix} \Theta(t_{12}) - \Theta(-t_{12}) & -\hat{\tau}^T \\ \hat{\tau} & -\Theta(t_{12}) + \Theta(-t_{12}) \end{pmatrix}, \qquad \hat{\Sigma}(t_1, t_2) = \frac{1}{q} \delta(t_{12}) \begin{pmatrix} 0 & -\hat{\tau}^T \\ \hat{\tau} & 0 \end{pmatrix}. \qquad (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

$$ \frac{L}{4q^2} \sum_{l,l's,s'} \int_0^T dt_1 dt_2 c_{ss'} \left( 2G_{ll'}^{ss'} \right)^q \delta(t_{12}) = 0. \qquad (3.13) $$

On the other hand, the third term will give us a non-trivial contribution

$$ -\frac{L}{2} \sum_{l,l's,s'} \int_0^T dt_1 dt_2 \Sigma_{ll'}^{ss'}(t_1, t_2) G_{ll'}^{ss'}(t_1, t_2) = -\frac{L}{2} \text{Tr}(\hat{\Sigma}^T \hat{G}) = -\frac{LTk}{2q}. \qquad (3.14) $$

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$:

$$ e^{S_1} = (\text{Pf}[\hat{\sigma}_z \partial_t - \hat{\Sigma}])^L = \text{Tr}(e^{-T\hat{H}})^L, \qquad (3.15) $$

where

$$ \hat{H} = -\frac{1}{2} (\hat{\chi}^+ i\hat{\chi}^-) \hat{\Sigma}_0 \begin{pmatrix} \hat{\chi}^+ \\ i\hat{\chi}^- \end{pmatrix}, \qquad (3.16) $$

and

$$ \hat{\Sigma}_0 = \frac{1}{q} \begin{pmatrix} 0 & -\hat{\tau}^T \\ \hat{\tau} & 0 \end{pmatrix}. \qquad (3.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

$$ (\mathrm{Pf}[\hat{\sigma}_z \partial_t - \hat{\Sigma}])^L = \mathrm{Tr}(e^{i \frac{T}{q} \hat{x}^- \hat{\tau} \hat{x}^+})^L. \quad (3.18) $$

A similar expression was calculated in [40], and here we use the analogous method. Let us observe that any arbitrary permutation can be expressed as a direct sum of disjoint cycles, $\hat{\tau} = \bigoplus_i \hat{\tau}_i^c$ which implies,

$$ \hat{\tau} = \begin{pmatrix} \hat{\tau}_1^c & 0 & 0 & 0 \\ 0 & \hat{\tau}_2^c & 0 & 0 \\ 0 & 0 & \dots & 0 \\ 0 & 0 & 0 & \hat{\tau}_m^c \end{pmatrix}, \quad (3.19) $$

then $\mathrm{Tr}(e^{i \frac{T}{q} \hat{x}^- \hat{\tau} \hat{x}^+}) = \mathrm{Tr}(e^{i \frac{T}{q} \sum_i \hat{x}_i^- \hat{\tau}_i^c \hat{x}_i^+})$, where the sum goes over the cycles. Also, we notice that each cycle can be represented in the canonical form, therefore we have

$$ \hat{\tau}_i^c = \begin{cases} \delta^{\alpha+1,\beta}, & n_c \text{ odd} \\ \operatorname{sgn}(\beta - \alpha)\delta^{\alpha+1,\beta}, & n_c \text{ even,} \end{cases} \quad (3.20) $$

where $n_c$ represents the length of a cycle, and $\delta^{n_c+1,\beta} = \delta^{1,\beta}$. Here, $\operatorname{sgn}(\beta - \alpha)$ 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 $\hat{\tau}$, 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:

$$ -T\hat{H}_c = i \frac{T}{q} \left( \hat{\chi}_\alpha^- \delta^{\alpha+1,\beta} \hat{\chi}_\beta^+ + (-1)^{n_c+1} \hat{\chi}_{n_c}^- \hat{\chi}_1^+ \right). \quad (3.21) $$

Introducing complex fermions

$$ c_{\alpha} = \frac{\hat{\chi}_{\alpha}^{+} + i\hat{\chi}_{\alpha}^{-}}{\sqrt{2}}, \qquad c_{\alpha}^{\dagger} = \frac{\hat{\chi}_{\alpha}^{+} - i\hat{\chi}_{\alpha}^{-}}{\sqrt{2}}, \quad (3.22) $$

and inserting this into the Hamiltonian we obtain

$$ -T\hat{H}(\hat{\tau}_c) = \frac{T}{2q} \sum_{\alpha} (c_{\alpha}^{\dagger} c_{\alpha+1} + c_{\alpha} c_{\alpha+1} + \mathrm{h.c.}) + \frac{T}{2q} ((-1)^{n_c+1} c_{n_c}^{\dagger} c_1 + \mathrm{h.c.}), \quad (3.23) $$

this is the Kitaev chain [41] with periodic/antiperiodic boundary conditions, eq. (3.20) and odd/even length of chain respectively. The diagonalisation of the chain is readerly done, giving

$$ \mathrm{Tr}(e^{-T\hat{H}(\hat{\tau}_c)}) \sim e^{\frac{n_c T}{2q}}, \quad T \to \infty. \quad (3.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 = \sum_c n_c$

$$ \mathrm{Tr}(e^{-T\hat{H}}) = \prod_c e^{\frac{n_c T}{2q}} = e^{\frac{kT}{2q}}. \quad (3.25) $$
Consequently we find the following expression for the Pfaffian:

$$e^{S_1} = (\mathrm{Pf}[\hat{\sigma}_z \partial_t - \Sigma])^L = e^{\frac{LkT}{2q}}, \qquad (3.26)$$

and the Frame potential of this solution can be expressed as:

$$F^{(k)}(\hat{\tau}) = \frac{1}{\mathcal{N}^{2k}} e^{S_1 + S_2 + S_3} = \frac{1}{\mathcal{N}^{2k}} e^{\frac{LkT}{2q} - \frac{LkT}{2q}} = \frac{1}{\mathcal{N}^{2k}}. \qquad (3.27)$$

We shall now perform the summation of all the possible saddle points, which is the space of permutation $\hat{\tau}$ 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$,

$$F^{(k)} = \sum_{\hat{\tau}} F^{(k)}(\hat{\tau}) = \frac{k!}{\mathcal{N}^{2k}}, \qquad (3.28)$$

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.

### 3.2.2 The Gaussian case $q = 2$

In the case where $q = 2$, the action has the continuous symmetry, $SO(k) \times 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:

$$\hat{G}(t_1, t_2) = \frac{e^{-\frac{|t_1 t_2|}{q}}}{2} \begin{pmatrix} \operatorname{sgn}(t_{12}) & -\hat{\theta}^T \\ \hat{\theta} & -\operatorname{sgn}(t_{12}) \end{pmatrix}, \qquad \hat{\Sigma}(t_1, t_2) = \frac{1}{q} \delta(t_{12}) \begin{pmatrix} 0 & -\hat{\theta}^T \\ \hat{\theta} & 0 \end{pmatrix}, \qquad (3.29)$$

where $\hat{\theta}$ 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:

$$\hat{\theta}_{2m} = \begin{pmatrix} \hat{R}_1 & 0 & 0 & 0 \\ 0 & \hat{R}_2 & 0 & 0 \\ 0 & 0 & \dots & 0 \\ 0 & 0 & 0 & \hat{R}_{m-1} \\ 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 \end{pmatrix}, \quad \text{and} \quad \hat{\theta}_{2m+1} = \begin{pmatrix} \hat{R}_1 & 0 & 0 & 0 & 0 \\ 0 & \hat{R}_2 & 0 & 0 & 0 \\ 0 & 0 & \dots & 0 & 0 \\ 0 & 0 & 0 & \hat{R}_m & 0 \\ 0 & 0 & 0 & 0 & \pm 1 \end{pmatrix}, \qquad (3.30)$$

where, $m$ is a positive integer, and $\hat{R}_i$ is an orthogonal matrix in two-dimensional space. Then, $\mathrm{Tr}(e^{i \frac{T}{q} \hat{\chi}^{-} \hat{\theta}^{+} \hat{\chi}^{+}}) = \mathrm{Tr}(e^{i \frac{T}{q} \sum_{i} \hat{\chi}_{i}^{-} \hat{R}_{i} \hat{\chi}_{i}^{+}})$. Here, we should also take into account that not all orthogonal matrices will preserve the parity of the initial state. Therefore, $\hat{R}_i$ 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:

$$-T\hat{H}_R = i\frac{T}{q}(\cos(\phi)\hat{\chi}_1^-\hat{\chi}_1^+ + \sin(\phi)\hat{\chi}_1^-\hat{\chi}_2^+ - \sin(\phi)\hat{\chi}_2^-\hat{\chi}_1^+ + \cos(\phi)\hat{\chi}_2^-\hat{\chi}_2^+) . \quad (3.31)$$

Here $\phi$ 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:

$$-T\hat{H}_R = \frac{T}{2q}(2\cos(\phi)(c_1^\dagger c_1 + c_2^\dagger c_2 - 1) + 2\sin(\phi)(c_1 c_2 - c_1^\dagger c_2^\dagger)) = \\ = \frac{T}{q}(c_1^\dagger c_2)(\cos\phi - \sin\phi - \sin\phi - \cos\phi) \begin{pmatrix} c_1 \\ c_2^\dagger \end{pmatrix},$$

after the diagonalisation the trace can be easily calculated. Using the well-known formula for the trace of quadratic density matrices $\mathrm{Tr}(e^{c_i^\dagger \Gamma_{ij} c_j}) = \det(1 + e^{\Gamma_{ij}})$, we obtain:

$$\mathrm{Tr}(e^{-T\hat{H}_R}) = \det \left( 1 + \exp \left( \frac{T}{q} \begin{pmatrix} 1 & 0 \\ 0 & -1 \end{pmatrix} \right) \right) = 2 + e^{\frac{T}{q}} + e^{-\frac{T}{q}} \rightarrow e^{\frac{T}{q}}, \quad T \rightarrow \infty. \quad (3.32)$$

Therefore the Frame potential on the particular solution which depends on an orthogonal matrix $\hat{\theta}$ is

$$F^{(k)}(\hat{\theta}) = \frac{1}{\mathcal{N}^{2k}} e^{\frac{LkT}{2q} - \frac{LkT}{2q}} = \frac{1}{\mathcal{N}^{2k}}. \quad (3.33)$$

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 $\frac{\log(L)}{L}$ corrections, while with $q > 2$ there are no zero-mass modes that can fluctuate around the saddle point solution.

### 3.2.3 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:

$$G_{ll'} = \begin{pmatrix} G_{ll'}^{++} & G_{ll'}^{+-} \\ G_{ll'}^{-+} & G_{ll'}^{--} \end{pmatrix}, \quad (3.34)$$

and $G_{ll'}^{++} = \frac{1}{L} \sum_{i=1}^{L} (\chi^+)_i^i (\chi^+)_i^{i'}$, 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 $\chi$, which contains $2k(2k-1)/2$ degrees of freedom. Let us check the commutation relations for this matrix

$$[G_{ll'}, G_{mm'}] = \frac{1}{L^2} \sum_i \sum_j [\chi_l^i \chi_{l'}^i, \chi_m^j \chi_{m'}^j], \qquad (3.35)$$

expanding the commutator and using anticommutation relations $\{\chi_l^i, \chi_m^j\} = \frac{1}{2}\delta^{ij}\delta^{lm}$ we have,

$$[G_{ll'}, G_{mm'}] = \frac{1}{2L} (\delta_{l'm} G_{lm'} - \delta_{lm'} G_{ml} + \delta_{lm} G_{m'l'} - \delta_{l'm} G_{lm}), \qquad (3.36)$$

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 [42], it was shown that $S_{ll'} = iG_{ll'}$ satisfy orthogonality relation in the large $L$ limit:

$$SS^T = \mathbf{0}, \qquad (3.37)$$

which for the matrix $\hat{G}$ this implies

$$\hat{G}\hat{G}^T = \mathbf{0}. \qquad (3.38)$$

Therefore this matrix should be not only skew-symmetric but also orthogonal. The space of these matrices isomorphic to the $SO(2k)/U(k)$ [43], which gives us $d = \frac{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 $\hat{G}$ we have

$$S = L \log(\mathrm{Pf}(\hat{\sigma}_z \partial_t - \hat{\Sigma})) - \frac{L}{4} \sum_{l,l's,s'} \int_0^T dt_1 dt_2 c_{ss'} (\sum_{ll'}^{ss'})^2 \delta(t_{12}). \qquad (3.39)$$

One can see that the problematic term here is $L \log(\mathrm{Pf}(\hat{\sigma}_z \partial_t - \hat{\Sigma}))$, which as we already noticed has a rotation symmetry in replica space $O_1 = e^{i\hat{A}_1\hat{\sigma}_z}$, $O_2 = e^{i\hat{A}_2\hat{\sigma}_z}$ with $\hat{A}_i$ skew-symmetric matrix of the size $k$ and by $\hat{A}\hat{\sigma}_z$ we suppose the tensor product. Notice that matrices $O_0 = e^{i\hat{A}}$ form a $SO(k)$ symmetry group, therefore the space of the symmetries is $SO(k) \times SO(k)$. Therefore,

$$\det(\hat{\sigma}_z \partial_t - O_i \hat{\Sigma} O_i^T) = \det(O_i O_i^T) \det(O_i^{-1} \hat{\sigma}_z \partial_t (O_i^T)^{-1} - \hat{\Sigma}), \qquad (3.40)$$

since matrices $O_i$ commute with $\hat{\sigma}_z$, we have

$$\det(\hat{\sigma}_z \partial_t - O_i \hat{\Sigma} O_i^T) = \det(\hat{\sigma}_z \partial_t - \hat{\Sigma}). \qquad (3.41)$$

Other possible transformations are not symmetries of the action, therefore the group is $B = SO(k) \times SO(k)$ and $\dim(B) = k(k-1)$. Remind that the saddle point solution has the form

$$\hat{G}(t_1, t_2) = \frac{e^{-\frac{|t_{12}|}{2}}}{2} \begin{pmatrix} \mathrm{sgn}(t_{12}) & -\hat{\theta}^T \\ \hat{\theta} & -\mathrm{sgn}(t_{12}) \end{pmatrix}, \qquad \hat{\Sigma} = \frac{1}{2} \delta(t_{12}) f(0) \begin{pmatrix} 0 & -\hat{\theta}^T \\ \hat{\theta} & 0 \end{pmatrix}, \qquad (3.42)$$
where $\hat{\theta}$ is an orthogonal matrix $k \times 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, $R_{ml} = \dim(B) - \dim(H) = k(k-1)/2$. Therefore, we find

$$R_m = d - R_{ml} = k(k-1) - k(k-1)/2 = k(k-1)/2, \quad (3.43)$$

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

$$\tilde{F}^{(k)}(T) = \frac{1}{\mathcal{N}^{2k}} e^{S_0} \left\langle \int D\delta G D\delta\Sigma e^{\frac{1}{2}S^{(2)}(\delta G, \delta\Sigma)} \right\rangle_{\hat{\theta}}, \quad (3.44)$$

where the on-shell action $S_0$ was found in previous sections and $\langle \dots \rangle_{\hat{\theta}}$ means the summation over all the solutions. Now let us consider the quadratic fluctuation term $S^{(2)}(\delta G, \delta\Sigma)$:

$$S^{(2)}(\delta G, \delta\Sigma) = \frac{\delta^2 S}{\delta\hat{G}\delta\hat{G}} \delta\hat{G}^2 + \frac{\delta^2 S}{\delta\hat{\Sigma}\delta\hat{\Sigma}} \delta\hat{\Sigma}^2 + \frac{\delta^2 S}{\delta\hat{G}\delta\hat{\Sigma}} \delta\hat{\Sigma}\delta\hat{G}, \quad (3.45)$$

doing variation we can find

$$\frac{\delta^2 S}{\delta\hat{G}\delta\hat{G}} \delta\hat{G}^2 = \frac{L}{2} \sum_{l,l's,s'} \int_0^T dt_1 dt_2 c^{ss'} (\delta G_{ll'}^{ss'})^2 \delta(t_{12}), \quad (3.46)$$

$$\frac{\delta^2 S}{\delta\hat{G}\delta\hat{\Sigma}} \delta\hat{\Sigma}\delta\hat{G} = -L \sum_{l,l's,s'} \int dt_1 dt_2 \delta\Sigma_{ll'}^{ss'}(t_1, t_2) \delta G_{ll'}^{ss'}(t_1, t_2). \quad (3.47)$$

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,

$$\frac{\delta^2 S}{\delta\hat{\Sigma}\delta\hat{\Sigma}} \delta\hat{\Sigma}^2 = -\frac{L}{2} \int dt_1 dt_2 G_{ll'}^{ss'}(t_1, t_2) \delta\Sigma_{l'm}^{s'\alpha}(t_1, t_2) G_{mm'}^{\alpha\alpha'}(t_1, t_2) \delta\Sigma_{m'l}^{\alpha's}(t_1, t_2), \quad (3.48)$$

where the field $\hat{G}$ is a solution of the saddle point equation, and we suppose the summation over the repeating indexes. Integrating out the fluctuations, we obtain:

$$\int D\delta\hat{G} D\delta\hat{\Sigma} e^{\frac{1}{2}S^{(2)}(\delta\hat{G}, \delta\hat{\Sigma})} \sim L^{-R_m}. \quad (3.49)$$

Here $R_m$ is a number of massive modes. We also need to calculate the normalisation term, which is a standard Gaussian integral

$$\tilde{F}^{(k)}(0) = \int D\hat{G} D\hat{\Sigma} e^{\left[ \frac{L}{4} \sum_{l,l's,s'} \int_0^T dt c_{ss'} (G_{ll'}^{ss'})^2 - \frac{L}{2} \int dt_1 dt_2 \Sigma_{ll'}^{ss'}(t_1, t_2) G_{ll'}^{ss'}(t_1, t_2) \right]} \sim \\ \sim L^{-\frac{d}{2}} \int D\delta\Sigma e^{\frac{L}{4} \sum_{l,l's,s'} \int_0^T dt c^{ss'} (\Sigma_{ll'}^{ss'})^2(t)} \sim L^{-d}, \quad (3.50)$$

and $d$ is a number of degrees of freedom. Putting it all together, we obtain for $q = 2$:

$$F_{q=2}^{(k)} = \frac{\tilde{F}^{(k)}(T)}{\tilde{F}^{(k)}(0)} \sim \frac{L^{d-R_m}}{\mathcal{N}^{2k}} \sim \frac{L^{R_{ml}}}{\mathcal{N}^{2k}}, \quad (3.51)$$
where $R_{ml} = \frac{k(k-1)}{2}$ is a number of massless modes, taking the logarithm of the Frame potential we can get:

$$ \lim_{L \to \infty} \left( \frac{1}{L} \log F_{q=2}^{(k)} \right) \sim -k \log(2) + \frac{k(k-1)}{2} \frac{1}{L} \log(L), \qquad (3.52) $$

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

$$ F_{q>2}^{(k)} \sim \frac{k!}{\mathcal{N}^{2k}}, \qquad (3.53) $$

together with,

$$ \lim_{L \to \infty} \left( \frac{1}{L} \log F_{q>2}^{(k)} \right) \sim -k \log(2), \qquad (3.54) $$

which indeed coincides with the first correction in system size for gHaar distributed Majorana fermions, see section A.1.

## 4 Dirac SYK model and Kelsdysh saddle points

We shall now consider the Brownian evolution given by the complex SYK Hamiltonian

$$ H(t) = \sum_{\substack{1 \le i_1 \le i_2 \dots \le i_{q/2} \le L \\ 1 \le i_{q/2+1} \le \dots \le i_q \le L}} h_{i_1 \dots i_q} \hat{c}_i^\dagger \dots \hat{c}_{i_{q/2}}^\dagger \hat{c}_{i_{q/2+1}} \dots \hat{c}_{i_q}. \qquad (4.1) $$

Here, $\hat{c}^\dagger$ and $\hat{c}$ are Dirac Fermions, and standard anti-commutation relations are applied $\{\hat{c}_k, \hat{c}^\dagger\} = \delta_{kl}$. As in the previous section, $h_{i_1 \dots i_q}$ are normally distributed random variables, but this time, they are complex and have the variance:

$$ \langle h_{i_1 \dots i_q}(t_1) h_{i_1 \dots i_q}^*(t_2) \rangle = \frac{2^q (q/2)!^2}{q^2 L^{q-1}} \delta(t_1 - t_2). \qquad (4.2) $$

We shall then repeat the analysis of the previous section. The model now possesses an extra symmetry compared to the Majorana case, which is the global U(1) charge conservation

$$ \hat{Q} = \sum_{i=1}^{L} \hat{c}_{i}^{\dagger} \hat{c}_{i}. \qquad (4.3) $$

We shall then show how this impacts the saddle point solutions at late time, both for $q > 2$ and $q = 2$.
## 4.1 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(\epsilon) = e^{-i\epsilon 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 = \frac{\langle Q \rangle}{L}$. The Frame potential can be expressed as

$$F^{(k)} = \frac{1}{\mathcal{N}^{2k}} \underbrace{\mathrm{Tr} \left[ \mathcal{T} \exp \left( -i \int_0^T dt H(t) \right) \right] \mathrm{Tr} \left[ \mathcal{T} \exp \left( i \int_0^T dt H(t) \right) \right] \dots}_{k \text{ times}} \quad (4.4)$$

Since each of the Hamiltonian conserves the total number of particles $\hat{Q} = \sum_i \hat{c}_i^\dagger \hat{c}_i$, we can split each trace into smaller traces over each charge sector, namely

$$\mathrm{Tr} \left[ \mathcal{T} \exp \left( \pm i \int_0^T dt H(t) \right) \right] = \sum_{Q^{\pm}} \mathrm{Tr} \left[ \mathcal{T} \exp \left( \pm i \int_0^T dt H(t) \right) P_{Q^{\pm}} \right], \quad (4.5)$$

with

$$P_{Q^{\pm}} = \int d\mu^{\pm} e^{\pm iL\mu^{\pm}} \int_{0}^{T} dt (\frac{1}{L} \sum_{i} \hat{c}_{i}^{\dagger} \hat{c}_{i-n^{\pm}}). \quad (4.6)$$

Namely, for each Keldysh contour and each replica we can introduce a chemical potential $\mu_l^s$, and we denote by $\mu$ the whole set of them. Writing down the path integral representation for each trace and by averaging over Hamiltonian realisations, we obtain

$$F^{(k)} = \frac{1}{\mathcal{N}^{2k}} \int d\eta(\mathbf{h}) \sum_{Q} \int d\mu \int \mathcal{D}\psi \mathcal{D}\bar{\psi} e^{S(\bar{\psi}, \psi)} \prod_{s,l} P_s(\mu_l^s), \quad (4.7)$$

with the action $S(\bar{\psi}, \psi)$ already introduced in eq. (2.19) and with the projectors now written as

$$P_s(\mu_l^s) = e^{i s T L \mu_l^s} (\bar{\psi}^{s,l}(0^+) \psi^{s',l'}(0) - n_l^s + \delta_{s,+} \delta_{l,l'}). \quad (4.8)$$

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:

$$S = L \log(\det(\hat{\sigma}_z \partial_t - \hat{\Sigma} + i\hat{\sigma}_z \hat{\mu})) + \frac{L}{2q^2} \sum_{l,l',s,s'} \int dt_1 dt_2 \delta(t_{12}) c_{ss'} (2G_{ll'}^{ss'})^{q/2} (-2G_{l'l}^{s's})^{q/2} \\ - L \sum_{l,l',s,s'} \int dt_1 dt_2 \Sigma_{ll'}^{ss'}(t_1,t_2) G_{l'l}^{ss'}(t_1,t_2) + iLT \mu_l^+(1-n_l^+) - iLT \mu_l^-(n_l^-). \quad (4.9)$$

Here $\hat{\mu} = \mu_l^s \delta_{ss'} \delta_{l'l}$ 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 $G_{ll'}^{ss'}(t) = \frac{1}{L} \langle \sum_i \psi_{i,l}^s(t) \bar{\psi}_{i,l'}^{s'}(0) \rangle$. We are interested in the equations of motion. Repeating the procedure from the section with Majorana fermions, one can find:

$$-(\hat{G}^{-1})^T(t_1, t_2) = \hat{\sigma}_z \partial_t - \hat{\Sigma}(t_1, t_2) + i\sigma_z \hat{\mu} \delta(t_{12}), \quad (4.10)$$

$$\Sigma^{ss'}(t_1, t_2) = \frac{c_{ss'}}{q} (2G^{ss'})^{\frac{q}{2}-1} (t_1, t_2) (-2G^{s's})^{\frac{q}{2}} (t_1, t_2) \delta(t_{12}), \quad (4.11)$$

$$G_{l,l'}^{++}(0^+) = \delta_{l,l'}(1-n_l^+), \quad G_{l,l'}^{--}(0^+) = -\delta_{l,l'}n_l^-. \quad (4.12)$$
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 $\hat{G}$ and $\hat{\Sigma}$ were real-valued. Now, they are allowed to be complex. In the case of single replica $k=1$, in the limit $t \to 0^\pm$ the Green function $G^{ss'}$ was found already in [44, 45] 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

$$F^{(k)}(T \to 0) = e^{-\frac{LTk}{q^2}}. \qquad (4.13)$$

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 $\hat{G}$, 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 B

$$\hat{G}_{\text{diag}} = \begin{pmatrix} f(t_{12})(\Theta(t_{12})(\hat{1} - \hat{n})) - \Theta(-t_{12})\hat{n} & -f(t_{12})\hat{n} \\ f(t_{12})(\hat{1} - \hat{n}) & f(t_{12})(\Theta(-t_{12})(\hat{1} - \hat{n}) - \Theta(t_{12})\hat{n}) \end{pmatrix}, \quad (4.14)$$

and

$$\hat{\Sigma}_{\text{diag}} = \frac{2^{q-1}}{q} \begin{pmatrix} \frac{1}{2}\Gamma(1-2\hat{n}) & -\Gamma(1-\hat{n}) \\ \hat{n}\Gamma & \frac{1}{2}\Gamma(1-2\hat{n}) \end{pmatrix} \delta(t_{12}), \quad (4.15)$$

with $f(t_{12}) = e^{-\frac{\Gamma|t_{12}|}{2}-i\hat{\mu}t_{12}}$ and $\hat{n} = \text{diag}(n_1, n_2, \dots, n_k)$ is a diagonal matrix in replica space. This matrix consists of the fermion densities from different replicas and $\Gamma = (\hat{n}(1-\hat{n}))^{\frac{q}{2}-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 $\hat{n}^- = \hat{\tau}\hat{n}^+\hat{\tau}^T$, due to the last two equations in eq. (4.10). This gives the solution:

$$\hat{G} = \begin{pmatrix} f(t_{12})(\Theta(t_{12})(\hat{1} - \hat{n})) - \Theta(-t_{12})\hat{n} & -f(t_{12})\hat{n}\hat{\tau}^T \\ \hat{\tau}f(t_{12})(\hat{1} - \hat{n}) & \hat{\tau}f(t_{12})(\Theta(-t_{12})(\hat{1} - \hat{n})) - \Theta(t_{12})\hat{n}\hat{\tau}^T \end{pmatrix}. \quad (4.16)$$

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:

$$e^{S_1} = (\det[\hat{\sigma}_z \partial_t - \hat{\Sigma}_d + i\sigma_z \hat{\mu}\delta(t_{12})])^L = \text{Tr}(e^{-T\hat{H}})^L, \quad (4.17)$$

where

$$H = (c_+^\dagger \quad c_-^\dagger) \begin{pmatrix} \hat{1} & 0 \\ 0 & \hat{i} \end{pmatrix} H_0 \begin{pmatrix} \hat{1} & 0 \\ 0 & \hat{i} \end{pmatrix} \begin{pmatrix} c_+ \\ c_- \end{pmatrix}, \quad (4.18)$$

$$H_0 = \begin{pmatrix} \frac{2^{q-1}}{2q} \Gamma(1-2\hat{n}) & -i\hat{\mu} & -\frac{2^{q-1}}{q} \Gamma(1-\hat{n}) \\ \frac{2^{q-1}}{q} \hat{n}\Gamma & \frac{2^{q-1}}{2q} \Gamma(1-2\hat{n}) & +i\hat{\mu} \end{pmatrix}. \quad (4.19)$$
Then we use the determinant formula

$$ \mathrm{Tr}(e^{-TH}) = \det \left( 1 + \exp \left( -T \begin{pmatrix} \hat{1} & 0 \\ 0 & \hat{i} \end{pmatrix} H_0 \begin{pmatrix} \hat{1} & 0 \\ 0 & \hat{i} \end{pmatrix} \right) \right). \quad (4.20) $$

Diagonalisation of the matrix exponent gives eigenvalues $e^{\pm \frac{T\Gamma_l}{2q} - iT\mu_l} = e^{\pm \frac{T(n_l(1-n_l))q}{2q} - iT\mu_l}$, where $\Gamma_l$ is a matrix element of the diagonal matrix $\Gamma$. Which gives the answer for the determinant calculation:

$$ \det(\sigma_z \partial_t - \Sigma_d + i\sigma_z \hat{\mu}\delta(t_{12}))^L \rightarrow e^{\frac{2^{q-1}LT}{2q} \sum_{l=1}^k \Gamma_l} - iLT \sum_l \mu_l. \quad (4.21) $$

Now let us move to the second term in the action:

$$ \begin{aligned} S_2 &= \frac{L}{2q^2} \sum_{l,l',s,s'} \int dt_1 dt_2 \delta(t_{12}) c_{ss'} (2G_{ll'}^{ss'})^{q/2} (-2G_{l'l}^{s's})^{q/2} = \\ &= \frac{L}{q^2} \sum_l \left( n_l^{\frac{q}{2}} (1-n_l)^{\frac{q}{2}} - \int dt \left( \Theta^q(t) (1-n_l)^{\frac{q}{2}} n_l^{\frac{q}{2}} + \Theta(-t)^q n_l^{\frac{q}{2}} (1-n_l)^{\frac{q}{2}} \right) \delta(t) \right) = 0, \end{aligned} \quad (4.22) $$

notice that both side limit for the sum of Heaviside functions is not well-defined, therefore we assume:

$$ \begin{aligned} \delta(t)(\Theta(t)a + \Theta(-t)b) &= \frac{1}{2} \left( \lim_{t \to 0^+} (\Theta(t)a + \Theta(-t)b) + \lim_{t \to 0^-} (\Theta(t)a + \Theta(-t)b) \right) \\ &= \frac{1}{2}(a+b). \end{aligned} \quad (4.23) $$

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,

$$ \begin{aligned} S_3 &= -L \sum_{l,l',s,s'} \int dt_1 dt_2 \Sigma_{ll'}^{ss'}(t_1, t_2) G_{ll'}^{ss'}(t_1, t_2) = \\ &= -L \int dt_1 dt_2 \left( \Sigma_{ll'}^{++} G_{ll'}^{++} + \Sigma_{ll'}^{--} G_{ll'}^{--} + \Sigma_{ll'}^{+-} G_{ll'}^{+-} + \Sigma_{ll'}^{-+} G_{ll'}^{-+} \right) = \\ &= -\frac{2^{q-1}}{2q} \sum_l (1-2n_l)^2 \Gamma_l - \frac{2^{q-1}}{q} \sum_l \Gamma_l (1-n_l)n_l = -\frac{2^{q-1}LT \sum_{l=1}^k \Gamma_l}{2q}. \end{aligned} \quad (4.24) $$

We find therefore the final expression for the Frame potential

$$ F^{(k)} = \sum_{\hat{\tau}} \int d\boldsymbol{\mu}^{+} \sum_{\boldsymbol{Q}^{+}} \frac{e^{S_1 + S_2 + S_3}}{\mathcal{N}^{2k}} = \frac{k! L^k}{\mathcal{N}^{2k}}, \quad (4.25) $$

where we used $\sum_{\boldsymbol{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$:

$$ 
\begin{aligned}
-(\hat{G}^{-1})^T(t_1, t_2) &= \hat{\sigma}_z \partial_t - \hat{\Sigma}(t_1, t_2) + i\sigma_z \hat{\mu}\delta(t_{12}), \\
\Sigma^{ss'}(t_1, t_2) &= \frac{c_{ss'}}{2}(-2G^{s's})(t_1, t_2)\delta(t_{12}), \\
G_l^{++}(0^+) &= 1 - n_l^+, \quad G_l^{--}(0^+) = -n_l^-.
\end{aligned} 
\quad (4.26) $$

Similar to the case with $q > 2$ we start with a replica diagonal solution with $n_l^+ = n_l^-$ different for each replica:

$$ \hat{G}_{\text{diag}} = \begin{pmatrix} f(t_{12})(\Theta(t_{12})(\hat{1} - \hat{n}) - \Theta(-t_{12})\hat{n}) & -f(t_{12})\hat{n} \\ f(t_{12})(\hat{1} - \hat{n}) & f(t_{12})(\Theta(-t_{12})(\hat{1} - \hat{n}) - \Theta(t_{12})\hat{n}) \end{pmatrix}, \quad (4.27) $$

where $f(t_{12}) = e^{-\frac{|t_{12}|}{2}} - i\hat{\mu}t_{12}$, since $\Gamma = \hat{1}$ in the case $q = 2$ and

$$ \hat{\Sigma}_{\text{diag}} = \begin{pmatrix} \frac{1}{2}(\hat{1} - 2\hat{n}) & -(\hat{1} - \hat{n}) \\ \hat{n} & \frac{1}{2}(\hat{1} - 2\hat{n}) \end{pmatrix} \delta(t_{12}). \quad (4.28) $$

Since our equations of motion have rotation symmetry in replica space we can consider the solution with an arbitrary unitary matrix $\hat{u}$:

$$ \hat{G} = \begin{pmatrix} f(t_{12})(\Theta(t_{12})(\hat{1} - \hat{n}) - \Theta(-t_{12})\hat{n}) & -f(t_{12})\hat{n}\hat{u}^{\dagger} \\ \hat{u}f(t_{12})(\hat{1} - \hat{n}) & \hat{u}f(t_{12})(\Theta(-t_{12})(\hat{1} - \hat{n}) - \Theta(t_{12})\hat{n})\hat{u}^{\dagger} \end{pmatrix}. \quad (4.29) $$

However notice that this solution implies boundary conditions $\hat{n}^- = \hat{u}\hat{n}^+\hat{u}^\dagger$, 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 $\hat{n} = n\delta_{ij} = n\hat{1}$. Notice that for this form of the solution, we have a generic rotated solution, with an arbitrary unitary matrix $\hat{u}$. It is clear that boundary conditions are simplified in this case, as $n^{-}\hat{1} = \hat{u}(n^{+}\hat{1})\hat{u}^{\dagger} = n^{+}\hat{1}$, and terms of this form are presented in the sum over the charge content. Therefore, the solution takes the form:

$$ \hat{G} = \begin{pmatrix} f(t_{12}) (\Theta(t_{12})(\hat{1} - n\hat{1}) - \Theta(-t_{12})n\hat{1}) & -f(t_{12})\hat{n}\hat{1} \\ f(t_{12})(\hat{1} - n\hat{1})\hat{u}^{\dagger} & f(t_{12})(\Theta(-t_{12})(\hat{1} - n\hat{1}) - \Theta(t_{12})n\hat{1}) \end{pmatrix}, \quad (4.30) $$

where $f(t_{12}) = e^{-\frac{|t_{12}|}{2}} - i\hat{\mu}\hat{1}t_{12}$, with $\hat{\mu}$ being just a number here. Here $\hat{u}$ is a unitary matrix in $k \times k$ space. Notice also that not all possible matrices $\hat{u}$ respect the parity of the initial state. Here we need to assume that $\det(\hat{u}) = 1$, in order to take into account the parity. Therefore, $\hat{u}$ 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 $\hat{n} \to n\hat{1}$.
### 4.2.1 Symmetries, and fluctuations

Let us start with the consideration of the Goldstone theorem for the case $q = 2$. The vector fields $\psi$ and $\bar{\psi}$ 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(\hat{\sigma}_z \partial - \hat{\Sigma} + i\hat{\sigma}_z \hat{\mu}))$ and other terms with chemical potential. Notice that in the case $\hat{n}^+ = \hat{n}^- = n\hat{1}$ and $\hat{\mu}^+ = \hat{\mu}^- = \mu\hat{1}$ the action again has rotation symmetries in the replica space $U_1 = e^{i\hat{A}_1\hat{\sigma}_z}$, $U_2 = e^{i\hat{A}_2\hat{1}}$, where $\hat{A}_i$ is a Hermitian matrix of size $k$. By $\hat{A}\hat{\sigma}_z$, we again mean the tensor product. Notice that matrices $U_0 = e^{i\hat{A}}$ form the $U(k)$ group. Therefore, the group of symmetries of the action is $B = U(k) \times U(k)$.

Remember that the saddle point solution has the form eq. (4.30) where $\hat{u}$ 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:

$$R_m = d - R_{ml} = 4k^2 - k^2 + 1 = 3k^2. \qquad (4.31)$$

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, \Sigma_0) + S^{(2)}(\delta\hat{G}, \delta\hat{\Sigma})$ 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

$$\tilde{F}^{(k)}(T) \sim \sum_{Q=1}^{L} \frac{1}{\mathcal{N}^{2k}} L^{-R_m}. \qquad (4.32)$$

The normalisation term can be calculated the same way as eq. (3.50). Which gives the scaling with the system size:

$$\tilde{F}^{(k)}(0) \sim L^{-d}. \qquad (4.33)$$

Therefore in the final expression for the Frame potential, we get:

$$F_{q=2}^{(k)} = \frac{\tilde{F}^{(k)}(T)}{\tilde{F}^{(k)}(0)} \sim \frac{L}{\mathcal{N}^{2k}} L^{R_{ml}}, \qquad (4.34)$$

and the logarithm is

$$\lim_{L \to \infty} \left( \frac{1}{L} \log F_{q=2}^{(k)} \right) \sim -2k \log(2) + k^2 \frac{1}{L} \log(L) + O\left(\frac{1}{L}\right), \qquad (4.35)$$

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:

$$\lim_{L \to \infty} \left( \frac{1}{L} \log F_{q>2}^{(k)} \right) \sim -2k \log(2) + O\left(\frac{1}{L}\right). \qquad (4.36)$$

Notice that this result also coincides with the gHaar integrated result in eq. (A.11).
## 5 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_{\text{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_{\text{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 [46–49]. 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.

## Acknowledgments

We thank Lorenzo Piroli, Guido Giachetti, and Adam Nahum for useful discussions. This work has been partially funded by the ERC Starting Grant 101042293 (HEPIQ) (J.D.N. and A.T.) and by the ANR JCJC grant ANR-21-CE47-0003 (TamEnt) (A.D.L.).
## A Haar-averaged Frame potential

As we discussed during our work, it is expected for the Frame potential to converge to its Haar value during Brownian evolution. In this appendix, we calculate the Haar values for the Frame potential. First, we consider the case when the evolution matrix $U$ is selected from the corresponding group using the Haar measure. The groups we consider here are the 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\dagger 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^{\frac{1}{2}h_{ij}\chi_i\chi_j}$, where $\chi_j$ are Majorana fields, then we can express $\text{Tr}(U) = \text{Pf}(1 + e^{h_{ij}})$ (as for example shown in [50]) where $e^{h_{ij}} = \hat{u}$ is an orthogonal matrix with the size $L \times L$ ($L$ is even). We are interested in the case when matrices $\hat{u}$ are random with respect to the Haar measure. Notice that we are averaging over matrices with $\text{det}(\hat{u}) = 1$, as $\hat{u} = e^{\frac{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 $\hat{u}$ are $\lambda_1 = e^{ib_1}$, $\lambda_1^* = e^{-ib_1}$, ... and their product is 1). This gives us the following averaging

$$ 
\begin{aligned}
F_{\text{gHaar}}^k &= \frac{1}{\mathcal{N}^{2k}} \int |\det(1 + \hat{u})|^k d\mu(\hat{u}) = \\
&= \frac{1}{\mathcal{N}^{2k}} \int \left| \prod_{j=1}^{L} (1 + e^{i\theta_j}) \right|^k \prod_{i<j} |2(\cos\theta_i - \cos\theta_j)|^2 \frac{d\theta_1 \dots d\theta_{L/2}}{\left(\frac{L}{2}\right)! (2\pi)^{\frac{L}{2}}} = \\
&= \frac{1}{\mathcal{N}^{2k}} \int \left| \prod_{j=1}^{L/2} (2 + 2\cos(\theta_j)) \right|^k \prod_{i<j} |2(\cos\theta_i - \cos\theta_j)|^2 \frac{d\theta_1 \dots d\theta_{L/2}}{\left(\frac{L}{2}\right)! (2\pi)^{\frac{L}{2}}}. \quad (\text{A.1})
\end{aligned}
$$

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 [51]). Doing change of variables $\cos(\theta) = t$ and calculating the Jacobian of this transformation we can get

$$ F_{\text{gHaar}}^k = \frac{1}{\mathcal{N}^{2k}} \frac{1}{F^0} \int_{-1}^{1} dt_1 \dots \int_{-1}^{1} dt_{L/2} \prod_{j=1}^{L/2} |2(1+t_j)|^k \prod_{i<j} |t_i - t_j|^2 \prod_{k=1}^{L/2} (1-t_k^2)^{-\frac{1}{2}}. \quad (\text{A.2}) $$

Now we need to calculate these integrals. For simplicity, we can introduce the function to name these integrals:

$$ F = 2^{\frac{Lk}{2}} \int_{-1}^{1} dt_1 \dots \int_{-1}^{1} dt_{L/2} \prod_{j=1}^{L/2} |1+t_j|^k \prod_{i<j} |t_i - t_j|^2 \prod_{k=1}^{L/2} (1-t_k^2)^{-\frac{1}{2}}. \quad (\text{A.3}) $$
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 [51])

$$F = 2^{\frac{Lk}{2}} 2^{\frac{L}{2} + \frac{L(k-1)}{2}} + \frac{L}{2}(\frac{L}{2}-1) S_{\frac{L}{2}}\left(k-\frac{1}{2}, -\frac{1}{2}, 1\right) = 2^{Lk} 2^{\frac{L}{2}(\frac{L}{2}-1)} S_{\frac{L}{2}}\left(k-\frac{1}{2}, -\frac{1}{2}, 1\right). \quad (A.4)$$

The Selberg integral has a well-known expression in terms of Gamma functions, which is:

$$S_L(\lambda_1, \lambda_2, \lambda) = \int_0^1 dt_1 \dots \int_0^1 dt_L \prod_{l=1}^L t_l^{\lambda_1} (1-t_l)^{\lambda_2} \prod_{k<j} |t_k - t_j|^{2\lambda} = \\ = \prod_{j=0}^{L-1} \frac{\Gamma(\lambda_1 + 1 + j\lambda)\Gamma(\lambda_2 + 1 + j\lambda)\Gamma(1 + (j+1)\lambda)}{\Gamma(\lambda_1 + \lambda_2 + 2 + (L+j-1)\lambda)\Gamma(\lambda+1)}. \quad (A.5)$$

Inserting this solution into the eq. (A.2) we can have the final expression

$$F_{\text{gHaar}}^k = \frac{2^{Lk} 2^{\frac{L}{2}(\frac{L}{2}-1)} S_{\frac{L}{2}}\left(k-\frac{1}{2}, -\frac{1}{2}, 1\right)}{2^{\frac{L}{2}(\frac{L}{2}-1)} S_{\frac{L}{2}}\left(-\frac{1}{2}, -\frac{1}{2}, 1\right)} = \frac{S_{\frac{L}{2}}\left(k-\frac{1}{2}, -\frac{1}{2}, 1\right)}{S_{\frac{L}{2}}\left(-\frac{1}{2}, -\frac{1}{2}, 1\right)} = \\ = \prod_{j=0}^{\frac{L}{2}-1} \frac{\Gamma\left(k+\frac{1}{2}+j\right)}{\Gamma\left(\frac{1}{2}+j\right)} \frac{\Gamma\left(\frac{L}{2}+j\right)}{\Gamma\left(k+\frac{N}{2}+j\right)} \sim \frac{c_k L^{k(k-1)/2}}{\mathcal{N}^{2k}} + \dots, \quad (A.6)$$

where $c_k = \frac{2}{2^{\frac{(k-1)(k-2)}{2}}} \prod_{i=0}^{k-1} \frac{1}{(2i-1)!!}$ and we define $(-1)!! = 1$. Therefore, for $k=1$ we have

$$F_{\text{gHaar}}^{1} = \frac{2}{2^{L}} = F_{\text{Haar}}^{1}, \quad (A.7)$$

which is equal to the Frame potential in the case when the evolution matrix was taken from the Haar ensemble of Orthogonal matrices. We are interested in the large $L$ limit expansion of the Frame potential for arbitrary $k$, therefore we can use a “Stirling-like” formula for the Barnes functions, which gives:

$$\lim_{L \to \infty} \left( \frac{1}{L} \log \left( F_{\text{gHaar}}^k \right) \right) \sim -k \log(2) + \frac{k(k-1)}{2} \frac{1}{L} \log(L). \quad (A.8)$$

## A.2 Haar average of frame potential for Dirac fermions

Now we focus on the Dirac fermions case. In this case, the large time Haar distribution for the evolution operator is split into different charge blocks, where each block is a unitary random matrix. The Frame potential then takes the form

$$F_{\text{Haar}}^{(k)} = \frac{1}{\mathcal{N}^{2k}} \left[ \prod_Q \int_{\text{Haar}} dU_Q \right] \left| \sum_Q \text{Tr} [U_Q] \right|^{2k}. \quad (A.9)$$

Given that only paired traces give a finite result, we obtain $k!$ pairs whose average is always equal to 1, as

$$\int_{\text{Haar}} dU |\text{Tr}[U]|^2 = 1, \quad (A.10)$$

therefore using $\sum_Q = L$, we obtain

$$\frac{1}{\mathcal{N}^{2k}} \int_{\text{Haar}} dU |\text{Tr}[U]|^{2k} = \frac{L^k}{\mathcal{N}^{2k}(k)!}. \quad (A.11)$$
### A.3 Gaussian-Haar averages for Dirac fermions

In the case of Gaussian evolution, we consider a unitary Gaussian operator $U$, i.e. $U = e^{i\hat{h}_{ij}}c_i^\dagger c_j$, where $c_j$ are Dirac fermions. Then we can express $\text{Tr}(U) = \text{det}(1 + \hat{u})$ where $e^{i\hat{h}_{ij}} = \hat{u}$ is a unitary matrix with the size $L \times L$ (with $L$ even). The matrices $\hat{u}$ are again random with respect to the Haar measure. This gives us the following average,

$$ 
\begin{aligned}
F_{\text{gHaar}}^{(k)} &= \frac{1}{\mathcal{N}^{2k}} \int |\det(1 + \hat{u})|^{2k} d\mu(\hat{u}) = \\
&= \frac{1}{\mathcal{N}^{2k}} \int \prod_{j=1}^{L} |(1 + e^{i\theta_j})|^{2k} \prod_{i<j} |e^{i\theta_i} - e^{i\theta_j}|^2 \frac{d\theta_1 \dots d\theta_N}{N!(2\pi)^N}.
\end{aligned}
\quad (A.12) $$

The measure in this expression comes from the fact that we integrate other unitary matrices (see chapter 2.6 in [51]). This integral is a Morris integral which can be expressed in terms of Gamma functions:

$$ 
\begin{aligned}
M_L(a, b, \lambda) &= \int_{-1/2}^{1/2} \prod_{j=1}^{L} e^{i\theta_i(a-b)/2} |(1+e^{i\theta_j})|^{a+b} \prod_{i<j} |e^{i\theta_i} - e^{i\theta_j}|^{2\lambda} d\theta_1 \dots d\theta_N = \\
&= \prod_{j=0}^{L-1} \frac{\Gamma(\lambda j + a + b + 1)\Gamma(\lambda j + 1 + \lambda)}{\Gamma(\lambda j + a + 1)\Gamma(\lambda j + b + 1)\Gamma(1 + \lambda)}.
\end{aligned}
\quad (A.13) $$

Inserting this into the original expression for the Frame potential we can get:

$$ F_{\text{gHaar}}^{k} = \frac{1}{\mathcal{N}^{2k}} \frac{M_L(k, k, 1)}{M_L(0, 0, 1)} = \frac{1}{\mathcal{N}^{2k}} \prod_{j=0}^{L-1} \frac{\Gamma(j + 2k + 1)\Gamma(j + 1)}{\Gamma(j + k + 1)^2}. \quad (A.14) $$

Again we are interested in the large $L$ limit expansion:

$$ F_{\text{Haar}}^{(k)} = \frac{\tilde{c}_k L^{k^2}}{\mathcal{N}^{2k}} + \dots \quad (A.15) $$

where we neglected subleading corrections in $L$, and $\tilde{c}_k = \frac{\text{sf}^2(k-1)}{\text{sf}(2k-1)}$ with $\text{sf}(n) = 1!2!\dots n!$ a super factorial. Comparing this to the full Haar average, we again find

$$ F_{\text{gHaar}}^{k} \geq F_{\text{Haar}}^{k}, \quad (A.16) $$

therefore the full Haar average is again smaller than the fermionic one.

### B Dirac SYK saddle point with fixed charge density

Let's start with replica-diagonal solution $Q_l^+ = Q_l^-$, we use the following ansatz

$$ \hat{G}_d = \frac{1}{2} \begin{pmatrix} f(\Theta(t_{12})(\hat{1} - \hat{n}) - \Theta(-t_{12})\hat{n}) & -f\hat{n} \\ f(\hat{1} - \hat{n}) & f(\Theta(-t_{12})(\hat{1} - \hat{n}) - \Theta(t_{12})\hat{n}) \end{pmatrix}, \quad (B.1) $$

where $f_{ll'}(t)$ is a time-dependent matrix in replica space with initial condition $f_{ll'}(0) = f(0)\delta_{ll'}$. First, we consider the equation

$$ \Sigma_{ll'}^{ss'}(t) = \frac{c_{ss'}}{q} (2G_{ll'}^{ss'})^{\frac{q}{2}-1} (t) (-2G_{ll'}^{s's})^{\frac{q}{2}} (-t) \delta(t), \quad (B.2) $$
which gives

$$ \hat{\Sigma}_d = \frac{f^{q-1}(0)}{q} \begin{pmatrix} \frac{1}{2}\Gamma(\hat{1} - 2\hat{n}) & -\Gamma(\hat{1} - \hat{n}) \\ \hat{n}\Gamma & \frac{1}{2}\Gamma(\hat{1} - 2\hat{n}) \end{pmatrix} \delta(t_{12}), \quad (B.3) $$

where $\hat{n} = \text{diag}(n_1, n_2, \dots, n_k)$ is a diagonal matrix in replica space, which consists of the fermion densities from different replicas, and $\Gamma = (\hat{n}(\hat{1} - \hat{n}))^{\frac{q}{2}-1}$ is also a diagonal matrix. Using the equation $-(\hat{G}^{-1})_{ll'}^{s's}(t) = \hat{\sigma}_z \partial_t - \hat{\Sigma}_{ll'}^{ss'}(t) + i\sigma_z \mu_l^s$ and Fourier transformation, we can find

$$ \hat{G}_d^{-1} = \begin{pmatrix} i\omega + \frac{1}{2}\tilde{\Gamma}(\hat{1} - 2\hat{n}) - i\hat{\mu}^{+} & \tilde{\Gamma}\hat{n} \\ -(\hat{1} - \hat{n})\tilde{\Gamma} & -i\omega + \frac{1}{2}\tilde{\Gamma}(\hat{1} - 2\hat{n}) + i\hat{\mu}^{-} \end{pmatrix} \quad (B.4) $$

where $\tilde{\Gamma} = \frac{f^{q-1}(0)}{q}\Gamma$. Here we can use the formula for the block matrix and assume $\mu_l^+ = \mu_l^- = \mu_l$ (since $Q_l^+ = Q_l^-$), therefore we find

$$ \begin{pmatrix} A & B \\ C & D \end{pmatrix}^{-1} = \begin{pmatrix} A^{-1} + A^{-1}B(D-CA^{-1}B)^{-1}CA^{-1} & -A^{-1}B(D-CA^{-1}B)^{-1} \\ -(D-CA^{-1}B)^{-1}CA^{-1} & (D-CA^{-1}B)^{-1} \end{pmatrix}, \quad (B.5) $$

which gives

$$ \hat{G}_d = \frac{1}{(\omega - \hat{\mu})^2 + \tilde{\Gamma}^2 \frac{1}{4}} \begin{pmatrix} -i\omega + \tilde{\Gamma}/2(1 - 2\hat{n}) + i\hat{\mu} & -\tilde{\Gamma}\hat{n} \\ (1 - \hat{n})\tilde{\Gamma} & i\omega + \frac{1}{2}\tilde{\Gamma}(1 - 2\hat{n}) - i\hat{\mu} \end{pmatrix}. \quad (B.6) $$

The inverse Fourier transform of each element of the matrix gives

$$ \hat{G}_d = \begin{pmatrix} e^{-\frac{\Gamma|t_{12}|}{2}} - i\hat{\mu}t_{12}(\Theta(t_{12})(\hat{1} - \hat{n}) - \Theta(-t_{12})\hat{n}) & -e^{-\frac{\Gamma|t_{12}|}{2}} - i\hat{\mu}t_{12}\hat{n} \\ e^{-\frac{\Gamma|t_{12}|}{2}} - i\hat{\mu}t_{12}(\hat{1} - \hat{n}) & e^{-\frac{\Gamma|t_{12}|}{2}} - i\hat{\mu}t_{12}(\Theta(-t_{12})(\hat{1} - \hat{n}) - \Theta(t_{12})\hat{n}) \end{pmatrix}. \quad (B.7) $$

**Open Access.** This article is distributed under the terms of the Creative Commons Attribution License (CC-BY4.0), which permits any use, distribution and reproduction in any medium, provided the original author(s) and source are credited.

## References

[1] P. Hayden and J. Preskill, *Black holes as mirrors: Quantum information in random subsystems*, JHEP 09 (2007) 120 [arXiv:0708.4025] [INSPIRE].

[2] S.H. Shenker and D. Stanford, *Black holes and the butterfly effect*, JHEP 03 (2014) 067 [arXiv:1306.0622] [INSPIRE].

[3] S.H. Shenker and D. Stanford, *Multiple Shocks*, JHEP 12 (2014) 046 [arXiv:1312.3296] [INSPIRE].

[4] J. Maldacena, S.H. Shenker and D. Stanford, *A bound on chaos*, JHEP 08 (2016) 106 [arXiv:1503.01409] [INSPIRE].

[5] A. Touil and S. Deffner, *Information Scrambling versus Decoherence — Two Competing Sinks for Entropy*, PRX Quantum 2 (2021) 010306 [arXiv:2008.05559] [INSPIRE].

[6] A. Almheiri et al., *Replica Wormholes and the Entropy of Hawking Radiation*, JHEP 05 (2020) 013 [arXiv:1911.12333] [INSPIRE].
[7] S. Xu and B. Swingle, Locality, Quantum Fluctuations, and Scrambling, Phys. Rev. X 9 (2019) 031048 [arXiv:1805.05376] [INSPIRE].

[8] S.-K. Jian, G. Bentsen and B. Swingle, Linear Growth of Circuit Complexity from Brownian Dynamics, arXiv:2206.14205.

[9] D. Gross, K. Audenaert and J. Eisert, Evenly distributed unitaries: On the structure of unitary designs, J. Math. Phys. 48 (2007) 052104 [INSPIRE].

[10] D.A. Roberts and B. Yoshida, Chaos and complexity by design, JHEP 04 (2017) 121 [arXiv:1610.04903] [INSPIRE].

[11] S.N. Hearth, M.O. Flynn, A. Chandran and C.R. Laumann, Unitary k-designs from random number-conserving quantum circuits, arXiv:2306.01035 [INSPIRE].

[12] A. Christopoulos, P. Le Doussal, D. Bernard and A. De Luca, Universal Out-of-Equilibrium Dynamics of 1D Critical Quantum Systems Perturbed by Noise Coupled to Energy, Phys. Rev. X 13 (2023) 011043 [arXiv:2110.15303] [INSPIRE].

[13] M. Srednicki, Chaos and Quantum Thermalization, Phys. Rev. E 50 (1994) 888 [cond-mat/9403051] [INSPIRE].

[14] L. D'Alessio, Y. Kafri, A. Polkovnikov and M. Rigol, From quantum chaos and eigenstate thermalization to statistical mechanics and thermodynamics, Adv. Phys. 65 (2016) 239 [arXiv:1509.06411] [INSPIRE].

[15] S. Pappalardi, L. Foini and J. Kurchan, Eigenstate Thermalization Hypothesis and Free Probability, Phys. Rev. Lett. 129 (2022) 170603 [arXiv:2204.11679] [INSPIRE].

[16] S. Hortikar and M. Srednicki, Random matrix elements and eigenfunctions in chaotic systems, Phys. Rev. E 57 (1998) 7313.

[17] J.M. Deutsch, Quantum statistical mechanics in a closed system, Phys. Rev. A 43 (1991) 2046 [INSPIRE].

[18] M. Lucas, L. Piroli, J. De Nardis and A. De Luca, Generalized deep thermalization for free fermions, Phys. Rev. A 107 (2023) 032215 [arXiv:2207.13628] [INSPIRE].

[19] J.S. Cotler et al., Emergent Quantum State Designs from Individual Many-Body Wave Functions, PRX Quantum 4 (2023) 010311 [arXiv:2103.03536] [INSPIRE].

[20] J. Choi et al., Preparing random states and benchmarking with many-body quantum chaos, Nature 613 (2023) 468 [arXiv:2103.03535] [INSPIRE].

[21] W.W. Ho and S. Choi, Exact emergent quantum state designs from quantum chaotic dynamics, Phys. Rev. Lett. 128 (2022) 060601 [arXiv:2109.07491] [INSPIRE].

[22] H. Shrotriya and W.W. Ho, Nonlocality of Deep Thermalization, arXiv:2305.08437 [INSPIRE].

[23] M. Ippoliti and W.W. Ho, Solvable model of deep thermalization with distinct design times, Quantum 6 (2022) 886 [arXiv:2208.10542] [INSPIRE].

[24] S. Sachdev and J. Ye, Gapless spin fluid ground state in a random, quantum Heisenberg magnet, Phys. Rev. Lett. 70 (1993) 3339 [cond-mat/9212030] [INSPIRE].

[25] P. Saad, S.H. Shenker and D. Stanford, A semiclassical ramp in SYK and in gravity, arXiv:1806.06840.

[26] S.-K. Jian and B. Swingle, Note on entropy dynamics in the Brownian SYK model, JHEP 03 (2021) 042 [arXiv:2011.08158] [INSPIRE].
[27] C. Sünderhauf et al., *Quantum chaos in the Brownian SYK model with large finite N: OTOCs and tripartite information*, JHEP 11 (2019) 038 [arXiv:1908.00775] [INSPIRE].

[28] P. Ribeiro, J. Vidal and R. Mosseri, *Exact spectrum of the Lipkin-Meshkov-Glick model in the thermodynamic limit and finite-size corrections*, Phys. Rev. E 78 (2008) 021106 [INSPIRE].

[29] X. Chen and T. Zhou, *Quantum chaos dynamics in long-range power law interaction systems*, Phys. Rev. B 100 (2019) 064305 [arXiv:1808.09812] [INSPIRE].

[30] D. Bernard, F.H.L. Essler, L. Hruza and M. Medenjak, *Dynamics of fluctuations in quantum simple exclusion processes*, SciPost Phys. 12 (2022) 042 [arXiv:2107.02662] [INSPIRE].

[31] D. Bernard and T. Jin, *Solution to the Quantum Symmetric Simple Exclusion Process: The Continuous Case*, Commun. Math. Phys. 384 (2021) 1141 [arXiv:2006.12222] [INSPIRE].

[32] T. Swann, D. Bernard and A. Nahum, *Spacetime picture for entanglement generation in noisy fermion chains*, arXiv:2302.12212.

[33] D. Bernard and L. Piroli, *Entanglement distribution in the quantum symmetric simple exclusion process*, Phys. Rev. E 104 (2021) 014146 [arXiv:2102.04745] [INSPIRE].

[34] Y. Gu, A. Kitaev, S. Sachdev and G. Tarnopolsky, *Notes on the complex Sachdev-Ye-Kitaev model*, JHEP 02 (2020) 157 [arXiv:1910.14099] [INSPIRE].

[35] Y. Gu, A. Lucas and X.-L. Qi, *Spread of entanglement in a Sachdev-Ye-Kitaev chain*, JHEP 2017 (2017) 1.

[36] G. Penington, S.H. Shenker, D. Stanford and Z. Yang, *Replica wormholes and the black hole interior*, JHEP 03 (2022) 205 [arXiv:1911.11977] [INSPIRE].

[37] R. Shankar, *Quantum Field Theory and Condensed Matter*, Cambridge University Press (2017).

[38] A.V. Lunkin, A.Y. Kitaev and M.V. Feigel'man, *Perturbed Sachdev-Ye-Kitaev Model: A polaron in the Hyperbolic Plane*, Phys. Rev. Lett. 125 (2020) 196602 [arXiv:2006.14535] [INSPIRE].

[39] M. Nakahara, *Geometry, topology and physics*, CRC Press (2003) [DOI:10.1201/9781315275826].

[40] S.-K. Jian and B. Swingle, *Phase transition in von Neumann entanglement entropy from replica symmetry breaking*, JHEP 11 (2023) 221.

[41] A.Y. Kitaev, *Unpaired Majorana fermions in quantum wires*, Phys.-Usp. 44 (2001) 131.

[42] M. Fava et al., *Nonlinear Sigma Models for Monitored Dynamics of Free Fermions*, Phys. Rev. X 13 (2023) 041045 [arXiv:2302.12820] [INSPIRE].

[43] A. Baker, *Matrix Groups An Introduction to Lie Group Theory*, Am. Math. Mon. 110 (2003) 446.

[44] P. Zhang, *Information scrambling and entanglement dynamics of complex Brownian Sachdev-Ye-Kitaev models*, JHEP 04 (2023) 105 [arXiv:2301.03189] [INSPIRE].

[45] X. Chen, Y. Gu and A. Lucas, *Many-body quantum dynamics slows down at low density*, SciPost Phys. 9 (2020) 071 [arXiv:2007.10352] [INSPIRE].

[46] S.-K. Jian and B. Swingle, *Phase transition in von Neumann entanglement entropy from replica symmetry breaking*, arXiv:2108.11973.

[47] M. Fava et al., *Nonlinear sigma models for monitored dynamics of free fermions*, arXiv:2302.12820.

[48] B. Skinner, J. Ruhman and A. Nahum, *Measurement-Induced Phase Transitions in the Dynamics of Entanglement*, Phys. Rev. X 9 (2019) 031009 [arXiv:1808.05953] [INSPIRE].
[49] S. Antonini et al., *Holographic measurement and bulk teleportation*, JHEP 12 (2022) 124 [arXiv:2209.12903] [INSPIRE].

[50] I. Klich, *A note on the full counting statistics of paired fermions*, J. Stat. Mech. 11 (2014) P11006.

[51] Forrester, *Log-gases and random matrices*, Princeton University Press (2010) [DOI:10.1515/9781400835416].