Note on entropy dynamics in the Brownian SYK model

We study the time evolution of R\'enyi entropy in a system of two coupled Brownian SYK clusters evolving from an initial product state. The R\'enyi entropy of one cluster grows linearly and then saturates to the coarse grained entropy. This Page curve is obtained by two different methods, a path integral saddle point analysis and an operator dynamics analysis. Using the Brownian character of the dynamics, we derive a master equation which controls the operator dynamics and gives the Page curve for purity. Insight into the physics of this complicated master equation is provided by a complementary path integral method: replica diagonal and non-diagonal saddles are responsible for the linear growth and saturation of R\'enyi entropy, respectively.

Since the Brownian SYK model is amenable to both saddle point methods and circuit techniques, in this note we report a calculation of the Page curve in a model consisting of two coupled Brownian SYK clusters using both saddle point methods and operator dynamics. The quantity we are interested in is the Rény entropy of one cluster after tracing out the other. To formulate a path integral representation of the Rényi entropy, the initial state is taken to be a tensor product of thermofield double (TFD) states in each subsystem obtained by doubling the Hilbert space to left and right sides 1 . This is the Brownian version of the setup in [5]: the entanglement between left and right sides is maximal, whilst it is initially zero between the two subsystems. We find analytic solutions to the saddle point equations and show that replica diagonal and non-diagonal solutions are responsible for early time growth and the later time saturation of the Rényi entropy, respectively.
On the other hand, the density matrix dynamics can be directly analyzed using an operator dynamics approach. For simplicity we take a tensor product of Kourkoulou-Maldacena states [45] in each subsystem as initial state. Again, there is no entanglement between the two subsystems initially. Then the system is evolved under the full Hamiltonian, and the Page curve is obtained from a corresponding master equation. We compare the results from these two methods, and find excellent agreement even thought the two approaches use slight different initial states. Complementing the saddle point method, the master equation knows about the microstate from the perspective of symmetry, i.e., the Fermi parity in our case. The Hilbert space factorizes into different Fermi parity sectors, leading to an order one correction to the coarse grained entropy. Note that while the operator dynamics approach gives access only to the second Rényi entropy, in the saddle point method, we are able to get solutions for Rényi entropy for arbitrary Rényi index n.
The rest of the paper is organized as follows. In Section 2 we discuss a saddle point analysis of the path integral representation of Rényi entropy. Both replica diagonal and replica non-diagonal solutions are obtained analytically and checked numerically. The Page curve is obtained using these solutions, with the replica diagonal solution being responsible for the linear growth and the replica non-diagonal solution leading to the saturation to the coarse grained entropy. In section 3 we study the Page curve using a different approach, using the operator dynamics of the Brownian SYK model. We derive the master equation governing the operator size distribution function. The initial linear growth and late time saturation can be obtained analytically from the master equation, which shows exact agreement with the saddle point analysis.
2 Rényi entropy dynamics from saddle points

Coupled Brownian SYK clusters
The time-dependent Hamiltonian of two coupled Brownian SYK clusters labelled by a = 1, 2 is given by where A = j 1 ...j |A| denotes an ascending list of length |A|, q is an even integer, and [ψ a ] A ≡ i |A|/2 ψ j 1 ,a ψ j 2 ,a ...ψ j |A| ,a , is a short-hand notation for an |A|-body interaction. The ψ j,a , j = 1, ..., N a , are Majorana fermions in subsystem a, and satisfy {ψ j,a , ψ j ,a } = δ jj δ aa . The summations in Eq. (2.1) are over all possible lists with the indicated number of fermions. J a A (t) and V A,B (t) are Brownian random interactions within and between the two subsystems, respectively. The interaction strength is drawn from a Gaussian distribution distribution with mean zero and variance given by 3) The over line denotes an average over the Gaussian distribution of couplings. The interaction strength has dimension one (the dimension of energy), so J and V also have dimension one, while the δ-function makes up another dimension one and also indicates that the couplings are Brownian variables, uncorrelated in time. Regarding the prefactor, the dependence on N a is chosen to facilitate the large-N limit and the dependence on q is chosen to facilitate the large-q expansion in Appendix D. In general, the coupling between two subsystem does not have to be the same q-body interaction as the interaction within each subsystems, but we make such a choice for simplicity.

Setup
To investigate the entropy dynamics, we consider a similar setup to Ref. [5]: starting from the tensor product of two thermofield double (TFD) states in each of the subsystems (a = 1, 2), we focus on the Rényi entropy of subsystem a = 1 by tracing out subsystem a = 2. Because Brownian random interactions do not conserve energy, we simply consider an infinite temperature TFD state, which is a maximally entangled state. To prepare such a state, we double the Hilbert space by introducting left (L) and right (R) copies of the fermions, ψ j,a,L and ψ j,a,R , for both subsystems a = 1, 2. Then the maximally entangled state and the initial density matrix are given by (ψ j,a,L + iψ j,a,R )|∞ = 0, ∀a = 1, 2, ∀j = 1, ..., N, ρ 0 = |∞ ∞|.
Consider a time evolution generated by the sum of left and right Hamiltonians. The random couplings are identical between the two sides, up to an overall coefficient, with H L (t) = H(t; ψ j,a,L ) and H R (t) = (−1) q/2 H(t; ψ j,a,R ). This choice implies that H R |∞ = H L |∞ . Hence, the reduced density matrix ρ 1 of the subsystem a = 1 (including both L and R pieces) at time t is where T denotes time ordering, and Tr a denotes the trace over subsystem a. The n-th Rényi entropy is e −(n−1)Sn = Tr 1 [ρ 1 (t) n ]. This joint left-right evolution is equivalent to a single sided evolution for twice the time, i.e., U (t) = T e −i 2t 0 dt H(t ) .
A path integral representation of the trace of the n-th power of the reduced density matrix is obtained via a standard replica trick using n copies of the system and twist fields to implement modified boundary conditions on subsystem 1. This formulation gives Tr , where Z (n) , n ≥ 2 is the replicated partition function with twist operators inserted, and Z (1) ≡ Z is the partition function of a single replica. The replicated partition function Z (n) = [Dψ]e −I can be implemented in a Keldysh contour with two twist operators at t = 0 and t = T , respectively [5]. The insertion of twist operators in the contour is shown in Fig. 1. The effective action for the path integral of the replicated systems is where s = ± stands for the forward and backward contour, α, β = 1, 2, ..., n are the replica indices, and g αβ + = δ αβ , g αβ − = δ α+1,β ≡ αβ is due to the twist operator.  Figure 1. The schematic plot of the Keldysh contour and twist operator. The blue and red contours represent subsystem a = 1 and a = 2, respectively. The black dots at two ends of contour represent twist operators. The red contour is the α replica while in s = − of the blue contour it changes to α − 1 due to the twist operator.
At this point, the Brownian random interactions are integrated out. Strictly speaking, what we calculate is the logarithm of the average of the replicated partition function, i.e., 1 1−n log Tr 1 [ρ n 1 ], instead of the average Rényi entropy 1 1−n log Tr 1 [ρ n 1 ] which involves an additional replica trick to obtain the averaged logarithm. Nevertheless, due to the large-N structure of the Brownian models, we expect the circuit-to-circuit fluctuation are suppressed [29,46] so that both quantities agree with each other at large N . After integrating out the Brownian variables and introducing the bilocal fields G and Σ, we arrive at the following effective action, where t 12 ≡ t 1 − t 2 , and c ++ = c −− = −1, c +− = c −+ = 1 is due to the Keldysh evolution. The summation over the replica indices and the contour indices is implicit.σ z denotes the Pauli matrix acting on contour space.

Saddle point solutions
For simplicity, we assume the two clusters have equal numbers of Majorana fermions N 1 = N 2 = N . To look for a replica diagonal solution, we can start by looking for a solution when the inter-cluster coupling V = 0. In this case, the problem reduces to n independent replicas. Moreover, because of the Brownian nature of the problem, the self-energy is local in time (2.13). Starting from the Green's functions ansatz , (2.14) one uses (2.13) to get whereΣ a (ω) = dt 12Σa (t 12 )e iωt 12 . Here the limit J T 1 is implicit, as we are interested in the long-time behaviors of Rényi entropy, so the Fourier transform becomes an integral. We also solve the Schwinger-Dyson equation (2.12, 2.13) numerically in Appendix B for finite T , and find excellent agreement with the analytic solution we give in the following.
Plugging the self-energy into (2.12) gives . (2.17) Comparing it to the ansatz, we find f (t 12 ) = e − J q |t 12 | , so the replica diagonal solution is Now we claim that the above function (2.18) is still a solution to the Schwinger-Dyson equation with finite coupling V > 0. The reason is two-fold: (1) for the replica-diagonal solution, the V term in (2.13) is only non-vanishing on intra Keldysh contour s = s term due to the twist operator, and (2) the self-energy only depends on the Green's function at t 12 = 0 which is vanishing on the intra Keldysh contour s = s term. So if one plugs (2.18) into the Schwinger-Dyson equation, the V term in (2.13) vanishes and the it solves the equation.
Besides the replica diagonal solution, the twist operator induces new replica non-diagonal solutions that are the analog of the wormhole solutions found in [5]. We assume the subsystem a = 2 still hosts the replica diagonal solution, G αβ 2 , Σ αβ 2 ∝ δ αβ , and find that the subsystem a = 1 supports a replica non-diagonal solution. To get a replica non-diagonal solution, the nontrivial part must come from the twist operator in (2.13), corresponding to the inter Keldysh contour correlation function that crosses the twist operator as shown in Fig. 1. The self-energy of subsystem a = 1 is given by For a diagonal solution Σ αβ 1,−+ ∝ δ αβ , the second term in (2.19) vanishes, and the equation reduces to the diagonal solution (2.18) as we have shown. However, the second term in (2.19) also suggests a replica non-diagonal solution, Σ αβ 1,−+ ∝ αβ .
This leads us to consider a replica non-diagonal ansatz for subsystem a = 1 (and a replica diagonal ansatz for subsystem a = 2), , , (2.20) where the replica indices are implicit, and˜ αβ ≡ sgn(α − β)δ α+1,β such that˜ T = 1. The sign is a choice of ordering between different replicas and is also justified from numerical solutions. A detailed derivation is given in Appendix A, where one finds that f 1 and f 2 are q |t 12 | . The replica non-diagonal solutions read Thus, we have found both replica diagonal and non-diagonal solutions. We check these solutions by numerically iterating the Schwinger-Dyson equation (2.12, 2.13). In the numerical calculations, we focus on n = 2, 3. As shown in Fig. 2, the agreement between the analytics and the numerical solutions is quite good. Note the small contributions near the twist operators in the replica non-diagonal case. The analytical result above corresponds to the limit of large T where these boundary contributions can be neglected, but they do give important contributions to the on-shell action as we discuss below.

Page curve from saddle points
We now show that the replica diagonal and non-diagonal solutions lead to the linear increase and the saturation of the Rényi entropy, respectively. We first give analytic results for Rényi entropy from the two saddle point solutions, and then numerically evaluate the onshell action to verify our analytic results.
Recall that the replica diagonal solution ( solution of replicated action without inserting the twist operator because the twist operator effect is proportional to V. As a result, the first line in (2.11) counts the total Hilbert space dimension. For each replica system, we have two subsystems a = 1, 2 each hosting N Majorana fermions (remember we have doubled the Hilbert space to prepare the maximally entangled initial state within each subsystems), thus the total Hilbert dimension is 2 nN . In the second line of (2.11), for replica diagonal solution the inter Kelysh contour component s = s is zero. The intra Keldysh contour component s = s leads to a linear increase of Renyi-n entropy 2 , i.e., where the first term in S (1) is cancel by the denominator which is the Hilbert space of n replicated systems Z n = 2 nN .
For the replica non-diagonal solution, we first show that it leads to a time-independent on-shell action 3 . The on-shell action is a function of J T and VT , I = I[J T, VT ], so its time derivative reads Plugging the replica non-diagonal solution (2.20) into the time derivative, we have where the vanishing prefactor is coming from summing over c ss . The difference between replica diagonal and non-diagonal solutions originates from the twist operator, and the nondiagonal solution has a nontrivial contribution from the inter Keldysh contour component, so the time derivative vanishes. Because the bulk of the on-shell action vanishes, its value is determined by boundary effects near the twist operators.
Explicitly, the onshell action for the replica non-diagonal solution is given by where the second line vanishes for the same reason as in (2.26), and we use (2.12) to get the third line. Hence, the on-shell action is determined by the Pfaffian of the Green's function.
It is convenient to approach the calculation of the Pfaffian by recalling that the a = 2 subsystem still hosts the replica diagonal solution, which is also a solution when there is no twist operator in (2.11). Consider the action without twist operators, where the replica indices are trivial because there are no twist operators. Similar to the previous calculation, one can show that the large-N solution of (2.30) for both subsystems is the same as the diagonal solution in subsystem a = 2 given in (2.20). We copy it here using the same symbol for convenience, The onshell action of (2.30) is simply − which is actually the logarithm of Hilbert space dimension, i.e., e −I (0) = Z n . Using this action as denominator, the Rényi entropy from replica non-diagonal solution can be written as whereĜ 1 andĜ 2 are given by (2.21). As discussed, (2.21) neglects the effect of twist operators near the boundary, which must be included to get the correct answer. Thus we calculate the Pfaffian using the full numerical solution where the effect of twist operators is automatically included. We evaluate the Pfaffian for n = 2, 3, ..., 10 in Appendix C and find that We expect that this result holds true for any integer n because it gives the correct answer for Brownian evolutions where the Rényi entropy is maximized at late time.
From the on-shell actions of the replica diagonal solution (2.22) and replica non-diagonal solution (2.32), we find that the n-th Rényi entropy at time T is which is the Page curve in the coupled Brownian SYK models and T * = n−1 n q 2 V log 2+O(1/N ). To convert the result to entropy per Majorana, we notice that the number of Majorana fermion in the doubled Hilbert space of subsystem a = 1 is 2N , so we have where t * = T * /2 is the Page time. It is interesting to note the Page time increases for increasing n, but remains finite for n → ∞. We also evaluate the Rényi entropy of the replica diagonal and non-diagonal solution numerically. The results for n = 2, 3 are shown in Fig. 3.

Finite time effects
So far we showed that the replica diagonal solution gives rise to linear entropy growth while the replica non-diagonal solution gives rise to entropy saturation at long times. When the nondiagonal saddle dominates, there were small contributions localized near the twist operator that lead to important effects. Here we discuss the effects of these contributions on the timescale for entropy saturation.
In fact, there are important times in the problem. The first is the time at which the two saddles exchange dominance. We refer to this as the Page time and note that, at large N , it is independent of N . The second is the time for the entropy to reach within a few bits of its maximal value. We refer to this time as the strong scrambling time. As we now show, this time scales like log N at large N .
The key point is that in the replica non-diagonal saddle, boundary effects contribute a term in the on-shell action of the form where b n grows no faster than polynomial in T , and λ is a local interaction scale, which in our case is given by λ = 2(J + V)/q at large q limit. After the Page time, the Renyi entropy is thus Hence, as N goes to infinity at fixed T , the entropy differs from its saturation value by an amount extensive in N . On the other hand, by taking T ∼ 1 λ log N , the correction term can be made order unity instead of order N . At these times, the entropy is therefore within a few bits of its saturation value. This form of the on-shell action follows from the exponential decay of correlations. In the limit of large time, the twist operator contribution is local and cannot depend on the temporal extent of the system. Since correlations are exponentially decaying in time, it follows that finite time effects must vanish exponentially fast, up to a polynomial prefactor. A detailed analysis of this physics is possible in the large q limit as discussed in Appendix D.

Purity from operator dynamics
Here we study the entropy dynamics from the perspective of operator dynamics, focusing on the purity, e −S 2 . This approach is complementary to the path integral approach, including offering easier access to finite N corrections.

Operator dynamics of the coupled Brownian SYK models
As before, we consider a Hilbert space made up of N 1 +N 2 Majorana fermions. The Majorana fermions satisfy {ψ i,a , ψ j,b } = δ ab δ ij , and form an orthonormal basis where A = j 1 ...j |A| is an ascending list of length |A| and [x] is the largest integer less than or equal to x. In particular, we use Γ A,0 = Γ A,∅ and Γ 0,A = Γ ∅,A to denote the basis locating solely in subsystem a = 1, 2.
We again take N 1 = N 2 for simplicity, so the Hilbert space dimension is 2 N . The inner product of operators is Tr We are interested in finding the master equation governing the time development of this probability distribution.
In terms of the basis, the Hamiltonian (2.1) can be rewritten as where the superscript Γ (q) implies the length of basis is q. The short-hand summation denotes the same summation as in (2.1). Due to the Brownian nature of the interactions, we can view the Hamiltonian as a random circuit. At each time step, the evolution operator is generated by U (dt) = e −iH(t)dt . When we discretize the time interval by a tiny time step dt T , it is easy to check the average over Brownian variables takes the following rule In the following, the over line denoting the disorder average is omitted for notational simplicity.
The evolution operator can be approximated to the first-order in dt by where C m n ≡ In the second line we have perform the disorder average. Assuming the operator O is Hermitian, the distribution at time t + dt via the definition (3.3) is given by where the summation over C, D is the same as in (3.4). We have performed disorder average in the second line.
The master equation of p m,m (t) can be derived straightforwardly. We leave the detailed derivation in Appendix E, and the result is where the first two lines is the out-going rate and the rest is the in-coming rate. It is straightforward but tedious to show that the following distribution is a stationary solution to the master equation which means the probability in any basis Γ A,B is the same, i.e., |c A,B | 2 = 2 −2N . This is consistent with the expectation of approaching an infinite temperature state in the Brownian evolution, where no any specific basis is preferred.
It is instructive also to check the symmetry of the master equation. First, after the disorder average the model has an SO(N ) × SO(N ) symmetry, and this is the reason that the master equation can be reduced p A,B → p m,m depending only on the length of the basis.
Second, depending on the parity of q/2, i.e., the interactions between two Brownian SYK models, the model has Z f 2 ×Z f 2 for even q/2, i.e., the Fermi parity is conserved separately in two subsystems, or Z f 2 for odd q/2, i.e., only the total Fermi parity is conserved. This leads to the result that p m,m couples only to p m+q/2−2k,m +q/2−2k in (3.14). If q/2 is even, then the master equations for distribution p m,m , (m, m ) ∈ (even, even), (even, odd), (odd, even), (odd, odd) decouple. If q/2 is odd, then the master equations for distribution p m,m , m + m ∈ even, odd decouple.

Purity evolution of a pure state
Now we relate the purity evolution to the operator dynamics. The discussion in the following is general for any initial density matrix ρ (in) in the Hilbert space span by 2N Majorana operators. The density matrix evolves so the reduced density matrix is The evolution of the operator wavefunction c A,B is captured by the master equation (3.14), so the purity dynamics is also dictated by the master equation. We should notice that the density matrix is not normalized with respect to 2 −N Tr[O † O] = 1. For simplicity, let us consider pure initial state so that Tr[ρ 2 ] = 1. We can look at normalized operator, So the purity is which is 2 N/2 times the probability of finding the operator in subsystem a = 1.

Setup
We consider the following state (not be confused with the TFD state considered in the previous section), |∞ , that is similar to the Kourkoulou-Maldacena state [45] (ψ 2j−1,a + iψ 2j,a )|∞ = 0, ∞|(ψ 2j−1,a − iψ 2j,a ) = 0, ∀a = 1, 2, ∀j = 1, ..., N. (3.26) So the initial density matrix is where in the second line we used the basis defined in (3.1). The probability distribution of the normalized density matrix operator at time zero is which implies the purity is one (equivalently, the second Rényi entropy is zero), This is consistent with the fact that the initial state is a product state between the two subsystems.

Page curve from the master equation
The system is prepared in a pure state with the initial distribution given by (3.29). Though it is not easy to solve the master equation exactly, we can obtain the final probability distribution. Considering the symmetry of the master equation, the final probability distribution depends on the parity of q/2. If q/2 is even, the final distribution is The deficit of log 4 is due to Z f 2 × Z f 2 symmetry, since the Hilbert space has four decoupled sectors. If the initial state (3.27) is prepared in the (even, even) sector, its maximal entropy is given by N 2 log 2 − log 4.
On the other hand, if q/2 is odd, the final distribution is The shortage of log 2 is due to the Z f 2 symmetry. The Hilbert space has two decoupled sectors, leading to a deficit of − log 2.
Thus we expect that under the Hamiltonian (3.4) evolution, the second Rényi entropy will increase from zero to almost the largest value N 2 log 2. We can get the increase rate, which is the outgoing rate from P (t) = 2 N/2 m p m,0 (t), i.e., the second line in (3.14). The initial outgoing rate is where in the second line we take the large-N limit with s = m/N fixed, and use the Gaussian distribution to approximate the Binomial distribution, i.e., Thus at t 1/V, the second Rényi entropy grows linearly, At late time, as we know that the second Rényi entropy will saturate at N 2 log 2, we have thus the following Page curve of second Rényi entropy per Majorana fermion, where t * = q 2 4V log 2 is the Page time. Even though we start from a different initial density matrix, we can still compare the result from the saddle point solution. We find that two results match exactly for second Rényi entropy per Majorana (2.36). This is the useful quantity since the Hilbert dimensions are different for the two cases.
We also numerically solve the master equation starting from (3.29) for q/2 = 2. The purity evolution as a function of time is plotted as the solid line in Fig. 4(a), explicitly showing the Page curve in the Brownian evolution. The dotted line and dash line are given by (3.41). The deviation from the linear increasing dotted line at early time is due to finite N effect. For the dashed line, we have included the finite N correction from Z f 2 × Z f 2 symmetry. Fig. 4(b) shows that after the Page time, the approach to the thermal value is given by an exponential function. This is consistent with a log N scrambling time as discussed in Sec. 2.5, where it is explained by the effect of twist operators.

Conclusion
We studied the Rényi entropy dynamics of coupled Brownian SYK clusters using both path integral and operator dynamics methods. While the Page curve has been observed in random circuit models before, we showed how the replica diagonal and non-diagonal saddle points give rise to the entanglement behavior. This structure is very similar to the replica wormhole scenario obtained in holographic calculations. We also discussed the scrambling time and Page time in Section 2.5 and Appendix D, where the crucial effect of the twist operator was discussed. One interesting future direction is to calculate the entanglement entropy directly, and consequently to reveal the role of entanglement islands in more generic quantum mechanical systems.

A Replica non-diagonal solution
According to the ansatz (2.20) and the Schwinger-Dyson equation (2.13), the self-energy iŝ Then (2.12) leads tô Comparing this to the ansatz (2.20), we find the replica non-diagonal solution

B The solution at finite time
We solve the Schwinger-Dyson equation numerically to compare it with the analytic solutions (2.18, 2.21). Starting from a noninteracting solution as an input, we iterate the Schwinger-Dyson equation (2.12, 2.13) until the result converges. This iteration was used in [47] and also in [5].
We show the results in

C Numerical calculation of the saddle point solution and the onshell action
For numerical convenience, here we adopt a different convention for the labeling of fields. We put both Keldysh contour indices and replica indices into the time argument 0 < s < 2nT (a similar convention is used in [48]): The forward contour for the α-th replica is s ∈ ((2α − 2)T, (2α − 1)T ), and the backward contour for the α-th replica is s ∈ ((2α − 1)T, 2αT ). We also introduce a sign factor to capture the forward and backward contour, In this convention we also adapt the action such that the interaction between two clusters is local at s. As an illustration, the contour convention for n = 2 is shown in Fig. 6. In this case, as seen from the figure, a replica diagonal solution in subsystem a = 1 will have nonvanishing correlation between s ∈ (0, 2t) and s ∈ (6t, 8t) and between s ∈ (2t, 4t) and s ∈ (4t, 6t). On the other hand, a nonvanishing correlation between s ∈ (0, 2t) and s ∈ (2t, 4t) and between s ∈ (4t, 6t) and s ∈ (6t, 8t) for subsystem a = 1 implies a replica non-diagonal solution. We will assume a = 2 has a diagonal solution in the following.
In terms of this convention, the action is where A = j 1 ...j |A| denotes an ascending list of length |A|, and [ψ a ] A ≡ i |A|/2 ψ j 1 ,a ψ j 2 ,a ...ψ j |A| ,a is a short-hand notation for |A|-body interaction. The summation is over all possible such lists from N a Majorana fermions. In general, the interaction strength is random variable with possible dependence on time s. The distributions of the interactions are defined by vanishing means and the following variances, where The effective action after averaging over random variables and introducing bilocal fields, i.e. the Green's function G a (s 1 , s 2 ) = 1 Na Na j=1 ψ j,a (s)ψ j,a (s ) and the self-energy Σ a (s, s ), reads The Schwinger-Dyson equation follows from the effective action is given bŷ For simplicity, we consider N 1 = N 2 = N . We numerically solve the Schwinger-Dyson equation (C.10, C.11) for n = 2, 3 and look for replica non-diagonal solution. In doing so, we use G 0,a (s, s ) = 1 2 sgn(s − s ) for times located at the same close time path. An illustration of the close time bath for n = 2 is given in Fig. 6. To get the replica non-diagonal solution, we start from an initial ansatz with small but non-zero non-diagonal correlations for subsystem a = 1 and a diagonal initial ansatz for subsystem a = 2. The results for n = 2 are shown in Fig. 7. As we discuss in above, Fig. 7(b) is a replica non-diagonal solution. It is intuitive to note from the figures that the only difference between the replica diagonal solution and the replica non-diagonal solution is those nonvanishing correlations at {4t < s 1 < 8t} ∩ {0 < s 2 < 4t} and {0 < t 1 < 4t} ∩ {4t < t 2 < 8t} sourced by the twist operators located at s = 0, 2t, 6t, 8t. We also get the results for n = 3 Rényi entropy, which are shown in Fig. 8(b), and there are six twist operators.
We also numerically calculate the Pfaffian in the calculation of onshell action and check the validity of (2.33). The result is plotted in Fig. 9 where we calculate log Pf(G −1 1 G 2 ) for the

D Non-diagonal solutions and twist operators
In this section, we discuss the effect of the twist operators. We will focus on the replica nondiagonal solution for the second Rény entropy n = 2 for simplicity, while the generalization to other Rényi entropy is straightforward. The first equation (C.10) couples functions nonlocally in time domain, while the second equation (C.11) is local. Using the large-q ansatz, G a = g a0 (1 + ga q ), the Schwinger-Dyson equation becomes, One advantage of the large-q equation of motion is that it becomes local in time variables. For simplicity, we will assume N 1 = N 2 = N . We also assume g 20 = 1 2 sgn(s 1 − s 2 ) when s 1 , s 2 are located at the same close time path and zero otherwise, while g 10 = 1 2 sgn(s 1 − s 2 ) to look for the non-diagonal solutions.
We can solve it in the regime 0 < s 1 < T and 0 < s 2 < T . For Brownian case, the large-q saddle point equation reads, The equation can be solved by realizing it is ∂ s 1 ∂ s 2 g a (s 1 , s 2 ) = 0 when s 1 = s 2 , and the δ function leads to a jump in the first derivative at s = s 1 = s 2 that can be solved easily. Supplementing with the boundary condition g a (s, s) = 0 and g a (s 1 , s 2 ) = g a (s 2 , s 1 ), the solutions are g 1 (s 1 , s 2 ) = g 2 (s 1 , s 2 ) = −(J + V)|s 1 − s 2 |, {0 < s 1 < T, 0 < s 2 < T }. (D.5) We can extend such calculations to the regime 0 < s 1 < 2T and 0 < s 2 < 2T , This is consistent with the analytic solution (2.21) and the numeric solution shown in Fig. (7). More generally, the solutions in (α − 1)2T < s 1 , s 2 < 2αT , α = 1, ..., n will be the same.
To investigate the nonvanishing correlation induced by the twist operator near the boundary, we first focus on the regime 3T < s 1 < 4T and 0 < s 2 < T . In this regime, there are two twist operators, where T 1 locates at (4T, 0) and T 2 locates (3T, T ) as also indicated by the nonvanishing correlations in Fig. 7. So the boundary conditions are ψ(T ) = ψ(3T ) and ψ(0) = −ψ(4T ). The minus sign is because of the Fermi operator. To simplify the notation, we shift s 1 → s 1 + 3T , so the regime is 0 < s 1 , s 2 < T . After the redefinition, the large-q equation of motion in this regime is (D.7) The absence of the V term is because the replica diagonal solution of subsystem a = 2 vanishes in this regime.
The solution is exponentially suppressed at large s, and at the large time, i.e., J T 1, the two twist operators are separated by a large distance T . So we can further assume in this regime the induced solution is separable as follows, where G T i denotes the induced solution by twist operator T i . At large-q limit, G T i ≈ (1 + g T i q ), and g T i satisfies large-q equation of motion (D.7). But they satisfy different boundary conditions because the two twist operators locate at different places, i.e., Let us first look at g T 1 . Owing to the delta function in the right-hand-side of (D.7), the solution is not differentiable at s 1 + s 2 = T , so we assume Taking the boundary conditions into consideration, the solution has the form, and f T 1 satisfies It is not hard to solve above differential equation, which leads to the solution, (D.14) And consequently, one can get the correlation function G T 1 ≈ 1 2 e g T 1 /q induced by the twist operator T 1 , To simplify the notation, we define the induced Green's function as 16) Intuitively, this function is exponentially suppressed away from s 1 = s 2 = 0 where the twist operator supposes to be located. So the solution at the regime is This is an approximate solution with accuracy O(e −J T ). The solutions in other regimes can be obtained in the same way, so we do not have to detail the calculation. Now we can discuss the effect of these twist operators to the onshell action. We mainly discuss the second Rényi entropy, but the results can be extended to n-th Rényi entropy by small modifications. Taking into account the twist operators, the replica non-diagonal solution is G 1 = G 2 + G T 1 + G T 2 , where G 2 denotes the diagonal solution in subsystem a = 2 (D.6), and G T 1 and G T 2 are induced solutions by the twist operators T 1 and T 2 , respectively. Notice G 2 and G T i have different domain of support. The onshell action is where in the second line, we expand the Tr log term and keep the lowest-order coupling between two twist operators G T 1 and G T 2 , because the factorized part contributes to the course grained entropy −N log 2 [48], and we neglect other subleading terms in the large-q limit. So including the parts from coupled induced Green's function, the contribution from the replica non-diagonal solution is (D.20) In getting above results, we neglect the second term in (D.16) which will not change the essential exponential factor. Then the second Rényi entropy from two saddle points reads After the Page time when the replica non-diagonal saddle point dominates, the Rényi entropy is actually not independent of time. The exponentially small overlaps between two twist operators mean that it takes times proportional to log N to fully scramble the information [27,49].
The large-q analysis of the twist operator can be extended to the regular SYK model. Here we calculate it at the infinite temperature for an illustration. The solutions at diagonal part is simple, yielding the solution where J 0 = √ J 2 + V 2 . So let us focus again on the regime 3T < s 1 < 4T and 0 < s 2 < T . Redefining s 1 → s 1 + 3T , the regime is 0 < s 1 , s 2 < T . The equation of motion now reads ∂ s 1 ∂ s 2 g 1 (s 1 , s 2 ) = −2J 2 e g 1 (s 1 ,s 2 ) . (D.23) The absence of the V term is because the replica diagonal solution of subsystem a = 2 vanishes in this regime. A general solution to above Liouville equation is g 1 (s 1 , s 2 ) = log h 1 (s 1 )h 2 (s 2 ) J 2 (h 1 (s 1 )−h 2 (s 2 )) . We expect the solution is exponentially suppressed at large s, and at the large time, i.e., J T 1, the two twist operators are separated by a large distance T . So we can further assume in this regime the induced solution is separable as follows, G 1 (s 1 , s 2 ) = G T 1 (s 1 , s 2 ) + G T 2 (s 1 , s 2 ), (D. 24) where G T i denotes the induced solution by twist operator T i . At large-q limit, G T i ≈ (1+ g T i q ), and g T i satisfies the Liouville equation. But they satisfy different boundary conditions because the two twist operators locate at different places, i.e., g T 1 (s 1 , 0) = 2 log 1 cosh J 0 (T − s 1 ) , g T 1 (T, s 2 ) = 2 log 1 cosh J 0 s 2 , (D. 25) g T 2 (s 1 , T ) = 2 log 1 cosh J 0 s 1 , g T 2 (0, s 2 ) = 2 log 1 cosh J 0 (T − s 2 ) . (D.26) After we take into account the boundary conditions, it is straightforward to get the following solutions, It will interesting to explore the effect of these twist operators in more details, which we leave as a future work.